RUNGE-KUTTA METHODS FOR SOLVING ORDINARY DIFFERENTIAL EQUATION

A PROJECT SUBMITTED TO THE
DEPARTMENT OF MATHEMATICS, FACULTY OF PHYSICAL SCIENCE,
KEBBI STATE UNIVERSITY OF SCIENCE AND TECHNOLOGY ALIERO.


IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE AWARD OF DEGREE OF BACHELOR OF SCIENCE (B. Sc. Hons.)
(B.SC. MATHEMATICS)
SEPTEMBER, 2019
DECLARATION
I, Anas Yakubu Bakura, hereby declare that this project Runge-Kutta Methods for Solving Ordinary Differential Equation has been carried out by me under the supervision of Mal. A malik, It has not been presented for award of any degree in any institution. All sources of information are specifically acknowledges by means of reference.
                                                               
..............................................                                                              ........................................
           Signature                                                                                                  Date









CERTIFICATION
This project entitled “Runge-Kutta Methods for Solving Ordinary Differential Equation” by Ibrahim Hajara Makarfi with registration number U11MT1030 meets the requirements governing the award of the degree of Bachelor of Science in Mathematics and is approved for its contribution to knowledge and literary representation.

-------------------------------- ---------------------------     
Mal. A malik               Date
Supervisor


-------------------------------- ---------------------------     
Prof. B. Sani                Date
Head of Department                                                                               

External Examiner
Name: -------------------------                                   ---------------------------     
Signature: ---------------------------                                                          Date





DEDICATION
I humbly dedicated this research work to my beloved parents (Alhaji Yakubu Halliru Bakura and Hajiya Balkisu Halliru) for the love and understanding right from in fancy, till this present stage of life. May Allah bless and grant them a long and healthy life. Ameen

ACKNOWLEDGEMENTS
My sincere appreciation goes to my husband (Alhaji Kabir Rabiu) for his support throughout my period of PROJECT, I remain indebted. I would like to extent my gratitude also to my sisters [ Halima, Binta, Aziza, Zainab, Munnira and Maryam] for their support both advisably and prayerful right from my childhood. My greatest gratitude also goes to my supervisor in person of Dr. I.A Fulatan for his untiring efforts, in spite of his ever crowded programmes not only to guide me through the task of writing this project but also to read and correct my scripts. I would like to acknowledge the input of all staff of the Department of Mathematics, Ahmadu Bello University Zaria, Thank you for answering my questions and for all the corrections on how to handle situations. Also thanks to my friends (especially my Tutor Nauwasi Usman) and colleagues at the Department of Mathematics, Ahmadu Bello University Zaria, for making my stay in the department memorable.










ABSTRACT
The Project work discussed the application of Runge-Kutta methods as a means of solving ordinary differential equation is demonstrated with the help of a specific fluid flow problem. The numerical results obtained are compared with the numerical results obtained by implicit and explicit finite difference methods. The error analysis and computational efficiency analysis are performed to test each method, Limitations of the Runge- Kutta method are also given.












Table of Contents
DECLARATION ii
CERTIFICATION iii
DEDICATION iv
ACKNOWLEDGEMENTS v
ABSTRACT vi
CHAPTER ONE 1
1.1 INTRODUCTION 1
1.2 SCOPE AND LIMITATION 1
1.3      DEFINITIONS OF MAJOR TERMS 2
1.3.1 Differential Equations: 2
1.3.2 Ordinary Differential Equation: 2
1.3.3      A solution of Ordinary Differential Equation: 3
1.3.4 Order of Differential Equation: 3
1.3.5 Degree of Differential Equation: 3
1.3.6 Initial Value Problem: 3
1.4 Errors in Numerical Computations 4
1.4.1 Types of Errors 4
CHAPTER TWO 5
LITERATURE REVIEW 5
2.1 INTRODUCTION 5
2.1.1 Newton Leibniz Years 5
2.1.2 Bernoulli Years 7
2.2 Runge-Kutta Method of solving ordinary differential equation 7
CHAPTER THREE 10
METHODOLOGY 10
3.1 INTRODUCTION 10
3.2 Method of Successive Approximation 10
3.2.1 Picard’s Method 11
3.2.2 Taylors Series Method 12
3.3 Single One Step Methods 13
3.3.1 Euler’s Method 13
3.3.2 Improved Euler’s method 15
3.3.3 Modified Euler’s Method 16
3.4 Multiple-step method 16
3.4.1 Milne-Thompson method 17
3.4.2 Adams-Bash Forth method 17
3.5 Step- Size Control in The Runge-Kutta Algorithm 17
3.6     Stability in the Runge - Kutta Algorithm 18
3.7 Third order Runge-Kutta method (Rk3) 18
3.8 Classic Runge-Kutta Method (Rk4) 21
3.9 Sixth Order Runge Kutta Method (Rk6) 26
3.10 Systems of ordinary Differential Equation 27
3.11 Second Order Initial Value Problem 29
CHAPTER FOUR 35
ANALYSIS AND DISCUSSION OF RESULTS 35
4.1 DISCUSSION 35
CHAPTER FIVE 44
SUMMARY, CONCLUSION AND RECOMMENDATION 44
5.1 SUMMARY 44
5.2 CONCLUSION 44
5.3 RECOMMENDATION 45
REFERENCES 47

 CHAPTER ONE
INTRODUCTION
This research work has its main objective aimed to outline solutions of ordinary differential Equations using one of the numerical methods, the Runge-Kutta methods. The Taylor's series will be used to solve the problems initially before applying the Runge-Kutta methods of Solutions.
As said above, Runge-Kutta method is a numerical solution method, and the development of numerical methods forms the subject of numerical analysis, a. branch of mathematics which has assumed major importance with many active research workers.
The Runge-Kutta methods of solutions are in different orders, but this research work will only be limited to the second, third and fourth orders, due to the cumbersome nature of these orders. In many practical problems involving differential equation a table of value y = y (x) satisfy a given initial conditions which is often wanted for a limited range of value of X. For instance we may want values of y(x) when x = x0 X h, v0 = 3h, x0 + 4h … … … , x0 + nh, when h = 0.5, 011..etc
Even when a solution y = y (x) of a differential equation can be written in terms of elementary functions it is easier to obtain this limited table of value by numerical method especially when the solution is unpleasant one. [Atkinson, (1989)].
1.2 SCOPE AND LIMITATION
The classical Runge-Kutta method is the most accurate numerical procedure used in obtaining approximation solution of initial value problem. And the greatest disadvantage seems to be that is rather difficult to estimate the error and further the method does not offer easy choking possibilities.
1.3      DEFINITIONS OF MAJOR TERMS
1.3.1 Differential Equations:
A differential equation is any equation containing differential coefficients, it involves derivatives (or differentials) of one or more dependent variables y1, y2….,yk  with respect to one or more independent variables x1, x2 ..., xk. The differential equation also involves derivatives of some unknown junctions.
For example:

1.3.2 Ordinary Differential Equation:
The ordinary differential equation is any function that involves only one independent variable that is only ordinary differential coefficient. It is a function of x, y and the derivatives of y up to any order such as f[x, y, dy, d3 y] = 0 for y (the dependent variable) in terms of x (The independent variable}-For example:
d2y/dx2 + 2dy/dx + y = 0                                                                             (1.5)
d2y/dx2 + 3xy = ℮x                                                                                                                                 (1.6)
1.3.3      A solution of Ordinary Differential Equation:
The solution of Ordinary Differential Equation y(n) = F(x, y,y1...,yn-1) is a function y(x) defined over a sub interval JCI which satisfies the above equation identically over the interval J. Any solution y(x) of the above equation should have derivatives of y at least up to order n in the internal J, for every x in J, the point (x, y(x), y'(x), ..., y(n-1)(x)) should live in the domain of definition of the function F.
1.3.4 Order of Differential Equation:
This is the order of the highest derivative which occurs in the equation in question. For Example:
y1 + x y = 3.  This is of order one or first order.
1.3.5 Degree of Differential Equation:
This is the exponent of the highest power of the highest-ordered derivative after the equation has been cleared of fractions and radicals in the dependent variables of its derivatives.
For example:

Is an equation is of order three while the degree is two.
1.3.6 Initial Value Problem:
The initial condition of any differential equation is known as the initial value problem.
1.4 Errors in Numerical Computations
An error is a derivative from accuracy or correctness. In numerical analysis, an error (E) is the difference between the true value(T) and the approximation value(A) of a mathematical problem i.e T-A, where T = True value and A = Approximation value. [Helfrick, (2005)].
1.4.1 Types of Errors
There are two main sources of the total in numerical approximations:
The global truncation error arises from the cumulative effect of six causes:
At each step we use an approximation formula to determine yn+1 (leading to a local truncation error).
The input data at each step are only approximation correct since in general (tn)yn.
Round-off error: Also cumulative, arises from using only finite number of digits. It can be shown that the global truncation error forth Euler method is proportional to h2 and for the Runge- Kutta method is proportional h4 respectively.
Initial errors: These are errors in the data collection e.g when a data is obtained from a physical or chemical experiment.
Absolute error: Is the error between two values and define as
Eabs = |x-u| where x donates exact value and x denotes  its approximation.
Relative error: Is the Absolute Error/ Actual Value
Percentage error = Relative x 100
CHAPTER TWO
LITERATURE REVIEW
2.1 INTRODUCTION
A widely used approach in the time integration of initial value problems for time dependent partial differential equations (PDEs) is the method of lines. This method transforms the PDE into a system of ordinary differential equations (ODEs) by discretization of the space variables and the uses and ODE solver for the time integration. Since ODEs originating spaced discretized PDEs have a special structure, ODE is appropriate. For example the well known fourth order Runge-Kutta method is highly inefficient if the PDE is parabolic, but it performs often quite satisfactory if the PDE is hyperbolic. [Atkinson, (1989)].
In Mathematics, differential equation traces the development of “differential equation” from Calentus, which itself was independently invented by English physicist Isaac Newton and German mathematician Gottfried Leibniza. The  differential equations, in concise form, from synopsis of the recent article “The history of differential equation, 1970-1950 “Differential equation began with Leibniz, the Bernouli brothers and others from the 1980s not long after Newton’s “fluxional equation” in the 1970s”
Differential equation differs from ordinary equation of mathematics, in addition to variables and constants they also contain derivatives of one or more of the variables involved.
2.1.1 Newton Leibniz Years
The exact chronological origin of the subject of differential equation is a bit of murky subject, for what seems to be a number of reason: One being secretiveness, two being private publication issues (private works published only decades letter), and three being the nature of the battle of mathematical and scientific discovery, which is a type of intellectual “war” (In the words of English Thomas Yong).
In circa 1971, English physicist Isaac Newton wrote his then unpublished the method of fluxion and infinite series (published in 1736), in which he classified first order differential equations, known to him as fluxional equations in to three (3) classes (using modern notation)
ODE
Class 1 Class 2 P.D.E
+ 
     
The first two classes contain ordinary derivatives of one or more dependent variables, with respect to a single independent variable and are known today as “ordinary differentials equation”, the third class involves the partial derivatives of one dependent variable and today are called “partial differential equations”.
In 1976 Newton solve his first differential equation. That same year, Leibniz introduced the term “differential equation” (aerations differential is, latin) or to denote a relationship between the differentials dx and dy of two variable x and y.
In 1963, Leibniz solved his first differential equation and that same year Newton published the results of previous differentials equations solution method a year that is said to mark the inception for the differential equation as district field in Mathematics
2.1.2 Bernoulli Years
Swiss mathematician, brothers Jacob Bernoulli (1654-1705) and Johann Benoulli (1667-1748), in Besel, Switzerland were among the first interceptors of Leibniz version of differentials calculus, to disprove Newton’s principle, on account that the brothers could not accept the theory, which Newton had proven, that the event and the planets rotate around the sum in elliptical orbits [3] the first work on subject of differential equation, supposedly was itallain mathematician Gabriele manfredi’s 1707 on the construction of first differential equation written between 1701 and 1707, published in Latin [4] the book was largely based or the med on the views of the Leibniz and the Bernoulli brothers. Most of the publication on differential equation and partial differential equation, in the years to follow. In the 18th century, seemed to expand on the version developed Leibniz a methodology, employed by those as Leonthard Euler, Daniel Bernoulli, Joseph Lagrange and Pierre Laplace.
A differential equation that codifies the condition of an exact differential function on this type, the prime examples (according to standard model) being state function, such as entropy ds, enthalpy dh, energy du etc, are differential function that are said to be path independent (in the context of a change of state of a body quantified by the Cycle integral, symbol ø)
2.2 Runge-Kutta Method of solving ordinary differential equation
In numerical analysis, the Runge Kutta methods are an important family of implicit and explicit iterative methods, which are used in temporal discretization for the approximation of solutions of ordinary differential equations. One hundred years ago the German mathematicians C. Runge-Kutta and M.W. Kutta developed three methods around 1900, the first of these three methods is the mid point rule adapted to ordinary differential equations while the second and third methods are different versions of the trapezoidal rule. The last of these methods suggests interactive computation of the stage values.
This work published in 1895, extended the approximation method of earlier to a more elaborate scheme which was capable of greater accuracy. The idea was to propagate the solution of an initial value problem toward a sequence of small-stops. In each stop, the rate of change of the solution is treated as constant and is found from the formula for the derivative evaluated at the beginning of the step.
For the equation y'(x) = f (x, y(x)), with given initial value y (x0) = y0, the first step is from the initial x0 to a slightly larger value x1, say. The approximate solution at this point is taken to be y1 = y0 + (x1 − x0) f (x0, y0). In general for a sequence of time values solution approximations y0, y1, y2… are given by
yn = yn-1 + (xn-xn-1) f (xn-1, yn-1)
The idea of Runge-Kutta was to base the approximate solution, not on this unsymmetrical and inaccurate Riemann rule, but on such improved formulas as the midpoint and trapezoidal rules. The requirement of evaluating the derivative at the midpoint and point or a step not yet completed, was achieved by first performing an enter type of calculation to obtain preliminary approximation to the solution at one of these points.
This characteristics feature of Runge – Kutta methods or evaluating the function  f  a number of times in working the way through each stop, was thus established,  the early work of Runge-Kutta and his immediate successors Heun and Kutta, laying the foundations for these methods and developing practical methods of increasing accuracy and efficiency.
Runge-Kutta methods between 1985-1925. He explored three main schemes; which are expressed as in his terminology together with the modern nation
The first of these three methods is the midpoint rule adapted to ordinary differential equations while the second and third methods are different versions of the trapezoidal rule. The last of these methods suggests interactive computation of the stage values.
However, more anural today would be the method since this hints at the implicit trapezoidal rule method in 1900, K. Heun took the other conditions as far as 4 and introduced amongst other methods the following of third order
The paper by W. Kutta which appeared in 1901, took the analysis of Runge-Kutta methods as far as order 5. He made a complete classification of order 4 methods and introduced the famous method in the work on 5th order methods. The first phase in the Runge-Kutta methods ended in the work of E.J. Mystrom. He took the analysis of the fifth order methods to its completion, but more importantly, he extended the use of Runge-Kutta method to second order differential equation systems. These systematize in dynamical problems and often be solved more efficiently when posed in their original form rather than as converted to an equivalent first order system, consider the second order system below




CHAPTER THREE

METHODOLOGY
3.1 INTRODUCTION
In the chapter, effort has been made to explain the method of solution for differential equation involving function of one variable. In this chapter extension will be made to a system (a collection) of differential equation consisting of n-equation in n-unknowns. Also this chapter examines Runge-Kutta method. RK4 for solving system of differential equations in each of this, we shall restrict ourselves to a system of first order and second order differential equation (initial value problem).
In general method that requires only knowledge of yn to determine yn+1 are called starting or single step method while that which make use of data at more than one point say yn, yn-1,yn-2 to determine yn+1 are called continuing or multi-step method. Numerical method falls into the following three classes. [Hairer, (1993)].
Method of successive approximation
Simple one step method.
Multi-step method
3.2 Method of Successive Approximation
This is a class method in which the unknown function; y is expressed in terms of X from which value of y can be obtained by direct substitution. Unclear its method, we may have the following method; i.e.
3.2.1 Picard’s Method
This method after Emily Picard a German Mathematician (who died in the early 19th Century) assume that the solution of the initial value problem. f (x, y), y (xo) = yo can be expressed as a power series of x. Picard’s method is a method that converts the differential equation into an equation involving integrals. Some differentials equation are difficult to solve but Picard’s method provides a numerical process by which a solution can be approximated i.e. if
 = f (x, y)                                                                                            (3.1)
Integrating equation (3.1)
            = ) dx
Therefore (3.1) becomes y(x) = yo+dx                                  (3.2)
The equation (3.2) which is an integral is called Picard’s function integrating
 y2(x) = y0−1 dx
yn (x) = y0 + 
This function y1(x), y2 (x), y3 (x)... yn(λ) are called Picard’s iterates converges.
Example: use Picard’s method to solve IVP dy/dx +2y = x,y(0) =1, h = 0.1
Solution
y0 = 0, f (x,y ) = dy/dx = x- 2y
y1(x) = y0 
= 1+ y0 
= 1-2x+ x2/2
y2(x) = y0
1 +dx
= 1 – 2x +   -
y3(x) = yo+
1-2x +  -+
yn (x) = y0 +
y1 (0.1) = 1- 2(0.1) – (0.1)2/2 =0.6973
y2 (0.2) = 1 – 2(0.2) - 5(0.2)2/2-(0.2)2/3 = 0.6973
y3 (0.3) = 1 – 2(0.3) + 5(0.2)2/2- 5(0.3)3/3 + (0.3)3/6 = 0.5881
And to calculate the error the exact value (or True value T) to the D.E
dy/dx + 2y = x y(0) = 1
y(x) = 5/4℮-2x + x/2 – 1/4
Which implies that y (0.1) = 0.8234, y(0.2) = 0.687 and y(0.3) = 0.5860
Therefore Absolute Error
Eabs = 0.5860-0.5881= 0.0021
Relative Error =0.0021/0.3860= 0.0035
3.2.2 Taylors Series Method
This method in named after English Mathematical; Brook Taylor (1685-1993) and the method says is dy/dx = f (x,y), y(xo) = yo is an IVP and if f(x,y) is analytic at x = xo then the solution of the equation can be expressed in its form;
y(x) = y(xo) + (x-xo)y' (xo) + (x-xo)2/2! y"(xo) … … … + (x-x0)n yn(x0)/n!
En+1 where En+1  ≤ │(x-xo)n+1yn+1(xo)/(n+1)!│, Eo ϵ (x, xo)
The solution of example can expressed as
y(x) = 1 + x y'(xo) + x2 y" yn+1(xo)/2! … … + xn yn/n!(xo) + En+1
This implies y (x) – 1 – 2x 5x2/2 -5x3/3 +5x4/4 … … … y(0.1) = 0.8234
y(0.2) = 0.6346 y(0.3) = 0.5867
3.3 Single One Step Methods
In this method, each successive value(s) yn+1 is obtained only from information from the immediate proceeding value(s) yn, some of the single step method include
Euler’s method
Improve Euler’s method
Modified Euler’s method
Runge-Kutta method
3.3.1 Euler’s Method
The Euler’s method named after Leonhard Euler is the simplest of one-step method and has a limited application because of its low accuracy.
The Euler’s method states that id dy/dx = f(x,y).
y(x0) = y0 is an IUP, with step size h, and if y1, y2, … … … yn corresponds to the value of x = xo,x1,x2 … … … xn where xn = x0 + xn
Then y1= y0+ hf (x0,y1)
y2 = y1+ hf (x1, y1)
yn = yn-1 h f (xn-1,yn-1)
Proof
Given the IVP dy/dx = f (x,y), y(xo) = y0
From the Taylor’s series expansion of y(x) about x= xo;(xo,h) + hy'(x), neglecting 2nd order and higher order, y'(x) = y(xo+h)- yo/h
y(xo+h) - y(xo)/h = f(x,y)
This implies that y(xo + h) = h(x0) + hf(x0,y0)
y1 = y0 +hf(x0, y0)
yn-1 = yn+ hf (xn,yn)
En+1 = │(x-x0)n+1yn+1(x0)/(n+1)!│, Eoϵ(x,x0)
Note that the error margin due to truncation of the Taylors series is relatively higher than in some other method as will be seen later in this section. Example:
Use Euler’s method to solve the IVP y1 + 2y = x, y(0), h = 0.1
Solution
y1 = f (x, y) = x-2y, y0 = 1, x0 = 0, h = 0.1
Now
y(0.1) = y0+hf(x0, y0)
y (0.1) = 1+0.1 (-2) = 0.8000
y (0.2) = 0.8+0.1(-1.5) = 0.6500
Also y (0.3) = 0.6 + 0.1 (-1.1) = 0.5400
Similarly for yn , y3,………. and so on to calculate the error. Exact value =0.5860
Eabs = 0.5860 - 0.5400 = 0.046
And relative Error = 0.046/0.5860 = 0.0784
3.3.2 Improved Euler’s method
The numerical method defined by the formula;
 yn+1 = yn+h/2 (f (x0, yn) - f (xn-1, yn-1))
Where yn+1 = yn+ hf (x0, yn)
This is commonly known as the improved Euler’s method or Heur is method.
Example
dy/dx = x-2y, y (0 )= 1  using improved Euler’s method to solve h = 0.1
Solution
f  (x, y) = x - 2y, yo = 1, x0 = 0, h = 0.1 for x1=0
y' (0.1) =1+0.1(-2) = 0.8000
y (0.1) = 1+0.05(-3.5) = 0.8250
y' (0.2) = 0.825 + 0.1(-1.55) = 0.6700
y (0.2) = 0.825 + 0.005(-2.69) = 0.6905
y' (0.3) = 0.6905 + 0.1(-1.181) = 0.5724
y (0.3) = 0.6905 + 0.05(-2.0258) = 0.5892
Similar for y4, y5...............and so on
3.3.3 Modified Euler’s Method
Another modification of the Euler’s method is to use the slope of the function at the estimated midpoint of (xn, yn) and (xn+1, yn+1) to approximate yn+1. Thus yn+1.=  yn+ fh (xn+h/2+h/2f (xn, yn))  this method is  called modified Euler’s method (or polygon method).
Example;
Using modified Euler’s method to solve the IVP
dy/dx = x-2y, y(0)=1 and h=0.1
Solution
y (0.1) = 1+0.1 (-1.75) = 0.8250
y(0.2) = 0.825+0.1(-3.45) = 0.6905
y(0.3) = 0.6905+0.1(-1.0129) = 0.5892
3.4 Multiple-step method
In this class of methods, the proceeding (next) value of a function is computed using four values of the function from different (previous) value such that yn, yn+1 from two previous steps is needed. These methods include;
Milne; Thompson method
Adams; Bash forth method
The above methods are called predictor formula which is used to predict (estimate) the required value of the function and the other formula (corrector) improve on the accuracy of the value obtained from the corrector formula.
Furthermore, the predictor-corrector formula can be derived from the numerical integration formula using Newton’s backward or forward  difference formula. This formula can be derived from the Milne-integration scheme and is represented by a pair of equation.
3.4.1 Milne-Thompson method
Milne-Thompson predictor corrector method is given as
Yn-1p = yn-3 + (2yn-21-yn-11-2yn1)
Corrector: yn-3 + (yn-11-4yn1-yn-11)
And starting values y1, y2 and y3 are obtained from the one single-step methods.
3.4.2 Adams-Bash Forth method
Adams bash forth predictor-corrector method is given as predictor :
Predictor:
Yn+1p-yn+(55yn-1-39yn-139yn-2+37n-2-9n-3)
Corrector:yn+1p=yn+(9y1n-1-19y1n-15y1n-1-y1n-2)anbe
3.5 Step- Size Control in The Runge-Kutta Algorithm
In order to choose a reasonable step size, we need some estimate of the error being committed in integrating across one-step. On the one hand, the step size should be small enough to achieve required accuracy (if possible), on the other hand, it should be as large as possible in order to keep rounding errors (a function of the number of arithmetic operations performed) under control and to avoid an excessive number of derivative evaluations. [Vetterling, (2007)].
3.6     Stability in the Runge - Kutta Algorithm
Stability is said to be another criterion for selecting an algorithm for the solution of a differential equation with given initial conditions.
This is a somewhat ambiguous term and appears in the literature with a variety of qualifying adjectives (inherent, partial, relative, weak, strong, absolute etc.). In general, a solution is said to be unstable if errors introduced at some stage in the calculations (for example, from erroneous initial conditions or local truncation or round-off errors) are propagated without bound throughout subsequent calculations.
Inherent instability is associated with the equation being solved and the initial condition specified, but does not depend on the particular algorithm being used. Depending on the equation being solved, its initial conditions, and the particular one step method being used, inherently unstable another form of instability, Partial instability, may be observed, even when the equation is not. [Dahlquist, (1963)].
3.7 Third order Runge-Kutta method (Rk3)
The method was developed by Carl-Runge (1856-1927) and M .W Kutta (1876-1944). In this method is an = IVP and if
 and




Contributing this  to  gives



The scheme is called Runge-Kutta of order three. It is more accurate than the improved Euler. In calculating a numerical solution to initial value problem using third-order Runge-Kutta method, one first computes the various  and then determiney n+!. Since each  depends on  the set of three  values must newly compute each time a new  is encountered.
Example;
Solve the
Solution:
For



Nown+1



For



Now;



For




Now;



Similarity, the same procedure to obtain  and so on.
3.8 Classic Runge-Kutta Method (Rk4)
There is different order of Runge-Kutta method in which the fourth-order Runge-Kutta method also called classical Runge-Kutta method.
This is one of the most widely used method and it is particularly suitable in cases when computation of higher derivatives is complicated. This method like the previous method can also be applied to higher order up by means of transforming the system to first order equation. [Tan and Chen, (2012)].
In this case the  can be successively as follows;

Where;




Derivation of RK4
Let  be a deferential equation with starting point  and interval length h.
Then we put;





Constants are determined in such a way that  becomes a good approximation of  as possible. By series expansion the following system is obtained.



)=



We have eight equations with ten unknowns, chosen two equation to be arbitrary; we assume ,       
And if we further chose  and hence

Thus, we have the final formula system





Example;
Solve the  using classical Runge-Kutta method (RK4).
Solution





Having the equation  with initial condition




(0.1,0.8225)=0.1545

Now the initial for the next computation are;






The initial for the next computation are








Hence the same procedure for y1, y2, y3, … … …, yn
3.9 Sixth Order Runge Kutta Method (Rk6)
We shall show only two example of higher order Runge-Kutta method. These may be of use in finding higher accuracy starting values for linear multistep algorithms. [Kaw and Kalu (2008)]. The example is the sixth order Kutta Nystrom method it is a six-stage method as:
;









3.10 Systems of ordinary Differential Equation
In this section we are working on the application of the 4th order R-K method systems of equations. Mathematical models of many applications in science engineering and social sciences such as problems related to financial management. Involves a system several 1st order equations. These system equations can be represented as follows:
=ym(
The equation should be solved for  over interval (a,b).
Consider a system of differential equations in two unknown functions  are functions of the independent variable. 
That is;
=








And where;




Example;
Solve the system of equation

Solution:


Using the formula above, we have the following










3.11 Second Order Initial Value Problem
An important property of system study is that the second order initial value problem;






Enhanced an alternative way to solve 2nd order IVP. A system of two (2) differential actions are obtained by substitution.
Using Runge-Kutta Method
Runge-Kutta method can be applied directly to differential equations higher order. For simple the equation  we put  and obtain the following system of first equation.


This is a special case of


Which integrate through

2
k3 = h f (x +1/2 h, y + ½ k2, 2 + ½ l2);
k4 =h f (x +h , y + k3,2 + l3); l4 =hg(x +h, y + k3, 2 +l3)
k =1/6(k1 +2k2 + 2k3 +k4); l = 1/6 (l1 +2/2 + 213 + l4)
The new value are ( x +h , y + k ,2 + 1 ).Which are used for the next computation
Example;
Solve the second order initial value problem; y" = x y' –y
y' = 0,y(o) = 1 with step length, h = 0.1. using fourth order Runge -Kutta method
Solution;
Given the equation y" = x y' – y
Consider y' = z = 0 (= f  (x ,y, z)),and y" = z' =x2 – y (= g (x, y, 2)).
Tanking h = 0.1 with initial condition: x = 0, y = 1, y = 0 (= 2).using the formula;
k1 = h f  ( x, yz ); l1 = h g ( x, y, 2)
k2  = h f  ( x + ½ h, y + ½ k1, 2 + ½  l1 ); l2 = h g ( x + ½ h, y + ½ k1 z + ½ l1 )
k3 = h f ( x + ½ h ,y + ½ k2, 2 + ½ l2); l3 =h g ( x +1/2 h, y + ½ k2, 2 + ½ l2)
k4 = h f ( x + h, y + k3, 2 + l3); l4 = h g ( x +h, y + k3, 2 + l3)
k = 1/6 ( k1 + 2k2 +2k3 +k4); l = 1/6 ( l2 + 212 +213 + l4)
Now
k1 = h f  (o,1,0 )=0
l1 = h g ( 0, 1, 0 ) = - 0.1000
k2 = h f ( 0.05,1,-0.05 ) = 0.0050
l2 = h f ( 0,05,1,-0.05) = -0.1002
k3= h f ( 0,15,0.9875,-0.15070) = -0.1010
k4 = h f  ( 0.2,0.9875,-0.2011) = - 0.0201
l = 1/6 (-0.6069)
and
k = 1/6 (-0.0901)= 0.0150
Then ( x +h, y + k, 2 + 1 ) ≈ ( 0.2,0.9800,-0.2012) are the initial value for the next step
 i.e x2 = 0.2
y (0,2) = -0.0642
y'(0.2) = z (0.1) = -0.2012
For compute the next step
k1 = h f  (0.2,0.9800,-0.2012) = -0.0201
k1 =h g (0.2,0.9800, - 0.2012) = -0.1020
k2 =h f  (0.25,0.9699, -0.3044) = -0.0252
l2 =h g ( 0.25,0.699, - 0.2522) = -0.1032
k3 =h f (0.25,0.9674, -0.3044) = -0.0304
l3 = h g(0.25,0.9674, -0.3044) = -0.1043
l4 = h f (0.3,0.9496, -0.3055) = -0.0305
l4 = h g (0.3,0.9496, -0.3055) = -0.1041
Then x + h, y +k, z-1)≈(0.3,0.9229,-0.3047) are the initial value for the next computation
i.e  x1 = 0.3
y (0.3) =0.9229
And y' (0,3) = -0.3047
Hence, the same procedure for y4 y5 -- -- -- -- -- -- and y (0,4), y (0,5)… … … … … …
Example (Using Euler’s Method)
Solve the second order initial value problem; y" = x y' –y
y' = 0,y (0) = 1 with step length h = 0.1 using Euler’s method
Solution
Given the equation y" = x y' – y
Consider y' = z = 0 (= f(x, y, z)), and y" = z' = x z – y (=g(x, y, z)).
Taking h = 0,1 with initial condition: x = 0, y = 1, y' = 0 (=z).using the formula:
yn+1 = y' n + h f  ( Xn', yn', 2n ) and
y'n+1  =  y'n + h f  ( xn, yn', zn)
y1 = 1 + 0.1f ( 0,1,0 ) =1,000
And y'1 = 0 + 0.1g (0, 1, 0 ) = ─ 0.1000
Then
( x + h, y, y'1) ≈ (0.1,1, - 0.1 ) are the initial values for the next step
i,e x = 0.1000 and y' (0.1) = ─ 0,01000
Now to compute for the next step
y2  = 1 + 0,1 f (0, 1, 1 – 0.1) =0.9900
And
 y'2 = 0.1 +0.1g ( 0, 1, 1 – 0.1 ) = -2010
Then
( xn +  h, y2, y'n) ≈ (0.2000,0.9900 – 0.2010) there the initial values for the next step.
Also to compute for the next step
y3 = 0.9900 + 0.1 f  (0.2000,0.9900,-0.2010) = -0.3040
Then
( xn  + h, y2, y'3) ≈ ( 0.3000, 0.9699, ─ 0.3040) are the initial values for the next computation i.e
x3 = 0.3000
y3 = 0.9699
and
y'3 = ─ 0.3040
Hence, the same procedure for y(0.4), (0.5) … y' (0.4), yi (0.5) …






CHAPTER FOUR
ANALYSIS AND DISCUSSION OF RESULTS
4.1 DISCUSSION
The result of Numerical methods via Euler, improved ever and Runge-Kutta method discussed in the previous chapter are given in the table below. [Lambert, (1999)].
Table 4.1: From the result in the above table we can see that Runge-Kutta method of 4th  order is more accurate as a scheme for solving IVP.
N
X
Euler’s
Improved Euler’s
3rd order Rk method
4th  order Rk method
y-exact

0
0.0
1.0000
1.0000
1.0000
1.0000
1.0000

1
0.1
0.8250
0.8250
0.8233
0.8234
0.8234

2
0.2
0.6500
0.6905
0.6878
0.6897
0.6879

3
0.3
0.5400
0.5892
0.5848
0.5860
0.5860

The use of this formula to compute the approximations  y1, y2, y3 ………. Successively constitutes the Runge-Kutta method. Note that yn+1 = yn+h.k takes the Euler form.
 If we write     k = 1/6(k1+2k2+2k3+k4)
For the approximate average slope on the interval [xn, xn+1],
We first apply the Runge-Kutta method to the illustrative initial value problem
 dy/dx = x + y, y (0) = 1
The exact solution of this of this problem is y(x) = 2℮x − x─1. To make a point we use h=0.5, a l\arger step size. So only two steps are required to go from x = 0 to x = 1.
In the first step we use the formulas
K1 = 0 + 1 = 1,
K2 = (0 + 0.25) + (1+ (0.25). (1)) = 1.5,
K3 = (0+0.25) + (1+ (0.25). (1.5)) =1.625,
K4 = (0.5)+ (1+(0.5).(1.625)) = 2.3125,
And then
y1= 1+0.6/6[1+2.(1.5)+2.(1.625)+2.3125] = 1.7969
Similarly the second step yields y = 3.4347
Table 4.1: applying the improved Euler method with step size h = 0.1. We see that even with the larger step size, the Runge -Kutta method gives four to five times the accuracy (in terms of relative percentage errors) of the improved Euler method.




Table 4.2 :  Runge-Kutta and improved Euler results for initial value problem dy/dx   = x +y, y (0) =1.
              Improved Euler method
Runge-Kutta

X                 y with h= 0.1   %  error
actual y
y with h= 0.5       % error                           


0.0                  1.0000           0.00%
1.000
1.0000                0.00%                               


0.5                  1.7949            0.14%
1.7974
1.7969                0.03%                               


1.0                  3.4282             0.24%
3.4366
3.4347                 0.05%                               



Shows the result obtained by applying the improved Euler and Runge-Kutta method to the problem dy/dx = x + y, y (0)=1 with the same step size h = 0.1. The relative error in the improved Euler value at x= 1 is about 0.24%, but for the Runge-Kutta value, it is 0.00012%. In this comparison the Runge-Kutta method is about 2000 times as accurate, but requires only twice as many function evaluation, as the improved Euler method.


Table 4.3: Runge-Kutta and improved Euler results for the initial value problem dy/dx = x +y y(0)= 1, with the step size h = 0.1
X
Improved Euler y
Runge-Kutta y
Actual y


0.1
1.1100
1.110342
1.110342


0.2
1.2421
1.242805
1.242806


0.3
1.3985
1.399717
1.399718


0.4
1.5818
1.583648
1.583649


0.5
1.7449
1.797441
1.797443


0.6
2.0409
2.044236
2.044238


0.7
2.3231
2.327503
2.3277505


0.8
2.6456
2.651079
2.651082


0.9
3.0124
3.019203
3.019206


1.0
3.4282
3.436559
3.436564



For the Runge-Kutta method results in a rapid decrease in the magnitude of errors when the step size h is reduced (except for the possibility that very small step sizes may result in unacceptable round off errors). It follows from the inequality that (on a fixed bounded interval) halving the step size decreases the absolute error by a factor of (1/2)4 =1/16. Consequently, the common practice of successively halving the step size until the computed results “stabilize” is particularly effective with the Runge-Kutta method.
In example 5, we saw that Euler’s method is not adequate to approximate the solution y(x) of the initial value problem, dy/dx =  x2 + y2, y (0) = 1
as x approaches the infinite discontinuity near x = 0.969811. Now we apply the Runge-Kutta results on the interval [0.0, 0.9], computed with step sizes h = 0.1, h = 0.025. There is still some difficulty near x = 0.9, but it seems safe to conclude from these data that y(0.5) = 2.0670.
Table 4.4: Mathematical Models and Numerical Methods. Approximating the solution of the initial value problem.
X
y with h = 0.1
y with h = 0.05
y with h = 0.025

0.1
1.1115
1.1115
1.1115

0.3
1.4397
1.4397
1.4397

0.5
2.0670
2.0670
2.0670

0.7
3.6522
3.6529
3.6529

0.9
14.0218
14.2712
14.3021


 We therefore begin anew and apply the Runge-Kutta method to the initialvalue problem
dy/dx = x2+ y2, y(0.5) = 2.0670


Table 4.5: Approximating the solution of the initial value problem.
X
y with h = 0.01
y with h = 0.005
Y with h = 0.0025

0.5
2.0670
2.0670
2.0670

0.6
2.6440
2.6440
2.6440

0.7
3.6529
3.6529
3.6529

0.8
5.8486
5.8486
5.8486

0.9
14.3048
14.3049
14.3049


Table 4.5 Shows results on the interval [0.5, 0.9], obtained with step sizes h = 0.01, h = 0.005, and h = 0.0025. We now conclude that y(0.9) = 14.3049.
Finally, Table 4.6 Show results on the interval [0.90, 0.95] for the initial value problem
dy/dx = x2 + y2, y (0.9) = 14.3049,
Obtained using step sizes h = 0. 002, h = 0.0005. Our approximate result is y (0.95) = 50.4723. The actual value of the solution at x = 0.95 = 50471867. Our slight overestimate results mainly from the fact that the four-place initial value is (in effect) the result of rounding up the actual  value y(0.9) = 14.304864; such errors are magnified considebly as we approach the vertical asymptote.



Table 4.6: Approximate Consider the seemingly innocuous initial value problem
X
Y with h = 0.002
y with h = 0.001
Y with h = 0.0005

0.90
14.3049
14.3049
14.3049

0.91
16.7024
16.7024
16.7024

0.92
20.0617
20.0617
20.0617

0.93
25.1073
25.1073
25.1073

0.94
33.5363
33.5363
33.5363

0.95
50.4722
50.4723
50.4723


  dy/dx = 5y - 6е-x, y(0) = 1                                                                                        (1)
Whose exact solution is y(x) = е-x. The table below shows the result obtained by applying the Runge-Kutta method on the interval [0, 4] with step sizes h= 0.2, h= 0.1, h= 0.05. Obviously attempts are spectacularly unsuccessful. Although y(x) = е-x → 0 as x → +∞. It appears that our numerical approximation is headed toward-∞ rather than zero.






Table 4.7 : Attempts to solve numerically the initial value problem.

Runge-Kutta y
With h=0.2
Runge-Kutta y
With h= 0.1
Runge-Kutta y
With h=0.05
Actual y

0.4
0.8
1.2
1.6
2.0
2.4
2.8
3.2
3.6
4.0
0.66880
0.43713
0.21099
-0.46019
-4.72142
-35.53415
-261.25023
-1,916.69935
-14059.35494
-103.126.5270
0.670020
0.44833
0.29376
0.14697
-0.27026
-2.90419
-22.05352
-163.25077
-1205.71249
-8903.12866
0.67031
0.44926
0.19802
0.10668
-0.12102
-1.50367
-11.51868
-85.38156
-631.03934
0.67032
0.44933
0.30199
0.20190
0.13534
0.09072
0.06081
0.04076
0.02732
0.01832


The explanation lies in the fact that the general solution of the equation dy/dx  = е-x + cе5x     
The particular solution 1 satisfying the initial condition y(0)= 1is obtained with C= 0.  But any departure, however small, from the exact  general  solution y(x)= е-x even if due only to round off error introduces in effect a non zero  value of C in (2).
Difficulties of the sort illustrated by example above some times are unavoidable, but one can at least hope to recognize such a problem when it appears. Approximate values whose order of magnitude varies with changing step size are a common indicator of such instability. These difficulties are discussed in numerical analysis textbooks and are the subject of current research in the field. 















CHAPTER FIVE
SUMMARY, CONCLUSION AND RECOMMENDATIONS
5.1 SUMMARY
The basic definition of general terms such as differential equation, initial value problem, boundary value problem, homogeneous differential equation, linear and non-linear equations, errors and types of errors in Numerical computations were all discussed in chapter one.
Also numerical method to find approximation solution to ordinary differential equation such as Euler’s method. Improved Euler method and example were all discussed in chapter two of this booklet.
The main discussion of Runge-Kutta method, derivation and limitation of classical Runge-Kutta method, also RK4 for the system of differential equation solving were discussed in chapter three of this booklet However, RK6 (Nystrom method) was also mentioned.
The findings (results) of this project is represented in a clearly understandable form, tables, figures (graphs and diagrams etc) and are appropriately numbered and labeled.
5.2 CONCLUSION
 Runge-Kutta methods are among the numerical methods used in solving Ordinary Differential Equations.
In this project, Runge-Kutta methods of solving ordinary Differential Equation, the use of Runge-Kutta of orders two, three and four to obtain solutions for the initial value problems of ordinary differential equation of order one was used. The Euler’s method is the easiest to apply but required many steps (a small step size) to achieve any reasonable degree of accuracy, the improved Euler’s method is more accurate (with the same step size) but involved more arithmetic calculations at each iteration. Also, first order Taylor's series approximation was used to help us evaluate the solutions. From the result of the calculation, it was observed that the result obtained using six decimal places for the Runge-Kutta methods were almost similar.
In comparison, the accuracy of the result from the computation improves as the order of the methods used increases, thus making the fourth order Runge-Kutta method the best.
Finally, general study of these methods indicates that Runge-Kutta method of order n for n˂0 requires exactly n-functional evaluation n˃4 and more than n-evaluation if n˃4, this Runge-Kutta of order 4 usually gives optimal result (solution).
5.3 RECOMMENDATIONS
This method discusses the numerical method of solving ordinary differential equation. In particular method of successive approximation. Single (one-step) method and multi-step method. it also discusses the scholars who contributed and formulated some of the method. By numerical method for solving initial value problems, we meant a procedure for finding approximation value y0, y1, y2, ………..yn of exact solution y at the point x1, x2……..xn. The first step is to estimate y1 from the initial condition given i.e. y1 = F (x0). After finding the value of y1 we determine y2 and so on,
In general method that requires only knowledge of yn to determine yn+1 are called starting or single step method while that which make use of data at more than one point say yn, yn-1,yn-2 to determine yn+1 are called continuing or multi-step method. Stroud. K.A. (1996). Numerical method falls into the following three classes
the method of solution for differential equation involving function of one variable. In this chapter extension will be made to a system (a collection) of differential equation consisting of n-equation in n-unknowns. Also this chapter examines Runge-Kutta method. RK4 for solving system of differential equations in each of this, we shall restrict ourselves to a system of first order and second order differential equation (initial value problem).











                                               
REFERENCES
Atkinson, K.A. (1989). An Introduction to Numerical Analysis (second edition). John Wiley and Sons, New York, Pp; 1 - 4 & 6 -11.
Dahlquist, G. (1963). A special Stability Problem for Linear Multistep methods. Norton and Company, New York. Pp; 20.
Helfrick, A.D. (2005). Modern Electronic Instrumentation and Measurement Techniques.  Clarendon Press, Oxford. Pp; 4-5
Hairer, et al, (1993). Solving Ordinary Differential Equations I: Nonstiff problems. Springer Verlag, New York. Pp; 12-20
Kaw, A. and Kalu, E. (2008). Numerical Method with application (first edition).Troy and Sons, inc.  Pp; 29-37.
Lambert, J.D. (1991). Numerical solution of ordinary differential systems. The initial value problem. John Wiley and Sons, New York. Pp 38-45.
Tan, D. and Chen, Z. (2012). On a general formula of fourth order Runge- Kutta Method, Journal of Mathematics Sciences and Mathematics Education. Vol 1 (1), 78-95.
Vetterling, W.T. (2007). Adaptive Stepsize Control for Runge-Kutta, Cambridge University Press. Pp; 20