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

9.2 FiniteDifferenceMethods 291 ThePDEisvalidatall spatialpointsx∈Ω,butwemayrelax thisconditionand demandthat it is fulfilledat the internalmesh pointsonly,x1,.. .,xN−1: ∂u(xi,t) ∂t =β∂ 2u(xi,t) ∂x2 +g(xi,t), i=1,.. .,N−1 . (9.5) Now, at any point xi we can approximate the second-order derivative by a finite difference: ∂2u(xi,t) ∂x2 ≈ u(xi+1,t)−2u(xi,t)+u(xi−1,t) Δx2 . (9.6) It is commonto introducea shortnotationui(t) foru(xi,t), i.e.,uapproximatedat some mesh point xi in space. With this new notation we can, after inserting (9.6) in (9.5),writeanapproximationto the PDEat meshpoint (xi,t) as dui(t) dt =βui+1(t)−2ui(t)+ui−1(t) Δx2 +gi(t), i=1,.. .,N−1 . (9.7) Note thatwe haveadopted thenotationgi(t) forg(xi,t) too. What is (9.7)? This is nothing but a system of ordinary differential equations in N−1unknownsu1(t),.. .,uN−1(t)! Inotherwords,withaidofthefinitedifference approximation(9.6), we have reduced the single PDE to a system of ODEs, which weknowhowtosolve. In the literature, this strategy is called the methodof lines. We need to look into the initial and boundary conditions as well. The initial condition u(x,0) = I(x) translates to an initial condition for every unknown function ui(t): ui(0) = I(xi), i = 0,.. .,N. At the boundary x = 0 we need an ODE in our ODE system, which must come from the boundary condition at this point. The boundary condition reads u(0,t) = s(t). We can derive an ODE from this equation by differentiating both sides: u′0(t) = s′(t). The ODE system abovecannotbeusedforu′0 since thatequation involvessomequantityu′−1 outside the domain. Instead, we use the equationu′0(t)= s′(t) derived from the boundary condition.Forthisparticularequationwealsoneedtomakesuretheinitialcondition isu0(0)= s(0) (otherwisenothingwill happen:we getu=283K forever). We remark that a separate ODE for the (known) boundary conditionu0 = s(t) is not strictly needed. We can just work with the ODE system foru1,.. .,uN, and in the ODE foru0, replaceu0(t) by s(t). However, these authors prefer to have an ODEforeverypointvalueui, i=0,.. .,N,whichrequiresformulatingtheknown boundaryatx= 0 asan ODE. The reason for including the boundaryvalues in the ODE system is that the solution of the system is then the complete solution at all mesh points,which is convenient, since special treatment of the boundaryvalues is thenavoided. The condition ∂u/∂x = 0 at x = L is a bit more complicated, but we can approximatethe spatialderivativebya centeredfinitedifference: ∂u ∂x ∣ ∣ ∣ ∣ i=N ≈ uN+1 −uN−1 2Δx =0 .