Computational Physics - Department of Physics

(Axel Boer) #1

10.2 Diffusion equation 311


θ
∆x^2

(

ui− 1 ,j− 2 ui,j+ui+ 1 ,j

)

+

1 −θ
∆x^2

(

ui+ 1 ,j− 1 − 2 ui,j− 1 +ui− 1 ,j− 1

)

=

1

∆t

(

ui,j−ui,j− 1

)

, (10.9)

which forθ= 0 yields the forward formula for the first derivative and the explicit scheme,
whileθ= 1 yields the backward formula and the implicit scheme. These two schemes are
called the backward and forward Euler schemes, respectively. Forθ= 1 / 2 we obtain a new
scheme after its inventors, Crank and Nicolson. This schemeyields a truncation in time which
goes likeO(∆t^2 )and it is stable for all possible combinations of∆tand∆x.
To derive the Crank-Nicolson equation, we start with the forward Euler scheme and Taylor
expandu(x,t+∆t),u(x+∆x,t)andu(x−∆x,t)


u(x+∆x,t) =u(x,t)+∂u∂(xx,t)∆x+∂

(^2) u(x,t)
2 ∂x^2 ∆x
(^2) +O(∆x (^3) ), (10.10)
u(x−∆x,t) =u(x,t)−∂u∂(xx,t)∆x+∂
(^2) u(x,t)
2 ∂x^2 ∆x
(^2) +O(∆x (^3) ),
u(x,t+∆t) =u(x,t)+∂u(∂xt,t)∆t+O(∆t^2 ).
With these Taylor expansions the approximations for the derivatives takes the form
[∂u(x,t)
∂t


]

approx
=∂u∂(xt,t)+O(∆t), (10.11)

[∂ (^2) u(x,t)
∂x^2


]

approx

=∂

(^2) u(x,t)
∂x^2 +O(∆x
(^2) ).
It is easy to convince oneself that the backward Euler methodmust have the same truncation
errors as the forward Euler scheme.
For the Crank-Nicolson scheme we also need to Taylor expandu(x+∆x,t+∆t)andu(x−
∆x,t+∆t)aroundt′=t+∆t/ 2.
u(x+∆x,t+∆t) =u(x,t′)+∂u(x,t
′)
∂x ∆x+
∂u(x,t′)
∂t
∆t
2 +
∂^2 u(x,t′)
2 ∂x^2 ∆x
(^2) +∂^2 u(x,t′)
2 ∂t^2
∆t^2
4 +
∂^2 u(x,t′)
∂x∂t
∆t
2 ∆x+O(∆t
(^3) )
u(x−∆x,t+∆t) =u(x,t′)−∂u(x,t
′)
∂x ∆x+
∂u(x,t′)
∂t
∆t
2 +
∂^2 u(x,t′)
2 ∂x^2 ∆x
(^2) +∂^2 u(x,t′)
2 ∂t^2
∆t^2
4 −
∂^2 u(x,t′)
∂x∂t
∆t
2 ∆x+O(∆t
(^3) )
u(x+∆x,t) =u(x,t′)+∂u(x,t
′)
∂x ∆x−
∂u(x,t′)
∂t
∆t
2 +
∂^2 u(x,t′)
2 ∂x^2 ∆x
(^2) +∂^2 u(x,t′)
2 ∂t^2
∆t^2
4 −
∂^2 u(x,t′)
∂x∂t
∆t
2 ∆x+O(∆t
(^3) )
u(x−∆x,t) =u(x,t′)−∂u(x,t
′)
∂x ∆x−
∂u(x,t′)
∂t
∆t
2 +
∂^2 u(x,t′)
2 ∂x^2 ∆x
(^2) +∂^2 u(x,t′)
2 ∂t^2
∆t^2
4 +
∂^2 u(x,t′)
∂x∂t
∆t
2 ∆x+O(∆t
(^3) )
u(x,t+∆t) =u(x,t′)+∂u(x,t
′)
∂t
∆t
2 +
∂^2 u(x,t′)
2 ∂t^2 ∆t
(^2) +O(∆t (^3) )
u(x,t) =u(x,t′)−∂u(x,t
′)
∂t
∆t
2 +
∂^2 u(x,t′)
2 ∂t^2 ∆t
(^2) +O(∆t (^3) )
We now insert these expansions in the approximations for thederivatives to find
[∂u(x,t′)
∂t


]

approx
=∂u(x,t

′)
∂t +O(∆t

(^2) ), (10.12)
[∂ (^2) u(x,t′)
∂x^2


]

approx=

∂^2 u(x,t′)
∂x^2 +O(∆x

(^2) ).
Bringing all these equations together results in Eq. (10.9)withθ= 1 / 2.
The following table summarizes the three methods.
Using our previous definition ofα=∆t/∆x^2 we can rewrite Eq. (10.9) as

Free download pdf