306
9 Solving Partial Differential Equations
b) The Backward Euler, Forward Euler, and Crank-Nicolson methods can be given
a unified implementation. For a linear ODE u = au this formulation is known
as the θ rule:
u n+1 − u n
Δt
= (1 − θ)au
n
+ θau
n+1 .
For θ = 0 we recover the Forward Euler method, θ = 1 gives the Backward
Euler scheme, and θ = 1/2 corresponds to the Crank-Nicolson method. The
approximation error in the θ rule is proportional to Δt, except for θ = 1/2
where it is proportional to Δt 2 . For θ ≥ 1/2 the method is stable for all Δt.
Apply the θ rule to the ODE system for a one-dimensional diffusion equation.
Identify the linear system to be solved.
c) Implement the θ rule with aid of the Odespy package. The relevant object name
is ThetaRule:
solver = odespy.ThetaRule(rhs, f_is_linear=True, jac=K, theta=0.5)
d) Consider the physical application from Sect. 9.2.4. Run this case with the θ rule
and θ = 1/2 for the following values of Δt: 0.001, 0.01, 0.05. Report what you
see.
Filename: rod_ThetaRule.py.
Remarks Despite the fact that the Crank-Nicolson method, or the θ rule with
θ = 1/2, is theoretically more accurate than the Backward Euler and Forward
Euler schemes, it may exhibit non-physical oscillations as in the present example
if the solution is very steep. The oscillations are damped in time, and decreases
with decreasing Δt. To avoid oscillations one must have Δt at maximum twice the
stability limit of the Forward Euler method. This is one reason why the Backward
Euler method (or a 2-step backward scheme, see Exercise 9.3) are popular for
diffusion equations with abrupt initial conditions.
Exercise 9.6: Compute the Diffusion of a Gaussian Peak
Solve the following diffusion problem:
∂u
∂t
= β
∂ 2 u
∂x 2 ,
x ∈ (−1, 1), t ∈ (0, T ]
(9.34)
u(x, 0) =
1
√
2πσ
exp
−
x 2
2σ 2
,
x ∈ [−1, 1],
(9.35)
∂
∂x
u(−1, t) = 0,
t ∈ (0, T ],
(9.36)
∂
∂x
u(1, t) = 0,
t ∈ (0, T ] .
(9.37)
9 Solving Partial Differential Equations
b) The Backward Euler, Forward Euler, and Crank-Nicolson methods can be given
a unified implementation. For a linear ODE u = au this formulation is known
as the θ rule:
u n+1 − u n
Δt
= (1 − θ)au
n
+ θau
n+1 .
For θ = 0 we recover the Forward Euler method, θ = 1 gives the Backward
Euler scheme, and θ = 1/2 corresponds to the Crank-Nicolson method. The
approximation error in the θ rule is proportional to Δt, except for θ = 1/2
where it is proportional to Δt 2 . For θ ≥ 1/2 the method is stable for all Δt.
Apply the θ rule to the ODE system for a one-dimensional diffusion equation.
Identify the linear system to be solved.
c) Implement the θ rule with aid of the Odespy package. The relevant object name
is ThetaRule:
solver = odespy.ThetaRule(rhs, f_is_linear=True, jac=K, theta=0.5)
d) Consider the physical application from Sect. 9.2.4. Run this case with the θ rule
and θ = 1/2 for the following values of Δt: 0.001, 0.01, 0.05. Report what you
see.
Filename: rod_ThetaRule.py.
Remarks Despite the fact that the Crank-Nicolson method, or the θ rule with
θ = 1/2, is theoretically more accurate than the Backward Euler and Forward
Euler schemes, it may exhibit non-physical oscillations as in the present example
if the solution is very steep. The oscillations are damped in time, and decreases
with decreasing Δt. To avoid oscillations one must have Δt at maximum twice the
stability limit of the Forward Euler method. This is one reason why the Backward
Euler method (or a 2-step backward scheme, see Exercise 9.3) are popular for
diffusion equations with abrupt initial conditions.
Exercise 9.6: Compute the Diffusion of a Gaussian Peak
Solve the following diffusion problem:
∂u
∂t
= β
∂ 2 u
∂x 2 ,
x ∈ (−1, 1), t ∈ (0, T ]
(9.34)
u(x, 0) =
1
√
2πσ
exp
−
x 2
2σ 2
,
x ∈ [−1, 1],
(9.35)
∂
∂x
u(−1, t) = 0,
t ∈ (0, T ],
(9.36)
∂
∂x
u(1, t) = 0,
t ∈ (0, T ] .
(9.37)