9.3 Exercises
303
K[N,N] = (beta/dx**2)*(-2)
return K
import odespy
solver = odespy.BackwardEuler(rhs, f_is_linear=True, jac=K)
solver = odespy.ThetaRule(rhs, f_is_linear=True, jac=K, theta=0.5)
solver.set_initial_condition(U_0)
T = 1*60*60
N_t = int(round(T/dt))
time_points = linspace(0, T, N_t+1)
u, t = solver.solve(time_points)
The file rod_BE.py has all the details and shows a movie of the solution. We can
run it with any Δt we want, its size just impacts the accuracy of the first steps.
Odespy solvers apply dense matrices!
Looking at the entries of the K matrix, we realize that there are at maximum
three entries different from zero in each row. Therefore, most of the entries
are zeroes. The Odespy solvers expect dense square matrices as input, here
with (N + 1) × (N + 1) elements. When solving the linear systems, a lot
of storage and work are spent on the zero entries in the matrix. It would be
much more efficient to store the matrix as a tridiagonal matrix and apply a
specialized Gaussian elimination solver for tridiagonal systems. Actually, this
reduces the work from the order N 3 to the order N.
In one-dimensional diffusion problems, the savings of using a tridiagonal
matrix are modest in practice, since the matrices are very small anyway. In
two- and three-dimensional PDE problems, however, one cannot afford dense
square matrices. Rather, one must resort to more efficient storage formats and
algorithms tailored to such formats, but this is beyond the scope of the present
text.
9.3 Exercises
Exercise 9.1: Simulate a Diffusion Equation by Hand
Consider the problem given by (9.9), (9.10) and (9.14). Set N = 2 and compute
u
0
i , u 1
i and u 2
i by hand for i = 0, 1, 2. Use these values to construct a test
function for checking that the implementation is correct. Copy useful functions
from test_diffusion_pde_exact_linear.py and make a new test function
test_diffusion_hand_calculation.
Filename: test_rod_hand_calculations.py.
Exercise 9.2: Compute Temperature Variations in the Ground
The surface temperature at the ground shows daily and seasonal oscillations. When
the temperature rises at the surface, heat is propagated into the ground, and the
coefficient β in the diffusion equation determines how fast this propagation is. It
takes some time before the temperature rises down in the ground. At the surface, the
303
K[N,N] = (beta/dx**2)*(-2)
return K
import odespy
solver = odespy.BackwardEuler(rhs, f_is_linear=True, jac=K)
solver = odespy.ThetaRule(rhs, f_is_linear=True, jac=K, theta=0.5)
solver.set_initial_condition(U_0)
T = 1*60*60
N_t = int(round(T/dt))
time_points = linspace(0, T, N_t+1)
u, t = solver.solve(time_points)
The file rod_BE.py has all the details and shows a movie of the solution. We can
run it with any Δt we want, its size just impacts the accuracy of the first steps.
Odespy solvers apply dense matrices!
Looking at the entries of the K matrix, we realize that there are at maximum
three entries different from zero in each row. Therefore, most of the entries
are zeroes. The Odespy solvers expect dense square matrices as input, here
with (N + 1) × (N + 1) elements. When solving the linear systems, a lot
of storage and work are spent on the zero entries in the matrix. It would be
much more efficient to store the matrix as a tridiagonal matrix and apply a
specialized Gaussian elimination solver for tridiagonal systems. Actually, this
reduces the work from the order N 3 to the order N.
In one-dimensional diffusion problems, the savings of using a tridiagonal
matrix are modest in practice, since the matrices are very small anyway. In
two- and three-dimensional PDE problems, however, one cannot afford dense
square matrices. Rather, one must resort to more efficient storage formats and
algorithms tailored to such formats, but this is beyond the scope of the present
text.
9.3 Exercises
Exercise 9.1: Simulate a Diffusion Equation by Hand
Consider the problem given by (9.9), (9.10) and (9.14). Set N = 2 and compute
u
0
i , u 1
i and u 2
i by hand for i = 0, 1, 2. Use these values to construct a test
function for checking that the implementation is correct. Copy useful functions
from test_diffusion_pde_exact_linear.py and make a new test function
test_diffusion_hand_calculation.
Filename: test_rod_hand_calculations.py.
Exercise 9.2: Compute Temperature Variations in the Ground
The surface temperature at the ground shows daily and seasonal oscillations. When
the temperature rises at the surface, heat is propagated into the ground, and the
coefficient β in the diffusion equation determines how fast this propagation is. It
takes some time before the temperature rises down in the ground. At the surface, the