Example: Heat flow Notes Let us first solve the heat equation analytically before discussing the numerical method to solve it. Example: Heat equation in one spatial dimension. ∂T ∂2T =κ 2, ∂t ∂x where κ is the thermal diffusivity. Open boundaries: T (x, t) defined on −∞ < x < +∞ and t ≥ 0. Also, require that T (x, t) → 0 as x ± ∞ Initial value problem: T (x, t = 0) = Φ(x) We will use the Fourier transform in order to solve this problem.

Solving Heat flow with an integral transform We have two independent variables x, t. Apply Fourier transform to spatial variables x at constant t: Z +∞ 1 F[T (x, t)] = √ T (x, t)e −ikx dx 2π −∞

Notes

What is the Fourier transform of ∂T /∂x ? Z +∞ ∂T −ikx 1 F[Tx ] = √ e dx 2π −∞ ∂x Integrating by parts, we obtain  −ikx ∞ Z +∞ e ik F[Tx ] = √ T (x, t) +√ T (x, t)e −ikx dx 2π 2π −∞ −∞ Since T (x, t) vanishes as x → ±∞, the first term vanishes and we have F[Tx ] = ikF[T ]

Solving Heat flow with an integral transform Notes The Fourier transform of ∂T /∂t is simply Z +∞ 1 ∂T −ikx F[Tt ] = √ e dx 2π −∞ ∂t ∂ = F[T (x, t)] ∂t using Leibniz rule. Similarly, ∂2 F[T (x, t)] ∂t 2 2 F[Txx ] = (ik) F[T (x, t)] = −k 2 F[T (x, t)] F[Ttt ] =

Solving Heat flow with an integral transform Notes

Now apply Fourier transform to the heat equation (at constant t): F[Tt ] = F[κTxx ] ∂ F[T ] = −k 2 κF[T ] ∂t Note: Fourier transform method for solving this problem is appropriate since T (x, t) vanishes at x → ±∞ according to the BC. Let us denote τ (k, t) = F[T (x, t)]

(1) (2)

Solving Heat flow with an integral transform Notes ∂τ (k, t) = −k 2 κτ (k, t) (3) ∂t Equation (3) is now a simple ordinary differential equation. Heat equation is much easier to solve in the Fourier domain! The solution is τ (k, t) = e −k

2 κt

τ (k, 0)

(4)

Still need to transform the initial condition T (x, 0) = Φ(x): F[Φ(x)] = F[T (x, t = 0)] = τ (k, 0)

(5)

Solving Heat flow with an integral transform Notes Combining eqs. (4) and (5) τ (k, t) = e −k

2 κt

(6)

F[Φ(x)]

In order to obtain solution in real space, apply the inverse Fourier transform: T (x, t) = F−1 [τ (k, t)] = F−1 [e −k

2 κt

F[Φ(x)]]

(7)

However, can use the convolution theorem on the right hand side. Recall F[f ⊗ g ] = F[f ]F[g ]

Therefore, f ⊗ g = F−1 [F[f ]F[g ]]

(8)

Notes

Now apply this to our solution (eq. (7)): 2 Let F[f ] = e −k κt and F[g ] = F[Φ(x)]. It follows that g = F−1 [F[Φ(x)]] = Φ(x) by definition. Also, f

2

= F−1 [e −k κt ] 1 2 = √ e −x /(4κt) 2κt

(9) (10)

In the last step we used the known inverse Fourier transform of a Gaussian (see table of Fourier transforms - lecture 2): 1 2 2 F[e −αx ] = √ e −k /(4α) 2α

Solving Heat flow with an integral transform Notes Since T (x, t) = f ⊗ g

(11)

according to eqs. (8) and (7) we have T (x, t) = = =

1 2 e −x /(4κt) ⊗ Φ(x) 2κt Z +∞ 1 1 2 √ √ e −(x−ξ) /(4κt) Φ(ξ)dξ 2π −∞ 2κt Z +∞ 1 2 √ e −(x−ξ) /(4κt) Φ(ξ)dξ 4πκt −∞

√

(12) (13) (14)

This is the fundamental solution of the heat equation with open boundaries for an initial condition T (x, t = 0) = Φ(x).

Green’s function Notes The fundamental solution is the the convolution of the initial conditions with the Green’s function. Once the Green’s function of a linear PDE is known (for given B.C.’s), it can be solved for arbitrary initial conditions via a convolution integral (eq. (14)). The Green’s function is the solution of the PDE for a delta (impulse) function. Let the initial condition be Φ(x) = δ(x). Then Z +∞ 1 2 T (x, t) = √ e −(x−ξ) /(4κt) δ(ξ)dξ 4πκt −∞ 1 2 = √ e −x /(4κt) 4πκt Note: T (x, t) → δ(x) as t → 0.

Heat diffusion from delta function impulse Notes

The delta function at t = 0 becomes spread out as time goes on.

Finite difference methods Notes

In the following, let’s see how to solve PDE’s using finite difference methods. For simplicity, assume u to be a function of two independent variables x and t, only. The continuous function u is then discretized on a grid in the x, t plane where the points are separated by ∆x and ∆t, respectively. Notation: u(x, t) → u(j∆x, n∆t), where j and n are integers.

Approximating the derivatives In order to discretize the PDE we need to approximate the derivatives appearing in the PDE by finite differences. There are three possible approximations for the first order (partial) derivative: ∂u(j∆x, n∆t)/∂x The forward (Euler) difference: n − un uj+1 j

∆x The backward difference: n ujn − uj−1

∆x The centered difference: n − un uj+1 j−1

2∆x

Notes

Approximating the derivatives The same applies for the temporal partial derivative, except that j is kept constant. e.g. the forward difference is

Notes

ujn+1 − ujn ∂u(j∆x, n∆t) ≈ ∂t ∆t The finite differences follow from the Taylor expansion of u(x) (at constant t): u(x + ∆x) =

u(x − ∆x) =

1 u(x) + u 0 (x)∆x + u 00 (x)(∆x)2 2 1 + u 000 (x)(∆x)3 + O(∆x)4 6 1 u(x) − u 0 (x)∆x + u 00 (x)(∆x)2 2 1 − u 000 (x)(∆x)3 + O(∆x)4 6

(15) (16) (17) (18)

where in the last line ∆x → −∆x

Approximating the derivatives Notes

From the Taylor series expansion we can find the local errors: The forward (Euler) difference (eq.(17)): u(x + ∆x) − u(x) + O(∆x) ∆x The backward difference (eq.(18)): u(x) − u(x − ∆x) + O(∆x) ∆x The centered difference(subtract eq.(18) from eq. (17)): u(x + ∆x) − u(x − ∆x) + O(∆x)2 2∆x Note that the global (accumulated) error may be different from the local one.

Approximating the derivatives Notes What about the second derivative? By adding the two Taylor expansions (17) and (18) and solving for u 00 (x) we obtain u 00 (x) =

u(x + ∆x) − 2u(x) + u(x − ∆x) + O(∆x)2 (∆x)2

(19)

In index notation, this reads for a partial second derivative at constant t: n − 2u n + u n uj+1 ∂u(j∆x, n∆t) j j−1 ≈ 2 ∂x (∆x)2

(20)

This is the second centered difference. Similar expression for partial time derivative.

Discretizing the heat equation Notes Let us use a finite difference scheme to evaluate the heat equation numerically in one spatial dimension. Tt = κTxx

(21)

Again, assume open boundaries with initial condition: T (x, t = 0) = Φ(x)

Apply forward difference to Tt and second centered difference to Txx : n − 2T n + T n Tjn+1 − Tjn Tj+1 j j−1 =κ (22) ∆t (∆x)2 where Tjn = T (j∆x, n∆t).

Discretizing the heat equation - errors Notes There are two types of errors in this finite difference scheme:

Truncation errors: The local truncation error for the second centered difference is ∆x 2 and for the time derivative ∆t. These errors may accumulate and given an overall global truncation error that may be bigger than the local one. Round off error: This occurs in real computations, since the computer only retains a certain number of digits for floating point numbers. e.g. in C, float variables have 8 digit and double have 16 digit precision. These round off errors also accumulate. Always use double variables for finite difference schemes!

Discretizing the heat equation Notes Solving for

Tjn+1

we obtain

n n Tjn+1 = (1 − 2s)Tjn + s(Tj+1 + Tj−1 )

where s =

(23)

κ∆t . (∆x)2

This is an explicit difference scheme, since the values of the (n + 1)th time step are explicitly given in terms of the value at earlier times.

Example 1 Notes Let’s look at a simple example. For simplicity, set s = 1. Then the difference scheme simplifies to n n Tjn+1 = Tj+1 − Tjn + Tj−1

(24)

Tjn=0

Suppose Φ(xj ) is a simple impulse: = 1 for a particular point j in space and zero everywhere else. Can integrate this by “marching in time”:

Example 1 Notes Something has gone wrong ! Expect the impulse to diffuse away. According to the maximum principle of diffusion PDE’s, the maximum of T (x, t) occurs either at the initial conditions t = 0 or at the boundaries. In our example, this means that T (x, t) < 1 for t > 0. Instead the difference equations has given us an approximation where the central point is 19 and growing. It turns out that for this particular difference scheme ∆t and ∆x cannot be chosen arbitrarily.

Example 2 Notes

Let’s consider another example of the heat equation with the following boundary conditions: Tt = Txx

for 0 < x < π, t > 0

T =0

at x = 0, π

The initial condition is  T (x, 0) = Φ(x) =

x π−x

for 0 < x ≤ π/2 for π/2 < x < π

(25)

First, let’s solve this analytically using separation of variables.

Example 2 Notes

If the PDE is separable, then its solution can be expressed as T (x, t) = A(x)B(t) This allows us to turn the PDE into two ordinary differential equations (ODE’s). Substitution yields ∂T ∂t ∂B(t) A(x) ∂t

= κ

∂2T ∂x 2

= κB(t)

(26) ∂ 2 A(x) ∂x 2

(27)

Dividing by (A(t)B(x)κ) we obtain 1 ∂B(t) 1 ∂ 2 A(x) = −k 2 = κB(t) ∂t A(x) ∂x 2 where −k 2 is the constant of separation. LHS and RHS are equal for all x and t, which are independent of each other.

Example 2 Notes The two ODE’s are 1 ∂B(t) κB(t) ∂t 1 ∂ 2 A(x) A(x) ∂x 2

= −k 2

(28)

= −k 2

(29)

The solutions are A(x) = a cos(kx) + b sin(kx) B(t) = ce −κk

2t

(30) (31)

where a, b, c are constants of integration. The boundary conditions T (x = 0, t) = T (x = π, t) = 0, require a = 0 and k to be an integer.

Example 2 Notes Therefore, T (x, t) = bk sin(kx)e −κk

2t

(32)

where k is an integer. Fourier analysis allows us to obtain the solution for our specified initial condition: T (x, t) =

∞ X

bk sin(kx)e −κk

2t

(33)

k=1

where

( bk =

(−1)(k+1)/2 πk 2

0

for odd k for even k

(34)

Example 2 Notes This plot shows the initial temperature distribution Φ(x) and T (x, t) for t = 3π 2 /80.

As t → ∞, T → 0 everywhere.

Numerically solving example 2 Notes

Now solve numerically using our difference scheme n n Tjn+1 = (1 − 2s)Tjn + s(Tj+1 + Tj−1 )

(35)

κ∆t with s = (∆x) 2 = 5/11 = 0.45 and ∆x = π/20, where we have discretized space with J + 1 = 21 points. The discrete boundary and initial conditions are

u0n = uJn = 0 and uj0 = Φ(j∆x) where j = 0, 1, 2, . . . J.

Numerically solving example 2 Notes

Agrees very well with analytical solution.

Numerically solving example 2 Notes However, if we choose s = 5/9 = 0.55, the finite difference scheme yields wild oscillations similar to what we obtained for the impulse IC earlier.

Stability of finite difference scheme Notes In the following we will derive the stability criterion for this difference scheme: Let’s solve the difference scheme by separation of variables: Tjn = Aj Bn

(36)

Substitution in the difference scheme (eq. (35)) yields Aj Bn+1 = (1 − 2s)Aj Bn + s(Aj+1 Bn + Aj−1 Bn )

(37)

Dividing by Aj Bn we get Aj+1 + Aj−1 Bn+1 = ξ = 1 − 2s + s Bn Aj

(38)

where ξ is the constant of separation.

Stability of finite difference scheme Notes

Thus, Bn+1 =ξ Bn

(39)

Bn = ξ n B0

(40)

which has the solution where B0 is some constant. Now solve the spatial part: 1 − 2s + s

Aj+1 + Aj−1 =ξ Aj

(41)

We wish to investigate the eigenmodes of the solution. Therefore, assume a plane wave solution with wave number k (ignore boundary conditions for this analysis): Aj = e ikj∆x

(42)

Stability of finite difference scheme Notes Substituting yields e ik(j+1)∆x + e ik(j−1)∆x = ξ e ikj∆x ik∆x −ik∆x 1 − 2s + s(e +e ) = ξ

1 − 2s + s

1 − 2s + 2s cos(k∆x) = ξ 1 − 2s(1 − cos(k∆x)) = ξ   2 k∆x 1 − 4s sin = ξ 2 For the difference scheme to be stable, require |ξ(k)| ≤ 1 for all k since the amplitude grows as ∝

ξn

(43)

(eq. (40))

Stability of finite difference scheme Notes Therefore,

  2 k∆x 1 − 4s sin ≤1 2

for all k. This leads to the following stability criterion: 4s ≤ 2 2κ∆t ≤ 1 (∆x)2 Physically, this means that the maximum allowed time step ∆t is the diffusion time across a cell of size ∆x (up to a numerical 2 factor). The diffusion time across some length λ goes as tλ ∼ λκ

von Neumann stability analysis Notes This is the von Neumann stability analysis: The von Neumann analysis looks at the evolution of the eigenmodes of the solution by substituting ujn = ξ n e ikj∆x into the difference equation, where ξ(k) is the amplification factor For the difference scheme to be stable, require |ξ(k)| ≤ 1 + O(∆t) for all k Otherwise, eigenmodes grow exponentially in time.

Notes

Notes

Notes