302
9 Solving Partial Differential Equations
A N,N−1 = −Δt
2β
Δx 2
(9.26)
A N,N = 1 + Δt
2β
Δx 2
(9.27)
If we want to apply general methods for systems of ODEs on the form u =
f (u, t), we can assume a linear f (u, t) = Ku. The coefficient matrix K is found
from the right-hand side of (9.16)–(9.18) to be
K 1,1 = 0
(9.28)
K i,i−1 =
β
Δx 2 , i = 2, . . . , N − 1
(9.29)
K i,i+1 =
β
Δx 2 , i = 2, . . . , N − 1
(9.30)
K i,i = −
2β
Δx 2 , i = 2, . . . , N − 1
(9.31)
K N,N−1 =
2β
Δx 2
(9.32)
K N,N = −
2β
Δx 2
(9.33)
We see that A = I − Δt K.
To implement the Backward Euler scheme, we can either fill a matrix and call a
linear solver, or we can apply Odespy. We follow the latter strategy. Implicit methods
in Odespy need the K matrix above, given as an argument jac (Jacobian of f ) in the
call to odespy.BackwardEuler. Here is the Python code for the right-hand side of
the ODE system (rhs) and the K matrix (K) as well as statements for initializing
and running the Odespy solver BackwardEuler (in the file rod_BE.py):
def rhs(u, t):
N = len(u) - 1
rhs = zeros(N+1)
rhs[0] = dsdt(t)
for i in range(1, N):
rhs[i] = (beta/dx**2)*(u[i+1] - 2*u[i] + u[i-1]) + \
g(x[i], t)
rhs[N] = (beta/dx**2)*(2*u[i-1] + 2*dx*dudx(t) -
2*u[i]) + g(x[N], t)
return rhs
def K(u, t):
N = len(u) - 1
K = zeros((N+1,N+1))
K[0,0] = 0
for i in range(1, N):
K[i,i-1] = beta/dx**2
K[i,i] = -2*beta/dx**2
K[i,i+1] = beta/dx**2
K[N,N-1] = (beta/dx**2)*2