304
9 Solving Partial Differential Equations
temperature has then fallen. We are interested in how the temperature varies down
in the ground because of temperature oscillations on the surface.
Assuming homogeneous horizontal properties of the ground, at least locally,
and no variations of the temperature at the surface at a fixed point of time, we
can neglect the horizontal variations of the temperature. Then a one-dimensional
diffusion equation governs the heat propagation along a vertical axis called x. The
surface corresponds to x = 0 and the x axis point downwards into the ground. There
is no source term in the equation (actually, if rocks in the ground are radioactive, they
emit heat and that can be modeled by a source term, but this effect is neglected here).
At some depth x = L we assume that the heat changes in x vanish, so ∂u/∂x = 0
is an appropriate boundary condition at x = L. We assume a simple sinusoidal
temperature variation at the surface:
u(0, t) = T 0 + T a sin
2π
P
t
,
where P is the period, taken here as 24 h (24 · 60 · 60 s). The β coefficient may be
set to 10 −6 m 2 /s. Time is then measured in seconds. Set appropriate values for T 0
and T a .
a) Show that the present problem has an analytical solution of the form
u(x, t) = A + Be
−rx sin(ωt − rx),
for appropriate values of A, B, r, and ω.
b) Solve this heat propagation problem numerically for some days and animate
the temperature. You may use the Forward Euler method in time. Plot both
the numerical and analytical solution. As initial condition for the numerical
solution, use the exact solution during program development, and when the
curves coincide in the animation for all times, your implementation works, and
you can then switch to a constant initial condition: u(x, 0) = T 0 . For this
latter initial condition, how many periods of oscillations are necessary before
there is a good (visual) match between the numerical and exact solution (despite
differences at t = 0)?
Filename: ground_temp.py.
Exercise 9.3: Compare Implicit Methods
An equally stable, but more accurate method than the Backward Euler scheme, is
the so-called 2-step backward scheme, which for an ODE u = f (u, t) can be
expressed by
3u n+1 − 4u n + u n−1
2Δt
= f (u
n+1 , t n+1 ) .
The Odespy package offers this method as odespy.Backward2Step. The purpose
of this exercise is to compare three methods and animate the three solutions:
1. The Backward Euler method with Δt = 0.001
2. The backward 2-step method with Δt = 0.001
9 Solving Partial Differential Equations
temperature has then fallen. We are interested in how the temperature varies down
in the ground because of temperature oscillations on the surface.
Assuming homogeneous horizontal properties of the ground, at least locally,
and no variations of the temperature at the surface at a fixed point of time, we
can neglect the horizontal variations of the temperature. Then a one-dimensional
diffusion equation governs the heat propagation along a vertical axis called x. The
surface corresponds to x = 0 and the x axis point downwards into the ground. There
is no source term in the equation (actually, if rocks in the ground are radioactive, they
emit heat and that can be modeled by a source term, but this effect is neglected here).
At some depth x = L we assume that the heat changes in x vanish, so ∂u/∂x = 0
is an appropriate boundary condition at x = L. We assume a simple sinusoidal
temperature variation at the surface:
u(0, t) = T 0 + T a sin
2π
P
t
,
where P is the period, taken here as 24 h (24 · 60 · 60 s). The β coefficient may be
set to 10 −6 m 2 /s. Time is then measured in seconds. Set appropriate values for T 0
and T a .
a) Show that the present problem has an analytical solution of the form
u(x, t) = A + Be
−rx sin(ωt − rx),
for appropriate values of A, B, r, and ω.
b) Solve this heat propagation problem numerically for some days and animate
the temperature. You may use the Forward Euler method in time. Plot both
the numerical and analytical solution. As initial condition for the numerical
solution, use the exact solution during program development, and when the
curves coincide in the animation for all times, your implementation works, and
you can then switch to a constant initial condition: u(x, 0) = T 0 . For this
latter initial condition, how many periods of oscillations are necessary before
there is a good (visual) match between the numerical and exact solution (despite
differences at t = 0)?
Filename: ground_temp.py.
Exercise 9.3: Compare Implicit Methods
An equally stable, but more accurate method than the Backward Euler scheme, is
the so-called 2-step backward scheme, which for an ODE u = f (u, t) can be
expressed by
3u n+1 − 4u n + u n−1
2Δt
= f (u
n+1 , t n+1 ) .
The Odespy package offers this method as odespy.Backward2Step. The purpose
of this exercise is to compare three methods and animate the three solutions:
1. The Backward Euler method with Δt = 0.001
2. The backward 2-step method with Δt = 0.001