Page - 167 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

5.1 FiniteDifferenceMethods 167 This is a matter of translating (5.9), (5.10), and (5.14) to Python code (in file test_diffusion_pde_exact_linear.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[N-1] + 2*dx*dudx(t) - 2*u[N]) + g(x[N], t) return rhs def u_exact(x, t): return (3*t + 2)*(x - L) def dudx(t): return (3*t + 2) def s(t): return u_exact(0, t) def dsdt(t): return 3*(-L) def g(x, t): return 3*(x-L) Note thatdudx(t) is the functionrepresentingthe parameter in (5.14).Alsonote that the rhs function relies on access to global variables beta, dx, L, and x, and global functionsdsdt,g, anddudx. Weexpect the solution tobe correct regardless ofN and t, sowecan choose a smallN ,N D4, and t D0:1.A test functionwithN D4goes like def test_diffusion_exact_linear(): global beta, dx, L, x # needed in rhs L = 1.5 beta = 0.5 N = 4 x = linspace(0, L, N+1) dx = x[1] - x[0] u = zeros(N+1) U_0 = zeros(N+1) U_0[0] = s(0) U_0[1:] = u_exact(x[1:], 0) dt = 0.1 print dt u, t = ode_FE(rhs, U_0, dt, T=1.2)