Rk4

Ain Shams University College of Engineering Department of Electrical Engineering 2nd Year Electrical System Dynamics and Control Components

RUNGE-KUTTA 4th ORDER METHOD Under Supervision of Professor Dr. A. A. Wahdan Mohammed Shadied Abbas 2 Year Electrical Engineering Group B Section : 9 Bench No. : 6 nd

Table of Contents : I – Differential Equations .................................................... 3 i- Definition ....................................................................... 3 ii- Ordinary Differential Equations ................................... 4 iii- Numerical Ordinary Differential Equations ................ 5 II – Runge-Kutta Methods ................................................... 6 i- Introduction .................................................................... 6 ii- Fourth Order Runge-Kutta Method .............................. 6 iii- Computer-Aided 4th Order Runge-Kutta Method ........ 9 1- The Programming Language Python ...................... 9 2- The Program Code .................................................. 9 3- A Sample Runtime ................................................ 14 III – Utilising Runge-Kutta Method In Solving Higher Order Differential Equations ............................................... 15 IV – References And Sources …........................................ 18

2

I – Differential Equations : i- Definition : A differential equation is a mathematical equation for an unknown function of one or several variables that relates the values of the function itself and its derivatives of various orders. Differential equations play a prominent role in engineering, physics, economics and other disciplines. A simplified real world example of a differential equation is modeling the acceleration of a ball falling through the air (considering only gravity and air resistance). The ball's acceleration towards the ground is the acceleration due to gravity minus the deceleration due to air resistance. Gravity is a constant but air resistance is proportional to the ball's velocity. This means the ball's acceleration is dependent on its velocity. Because acceleration is the derivative of velocity, solving this problem requires a differential equation. Differential equations are mathematically studied from several different perspectives, mostly concerned with their solutions, functions that make the equation hold true. Only the simplest differential equations admit solutions given by explicit formulas. Many properties of solutions of a given differential equation may be determined without finding their exact form. If a self-contained formula for the solution is not available, the solution may be numerically approximated using computers. The theory of dynamical systems puts emphasis on qualitative analysis of systems described by differential equations, while many numerical methods have been developed to determine solutions with a given degree of accuracy.

3

ii- Ordinary Differential Equations : In mathematics, an ordinary differential equation (or ODE) is a relation that contains functions of only one independent variable, and one or more of its derivatives with respect to that variable. A simple example is Newton's second law of motion, which leads to the differential equation: ..

m x = F t

For the motion of a particle of constant mass m. In general, the force F depends upon the position of the particle x(t) at time t, and thus the unknown function x(t) appears on both sides of the differential equation, as is indicated in the notation F(x(t)). Ordinary differential equations arise in many different contexts including geometry, mechanics, astronomy and population modeling. Many famous mathematicians have studied differential equations and contributed to the field, including Newton, Leibniz, the Bernoulli family, Riccati, Clairaut, d'Alembert and Euler. Much study has been devoted to the solution of ordinary differential equations. In the case where the equation is linear, it can be solved by analytical methods. Unfortunately, most of the interesting differential equations are non-linear and, with a few exceptions, cannot be solved exactly. Approximate solutions are arrived at using computer approximations (see numerical ordinary differential equations).

4

iii- Numerical Ordinary Differential Equations : Numerical ordinary differential equations is the part of numerical analysis which studies the numerical solution of ordinary differential equations (ODEs). This field is also known under the name numerical integration, but some people reserve this term for the computation of integrals. Many differential equations cannot be solved analytically, in which case we have to satisfy ourselves with an approximation to the solution. The algorithms studied here can be used to compute such an approximation. An alternative method is to use techniques from calculus to obtain a series expansion of the solution. Ordinary differential equations occur in many scientific disciplines, for instance in mechanics, chemistry, biology, and economics. In addition, some methods in numerical partial differential equations convert the partial differential equation into an ordinary differential equation, which must then be solved. Some elementary methods were made to address this problem numerically such as The Euler Method and The Backward Euler Method. After that, other methods have emerged (which are generally more accurate and efficient).

Fig.(1) : Illustration of numerical integration for the differential equation y' = y,y(0) = 1. Blue: the Euler method, Green: the midpoint method and Red: the exact solution, y = et. The step size is h = 1.0. Figure (1)

5

II – Runge-Kutta Methods : i- Introduction : A method of numerically integrating ordinary differential equations by using a trial step at the midpoint of an interval to cancel out lower-order error terms. The second-order formula is k 1 =h f  x n , y n  k 2=h f  x n1/ 2 h , y n 1/ 2 k 1  3 y n1 = y n k 2 O h 

ii- Fourth Order Runge-Kutta Method : The fourth order Runge-Kutta method has been derived from Taylor formula where 2 y  x n1 = y n h f n h / 2 ! [ f x n  f n  f y n ] so it's possible to develop formula equivalent to third, fourth, fifth or even higher order Taylor formula, which do not involve derivatives of f . The derivation of a similar formula in the case of five term Taylor series formula h 2 ' ' h3 ' ' ' h4 4 y n1 = y n h y n  y n  y n  y n 2! 3! 4! '

where derivatives y'' , y''' and y(4) are obtained by successively differentiating y'=f(x,y) is too lengthy and tedious but not difficult.

6

The formula obtained after the derivation is 1 y n1 = y n  w 1 2w 2 2w 3 w 4  6 where ; w 1 =hf  x n , y n  1 1 w 2 =hf  x n  h , y n  w 1  2 2 1 1 w3=h f  x n  h , y n  w 2  2 2 w 4=hf  x n h , y n w 3 

This is the classical Runge-Kutta formula. It's equivalent to five term Taylor series above formula and hence it's fourth order method.

Example : Apply the Runge-Kutta method to the IVP ' y =1− x4y ; y 0=1 to find y(0.2). Take h = 0.1 and use six digits.

Solution : n 0

x 0.00 0.05 0.05 0.10

y 1.00000 1.25000 1.29750 1.30700

w=h(1-x+4y) 0.500000 0.595000 0.614000 0.612800

1 y1 =1 2.5228=1.420466 6

Δ y 0.50000 1.19000 0.22800 0.612800

2.5228

7

n 1

x 0.10 0.15 0.15 0.20

y 1.42046 1.74955 1.81282 1.82553

w=h(1-x+4y) 0.65818 0.78482 0.81014 0.81021

1 y 2 =1.420466 4.45832=2.44469 6

Δ y 0.65818 1.56964 1.42029 0.81021

4.45832

So, y(0.2) = 2.44469

8

iii- Computer-Aided 4th Order Runge-Kutta Method : 1- The Programming Language Python : Python is a general-purpose high-level programming language that firstly appeared in 1991. Its design philosophy emphasizes code readability. Python's core syntax and semantics are minimalistic, while the standard library is large and comprehensive. Its use of whitespace as block delimiters is unusual among popular programming languages. Python supports multiple programming paradigms (primarily object oriented, imperative, and functional) and features a fully dynamic type system and automatic memory management, similar to Perl, Ruby, Scheme, and Tcl. Like other dynamic languages, Python is often used as a scripting language. The language has an open, community-based development model managed by the non-profit Python Software Foundation, which maintains the de facto standard definition of the language in CPython, the reference implementation. Python runs on Windows, Linux/Unix, Mac OS X, OS/2, Amiga, Palm Handhelds, and Nokia mobile phones. Python has also been ported to the Java and .NET virtual machines. For more information about Python, you may refer to "http://python.org/"

2- The Program Code : 9

H:\Users\Mohammed\Desktop\Runge\RK4.py 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

# In Th e N ame OF Al lah , M ost Gr aci ous , M ost Me rci ful # Th e f oll owi ng cod e r epr ese nts a 1st or der OD E s olv er usi ng Run ge Kut ta 4 me tho d # Th e e qua tio n s hou ld be wri tte n a s y ' = f( x,y ) a nd it sho uld be va lid ac cor din g t o p yth on' s s ynt ax # Wr itt en by Moh amm ed Abb as "ab bas @ie ee. org " u sin g p yth on scr ipt ing la ngu age v2 .6. 1 # La st Edi t : 9: 52 PM , 1 6th of Ap ril 20 09

d ef RK4 ():

# u ser 's inp ut wit h p rob lem gi ven s : pri nt "Nu mer ica l 1 st ord er ODE so lve r u sin g R ung e-K utt a 4 me tho d" pri nt 'Al gor ith m w rit ten by Mo ham med Ab bas "a bba s@i eee .or g"\ nUs ing py tho n sc rip tin g l ang uag e v 2. 6.1 ' pri nt "\n Imp ort ant no tes :\ n1- Wh en ent eri ng the eq uat ion , p lea se mak e su re you r u se * f or mul tip lic ati on and ** fo r e xpo nen tia tio n ex : x ^2 sho uld be wr itt en as x** 2. \n2 - U se par ent hes es () to ski p or der of pr ece den ce. "

17

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

# e qua tio n f (x, y) is ent ere d a s a st rin g : eqn = raw _in put ("\ nPl eas e e nte r t he equ ati on as y' = " ) # i nit ial co ndi tio ns, st ep siz e a nd the po int of so lut ion : x0 = i npu t(" Ple ase en ter x0 = ") y0 = i npu t(" Ple ase en ter y( x0) = ") h = in put ("P lea se def ine th e s tep si ze h = ")

s = in put ("P lea se spe cif y

the po int at wh ich yo u n eed th e s olu tio n y (x) >

> ") 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

# n re pre sen ts how ma ny tim es the it era tio n i s d one us ing th e f oll owi ng f or loo p t ech niq ue : n = (s *1. 0)/ h n = in t(n ) # f our fl oat ing po int ar ray s f or eac h c olu mn of the it era tio n i nit ial ly z ero ed : xA = [ 0.0 ,0. 0,0 .0, 0.0 ] yA = [ 0.0 ,0. 0,0 .0, 0.0 ] wA = [ 0.0 ,0. 0,0 .0, 0.0 ] dyA = [0. 0,0 .0, 0.0 ,0. 0]

H:\Users\Mohammed\Desktop\Runge\RK4.py 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64

65 66 67 68 69 70

71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98

# t his is th e b egi nni ng of the it era tio n p roc ess , ' for ' l oop sh oul d g o a s m any ti mes as n is : for i in ran ge( n): # her e w e s imp ly tak e t he val ues of in iti al con dit ion va ria ble s x 0,y 0 ent ere d b y u ser : xA [0] =x0 yA [0] =y0 #c hec k #p rin t x A[0 ],y A[0 ] # the n w e a ssi gn the m t o v ari abl es cal led x and y to hel p e val () fun cti on eva lua te the en ter ed equ ati on as pyt hon 's val id exp res sio n, it the n h app ens ea ch and ev ery ti me we' d l ike to us e e val () fun cti on wit h o the r v alu es :

x= xA[ 0] y= yA[ 0] # her e w e u se eva l() fu nti on to eva lua te the en ter ed exp res sio n w ith x0, y0 val ues to ge t w 1(0 ) a nd the n w e u se w1( 0) to get th e n ext val ue of y : wA [0] = h * e val (eq n) yA [1] = y A[0 ] + (0 .5* wA[ 0]) #c hec k #p rin t w A[0 ],y A[1 ] # her e's ho w t he nex t v alu e o f x is ob tai ned an d h ow we use th e pre vio us tec hni que ag ain : xA [1] = x 0+( 0.5 *h) x= xA[ 1] y= yA[ 1] #c hec k #p rin t x A[1 ] wA [1] = h * e val (eq n) yA [2] = y A[0 ] + (0 .5* wA[ 1]) #c hec k #p rin t w A[1 ],y A[2 ]

H:\Users\Mohammed\Desktop\Runge\RK4.py 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151

152

xA [2] = x 0+( 0.5 *h) x= xA[ 2] y= yA[ 2] #c hec k #p rin t x A[2 ] wA [2] = h * e val (eq n) yA [3] = y A[0 ] + (0 .5* wA[ 2]) #c hec k #p rin t w A[2 ],y A[3 ] xA [3] = x 0+h x= xA[ 3] y= yA[ 3] #c hec k #p rin t x A[3 ]

wA [3] = h * e val (eq n) #c hec k #p rin t w A[3 ] # fin all y w e c alc ula te the va lue of dy fo r t he ite rat ion : dy A[0 ] dy A[1 ] dy A[2 ] dy A[3 ]

= = = =

wA [0] 2* wA[ 1] 2* wA[ 2] wA [3]

#c hec k #p rin t d yA[ 0], dyA [1] ,dy A[2 ],d yA[ 3] su m = dy A[0 ]+d yA[ 1]+ dyA [2] +dy A[3 ] # thi s i s h ow we get yn +1 : re sul t = yA [0] + ( sum /6. 0) # bef ore th e i ter ati on end s, we' d l ike to us e s ome va lue fr om it for the ne xt one : x0 =xA [3] y0 =re sul t

# x n+1 # y n+1

#c hec k #p rin t x 0,y 0 # f unc tio n r etu rns th e f ina l r esu lt of the lo op whi ch' s t he sol uti on of t he DE at the po int of in ter est an d t he val ue of thi s p oin t : # u nli ke C a nd som e o the r l ang uag es, py tho n i s c apa ble of re tur nin g m ult ipl e s tuf f n ot jus t o ne var iab le wit hou t t he exp lic it use of p oin ter s.

H:\Users\Mohammed\Desktop\Runge\RK4.py 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179

ret urn re sul t,s # e nd of RK4 () fun cti on

# he re we sim ply ca ll the fu nct ion to re cei ve the de sir ed val ues : f ina l,s = RK4 () # th en we rou nd the so lut ion to 5 dig its : f ina l = ro und (fi nal ,5) # th e o utp ut : p rin t " \ny (%s ) = %f " %(s ,fi nal ) # sc rip t t erm ina tio n : e sc = r aw_ inp ut( "\n \nP res s a ny key to en d . .." )

3- A Sample Runtime : A sample successful runtime of the code clearly represents the previous example. Numerical 1st order ODE solver using Runge-Kutta 4 method Algorithm written by Mohammed Abbas "[email protected]" Using python scripting language v 2.6.1 Important notes : 1- When entering the equation, please make sure your use * for multiplication and ** for exponentiation ex: x^2 should be written as x**2. 2- Use parentheses () to skip order of precedence. Please Please Please Please Please

enter the equation as y' = 1-x+(4*y) enter x0 = 0 enter y(x0) = 1 define the step size h = 0.1 specify the point at which you need the solution y(x) >> 0.2

y(0.2) = 2.444690 Press any key to end ...

14

III – Utilizing Runge-Kutta Method In Solving Higher Order Differential Equations : Any differential equation of order n can be written as a system of n first-order differential equations. Given an explicit ordinary differential equation of order n and dimension 1, '

''

F  x , y , y , y ,..... , y

n−1

= y

n

we define a new family of unknown functions y n := y

n−1

We can then rewrite the original differential equation as a system of differential equations with order 1 and dimension n. '

y1 = y 2 '

y 2= y 3 . . . '

y n−1 = y n '

y n =F  y n ,.... , y 1 , x

which can be written concisely in vector notation as '

y = F  x , y ;

with

y := y ,..... , y n 

15

Suppose we obtained a system from a higher order DE to be as follows : ' ' x = f t , x , y , y = g t , x , y ; x 0= x 0 , y 0= y0 then the Runge-Kutta formulas for the system will be : 1 x n1= x n  v 12 v 22 v 3 v 4  6 1 y n1 = y n  w 1 2w 2 2w 3w 4  6

where : v1 =hf t n , x n , y n  w 1 =hg t n , x n , y n  1 1 1 v 2 =hf t n h , x n v1 , y n  w 1  2 2 2 1 1 1 w 2 =hg t n  h , x n  v 1 , y n w 1  2 2 2 1 1 1 v 3=hf t n  h , x n  v 2 , y n w 2  2 2 2 1 1 1 w 3 =hg t n  h , x n  v 2 , y n  w 2  2 2 2 v 4 =hf t n h , x n v 3 , y n w 3  w 3 =hg t n h , x n v 3 , y n w 3 

As the original formula error is O h5  , so the accumulated formula error is O h 4  . This method is sufficiently accurate.

16

Example : Use Runge-Kutta method with step size h = 0.25 to compute x(0.5) for the IVP x ' ' −tx ' − x=0 ; x 0=0 , x ' 0=1 Use five digits.

Solution : The equivalent system is : '

x = y= f ' y = xty=g x 0=0 , y 0=1

The iteration goes as follows : n 0

n 1

t 0.000 0.125 0.125 0.250

x y V=hf 0.00000 1.0000 0.25000 0.12500 1.0000 0.25000 0.12500 1.0313 0.25783 0.25783 1.0635 0.26588 x 1=00.25526=0.25526 y1 =10.063815=1.0638

W=hg 0.00000 0.06250 0.06347 0.13093

t x y V=hf W=hg 0.250 0.25526 1.0638 0.26595 0.13030 0.375 0.38824 1.1290 0.28225 0.20290 0.375 0.39639 1.1653 0.29133 0.20848 0.500 0.54659 1.2723 0.31808 0.29569 x 0.5≈ x 2 =0.255260.28853=0.54379 ' x 0.5= y 0.5≈ y 2 =1.06380.20813=1.2719

Δx 0.25000 0.50000 0.51566 0.26588

Δy 0.00000 0.12500 0.12696 0.13093

0.28853 0.063815

Δx 0.26595 0.56450 0.58266 0.31808

Δy 0.13030 0.40580 0.41696 0.29569

0.28853 0.20813

17

Although we took n = 2 only, the total error made in computing x(0.5) is 0.00004 which is much better than any of the second order methods with n = 5. The code can be moderately modified to perform for this method.

IV – References And Sources : 1- Modeling and Analysis of Dynamic Systems Publisher: Wiley; 3rd edition (August 6, 2001) ISBN-13: 978-0471394426 By Charles M. Close, Dean K. Frederick and Jonathan C. Newell 2- Engineering Mathematics (3) Part II, ASU 2nd Year Electrical Publisher: Dar Teba Press By Professor Dr. Reda El Barkouky 3- Learning Python 3rd Edition Publisher: O'Reilly Media, Inc.; 3rd edition (July 16, 2008) ISBN-13: 978-0596513986 By Mark Lutz 4- The Wikipedia, http://wikipedia.org/ 5- Wolfram Mathworld, http://mathworld.wolfram.com/

18