Numerical Methods for Differential Equations with Python

1 JohnS ButlerNumerical Methods forDifferential Equations withPythonC O N T E N T Si initial value problems61 Numerical solutions to initial value approximation of of Forward Euler for one about ordinary Differential s Sheet222 higher order order Taylor Methods233 runge kutta of Second order runge kutta second order : Midpoint 2nd order runge Kuttaa0= : Heun s order runge kutta kutta fourth s 4th order runge choice of method and Sheet2 374 multi-step of a explicit multistep Derivation of a explicit method three step four step of the implicit multi-step of Adam s step-size multi-step Sheet3 535 consistency,convergence and Step Sheet4 Value Problem Review Questions642 Contents3ii Numerical solutions to boundary value prob-lems696 boundary value of order Value theorems about boundary value Shooting Shooting method for non-linear Difference method80iii Numerical solutions to partial Differential equa-tions847 partial Differential Operators898 parabolic Heat explicit method for the heat implicit (BTCS) method for the Heat implicit (BTCS)

2 For the Heat Nicholson Implicit Crank-Nicholson solution of the Theta General Matrix Boundary Derivative Boundary Truncation Error and and and by the Fourier Series method (von Neumann smethod) for the explicit FTCS for the implicit BTCS for the Crank Nicholson Equations Nicholson Methods1269 elliptic pde s five point approximation of the representation of the five point Matrix form of the discrete :Homogeneous equation with : non-homogeneous equation withzero : Inhomogeneous equation with non-zero and Equations Questions15010 hyperbolic Equations Wave Difference Method for Hyperbolic of the scalar of the Finite Difference Freidrich Lewy Neumann stability for the Forward Neumann stability for the Lax-Friedrich15711 variational Methods -Galerkin bounds of Finite Element methods163A C R O N Y M SIVPI nitial Value ProblemsBVPB oundary Value ProblemsODEO rdinary Differential EquationsPDEP artial Differential EquationsRKRunge Kutta5 Part II N I T I A L VA L U E P R O B L E M S1N U M E R I C A L S O L U T I O N S T O I N I T I A L VA L U EP R O B L E M SDifferential Equations have numerous

3 Applications to describe dy-namics from physics to biology to value problems are subset of ordinary Differential Equation(ODE s) with the formy =f(x)(1)fis a function. The general solution to (1) isy= f(x)dx+c,containing an arbitrary constantc. In order to determine the solutionuniquely it is necessary to impose an initial condition,y(x0) =y0.(2)Example1 Simple ExampleThe Differential equation describes the rate of changeof an oscillating input. The general solution of theequationy =sin(x)(3)is,y= cos(x) +c,with the initial condition,y(0) =2,then it is easy to findc= the desired solution is,y=2 cos(x).The more general ordinary Differential Equation is of the formy =f(x,y),(4)7numerical solutions to initial value problems8is approached in a similar us considery =a(x)y(x) +b(x),The given functionsa(x)andb(x)are assumed continuous for thisequationf(x,z) =a(x)z(x) +b(x),and the general solution can be found using the method of ExampleDifferential Equations of the formy (x) = y(x) +b(x),x x0,(5)where is a given constant andb(x)is a continuousintegrable function has a unique analytic.

4 Multiply-ing the equation (5) by the integrating factore x, wecan reformulated(e xy(x))dx=e xb(x).Integrating both sides fromx0toxwe obtain(e xy(x)) =c+ xx0e tb(t)dt,so the general solution isy(x) =ce x+ xx0e (x t)b(t)dt,withcan arbitrary constantc=e x0y(x0).For a great number of Initial Value Problems there is no knownexact (analytic) solution as the Equations are non-linear, for exam-pley =exy4, or discontinuous or stochastic. There for a numericalmethod is used to approximate the approximation of Numerical approximation of of Forward Euler for one stepThe left hand side of a initial value problemd fdxcan be approximatedbyTaylors theoremexpand about a pointx0giving:f(x1) =f(x0) + (x1 x0)f (x0) + ,(6)where is the truncation error, =(x1 x0)22!

5 F ( ), [x0,x1].(7)Rearranging and lettingh=x1 x0the equation becomesf (x0) =f(x1) f(x0)h h2f ( ).The forward Euler method can also be derived using a variation onthe Lagrange interpolation formula called the divided functionf(x)can be approximated by a polynomial of degreePn(x)and an error term,f(x) =Pn(x) +error,=f(x0) +f[x0,x1](x x0) +f[x0,x1,x2](x x0)(x x1),+..+f[x0, ..,xn] n 1i=0(x xi) +error,wheref[x0,x1] =f(x1) f(x0)x1 x0,f[x0,x1,x2] =f[x1,x2] f[x0,x1]x2 x0,f[x0,x1, ..,xn] =f[x1,x2, ..,xn] f[x0,x1, ..,xn 1]xn x0,DifferentiatingPn(x)P n(x) =f[x0,x1] +f[x0,x1,x2]{(x x0) + (x x1)},+..+f[x0, ..,xn]n 1 i=0(x x0)..(x xn 1)(x xi),and the error becomeserror= (x x0).

6 (x xn)fn+1( )(n+1)!. approximation of Differentiation10 Applying this to define our first derivative, we havef (x) =f[x0,x1] =f(x1) f(x0)x1 x0,this leads us other formulas for computing the derivativesf (x) =f(x1) f(x0)x1 x0+O(h),Euler,f (x) =f(x1) f(x 1)x1 x 1+O(h2), the same method we can get out computational estimates forthe2nd derivativef (x0) =f2 2f1+f0h2+O(h2),f (x0) =f1 2f0+f 1h2+O(h2), numerically solve the first order ordinary Differ-ential Equation (4)y =f(x,y),a x b,the derivativey is approximated bywi+1 wixi+1 xi=wi+1 wih,wherewiis the Numerical approximation Differential Equation is converted to a discretedifference equation with steps of sizeh,wi+1 wih=f(xi,w).

7 Rearranging the difference equation gives the equa-tionwi+1=wi+h f(xi,w),which can be used to approximate the solution atwi+1given information aboutyat example ODE y =sin(x) approximation of Differentiation11 Example4 Applying the Euler formula to the first order equationwith an oscillating input (3)y =sin(x),0 x equation can be approximated using the forwardEuler aswi+1 wih=sin(xi).Rearranging the equation gives the discrete differenceequation with the unknowns on the left and the knowvalues of the rightwi+1=wi+hsin(xi).The Python code bellow implements this differenceequation. The output of the code is shown in # Numerical s o l u t i o n o f a C o s i n e d i f f e r e n t i a le q u a t i o n2i m p o r t numpy a s np3i m p o r t math4i m p o r t m a t p l o t l i b.

8 P y p l o t a s p l t56h = =08b=10910N= i n t ( b a/h )11w=np . z e r o s (N)12x=np . z e r o s (N)13A n a l y t i cS o l u t i o n =np . z e r o s (N)1415# I n i t i a l C o n d i t i o n s16w[0] = [0] =018A n a l y t i cS o l u t i o n [0] = o r i i n r a n g e (1,N) :20w[ i ] =w[ i 1]+h math . s i n ( x [ i 1] )21x [ i ] = x [ i 1]+h22A n a l y t i cS o l u t i o n [ i ]= math . c o s ( x [ i ] )2324f i g = p l t . f i g u r e ( f i g s i z e = (8,4) )2526# l e f t hand p l o t27ax = f i g . a d ds u b p l o t (1,3,1)28p l t . p l o t ( x , w, c o l o r = r e d )29# ax . l e g e n d ( l o c = b e s t )30p l t . t i t l e ( Numerical S o l u t i o n ) approximation of Differentiation123132# r i g h t hand p l o t33ax = f i g.

9 A d ds u b p l o t (1,3,2)34p l t . p l o t ( x , A n a l y t i cS o l u t i o n , c o l o r = b l u e )35p l t . t i t l e ( A n a l y t i c S o l u t i o n )3637# ax . l e g e n d ( l o c = b e s t )38ax = f i g . a d ds u b p l o t (1,3,3)39p l t . p l o t ( x , A n a l y t i cS o l u t i o n w, c o l o r = b l u e )40p l t . t i t l e ( E r r o r )4142# t i t l e , e x p l a n a t o r y t e x t and s a v e43f i g . s u p t i t l e ( S i n e S o l u t i o n , f o n t s i z e =20)44p l t . t i g h tl a y o u t ( )45p l t . s u b p l o t sa d j u s t ( t o p = ) : Python Numerical and Analytical Solution of : Python output: Numerical (left), Analytic (middle) anderror(right) fory =sin(x)Equation3with h= example problem population growth y = population growth can be describe as a firstorder Differential equation of the form:y = y.

10 (8)This has an exact solution ofy=Ce the initial condition of conditiony(0) = approximation of Differentiation13and a rate of change of = analytic solution isy= the Euler formula to the first order equation(8)y = approximated bywi+1 wih= the equation gives the difference equa-tionwi+1=wi+h( ).The Python code below and the output is plotted # Numerical s o l u t i o n o f a d i f f e r e n t i a l e q u a t i o n2i m p o r t numpy a s np3i m p o r t math4i m p o r t m a t p l o t l i b . p y p l o t a s p l t56h = a u = =09b=101011N= i n t ( ( b a ) /h )12w=np . z e r o s (N)13x=np . z e r o s (N)14A n a l y t i cS o l u t i o n =np.