Page - 293 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

9.2 FiniteDifferenceMethods 293 9.2.2 ConstructionofaTestProblemwithKnownDiscrete Solution At this point, it is tempting to implement a real physical case and run it. However, PDEsconstituteanon-trivial topicwheremathematicalandprogrammingmistakes come easy. A better start is therefore to address a carefully designed test example wherewecancheckthat themethodworks.Themostattractiveexamplesfor testing implementationsare those without approximationerrors, because we know exactly whatnumberstheprogramshouldproduce.It turnsout thatsolutionsu(x,t) thatare linear in time and in space can be exactly reproduced by most numerical methods forPDEs. Acandidatesolutionmightbe u(x,t)= (3t+2)(x−L). Inserting thisu in thegoverningequationgives 3(x−L)=0+g(x,t) ⇒ g(x,t)=3(x−L). What about the boundary conditions? We realize that ∂u/∂x = 3t + 2 for x = L, which breaks the assumption of ∂u/∂x = 0 at x = L in the formulation of the numerical method above. Moreover, u(0,t) = −L(3t + 2), so we must set s(t)=−L(3t+2)and s′(t)=−3L. Finally, the initial condition dictates I(x)= 2(x−L), but recall that we musthaveu0 = s(0), andui = I(xi), i=1,.. .,N: it is important thatu0 starts out at the right value dictated by s(t) in case I(0) is not equal thisvalue. First we need to generalizeourmethod to handle∂u/∂x=γ =0atx=L. We thenhave uN+1(t)−uN−1(t) 2Δx =γ ⇒ uN+1 =uN−1 +2γΔx, which inserted in (9.7)gives duN(t) dt =β2uN−1(t)+2γΔx−2uN(t) Δx2 +gN(t). (9.14) 9.2.3 Implementation:ForwardEulerMethod In particular, we may use the Forward Euler method as implemented in the general function ode_FE in the module ode_system_FE from Sect.8.3.6. The ode_FE function needs a specification of the right-hand side of the ODE system. This is a matter of translating (9.9), (9.10), and (9.14) to Python code (in file test_diffusion_pde_exact_linear.py): def rhs(u, t): N = len(u) - 1 rhs = np.zeros(N+1) rhs[0] = dsdt(t)