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) 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)

2 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 applications to describe

3 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)..(x xn)fn+1( )(n+1)!

6 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).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).

7 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 . 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.

8 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 . 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.

9 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.(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( ).

10 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 . z e r o s (N)1516N u m e r i c a lS o l u t i o n [0] =117x [0] =018w[0] =11920f o r i i n r a n g e (1,N) :21w[ i ] =w[ i 1]+dx ( t a u ) w[ i 1]22x [ i ] = x [ i 1]+dx23A n a l y t i cS o l u t i o n [ i ] = math . exp ( t a u x [ i ] )242526f i g = p l t . f i g u r e ( f i g s i z e = (8,4) )27# l e f t hand p l o t28ax = f i g.