Seite - 302 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Band Second Edition

302 9 SolvingPartialDifferentialEquations AN,N−1 =−Δt 2β Δx2 (9.26) AN,N =1+Δt 2β Δx2 (9.27) If we want to apply general methods for systems of ODEs on the form u′ = f(u,t), we can assume a linearf(u,t) =Ku. The coefficient matrixK is found fromthe right-handsideof (9.16)–(9.18)to be K1,1 =0 (9.28) Ki,i−1 = β Δx2 , i=2,.. .,N−1 (9.29) Ki,i+1 = β Δx2 , i=2,.. .,N−1 (9.30) Ki,i =− 2β Δx2 , i=2,.. .,N−1 (9.31) KN,N−1 = 2β Δx2 (9.32) KN,N =− 2β Δx2 (9.33) We see thatA= I−ΔtK. To implement the Backward Euler scheme, we can either fill a matrix and call a linearsolver,orwecanapplyOdespy.Wefollowthelatterstrategy.Implicitmethods inOdespyneedtheKmatrixabove,givenasanargumentjac(Jacobianoff) in the call toodespy.BackwardEuler.Here is thePythoncodefor theright-handsideof the ODE system (rhs) and the K matrix (K) as well as statements for initializing andrunningtheOdespysolverBackwardEuler(in thefilerod_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