Seite - 176 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

176 5 SolvingPartialDifferentialEquations In the general case (5.16)–(5.18), the coefficientmatrix is an .N C1/ .N C1/ matrixwithzeroentries, except for A1;1 D1 (5.22) Ai;i 1 D t ˇ x2 ; i D2;:: :;N 1 (5.23) Ai;iC1 D t ˇ x2 ; i D2;:: :;N 1 (5.24) Ai;i D1C2 t ˇ x2 ; i D2;:: :;N 1 (5.25) AN;N 1 D t 2ˇ x2 (5.26) AN;N D1C t 2ˇ x2 (5.27) If wewant to apply general methods for systems of ODEs on the formu0 D f.u;t/, we can assumea linearf.u;t/ DKu. The coefficientmatrixK is found fromthe right-handsideof (5.16)–(5.18) tobe K1;1 D0 (5.28) Ki;i 1 D ˇ x2 ; i D2;:: :;N 1 (5.29) Ki;iC1 D ˇ x2 ; i D2;:: :;N 1 (5.30) Ki;i D 2ˇ x2 ; i D2;:: :;N 1 (5.31) KN;N 1 D 2ˇ x2 (5.32) KN;N D 2ˇ x2 (5.33) Wesee thatADI tK. To implement theBackwardEuler scheme,we can either fill amatrix and call a linear solver, or we can apply Odespy. We follow the latter strategy. Implicit methods inOdespyneed theKmatrix above, givenas an argumentjac (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 theK matrix (K) as well as state- ments for initializing and running theOdespy 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