Solving System of ODEs using Three Forms of RK4
by admin in Differential Equations , Math, Statistics, and Optimization , MATLAB Family on April 29, 2019This code solves multi ordinary differential equations (ODEs) using three forms of fourth order Runge Kutta (RK4). Two types of charts are presented the first is the numerical solution of the given equations and is about three charts each is for specific RK4 form. The second chart type presents the residuals.
Code Outputs:
 Three charts of the numerical solutions for three RK4 forms.
 Residual chart of the three RK4 forms.
 Final residual value of each RK4 form.
Final Residual Value of First Form = 0.042018
Final Residual Value of Second Form = 0.042018
Final Residual Value of Third Form = 0.042018
Input Requirements:
 ODEs you want to solve.
 Initial Conditions.
 Initial and Final Integration Interval.
 Number of Interval Divisions.
Editing the Code:
The code has been built in easy way to edit it. First add your equations in the Fun.m file the type the initial conditions in the main function and run the code.
You can uncomment any RK4 form, and present the other results.
About the Method:
Runge–Kutta methods are a family of implicit and explicit iterative methods, which include the wellknown routine called the Euler Method, used in temporal discretization for the approximate solutions of ordinary differential equations. These methods were developed around 1900 by the German mathematicians Carl Runge and Wilhelm Kutta.
List of Runge–Kutta methods:
which take the form
Each method listed on this page is defined by its Butcher tableau, which puts the coefficients of the method in a table as follows:
The most widely known member of the Runge–Kutta family is generally referred to as ‘RK4‘, ‘classical Runge–Kutta method‘ or simply as ‘the Runge–Kutta method‘.
Let an initial value problem be specified as follows:
Here y is an unknown function (scalar or vector) of time t, which we would like to approximate; we are told that y’, the rate at which y changes, is a function of t and of y itself. At the initial time t0 the corresponding y value is y0. The function f and the data t0, y0 are given.
Now pick a stepsize h > 0 and define
for n = 0, 1, 2, 3, …, using
 (Note: the above equations have different but equivalent definitions in different texts).
Here yn+1 is the RK4 approximation of y(tn+1), and the next value (yn+1) is determined by the present value (yn) plus the weighted average of four increments, where each increment is the product of the size of the interval, h, and an estimated slope specified by function f on the righthand side of the differential equation.
 k1 is the increment based on the slope at the beginning of the interval, using y (Euler’s method);
 k2 is the increment based on the slope at the midpoint of the interval, using y and k1;
 k3 is again the increment based on the slope at the midpoint, but now using y and k2;
 k4 is the increment based on the slope at the end of the interval, using y and k3.
In averaging the four increments, greater weight is given to the increments at the midpoint. If f is independent of y, so that the differential equation is equivalent to a simple integral, then RK4 is Simpson’s rule.
The RK4 method is a fourthorder method, meaning that the local truncation error is on the order of O(h^5), while the total accumulated error is on the order of O(h^4).
In many practical applications the function f is independent of t (so called autonomous system, or timeinvariant system, especially in physics), and their increments are not computed at all and not passed to function f , with only the final formula for tn+1 used.
Values of Butcher Tableau:
The first standard form of Runge Kutta has the following Butcher tableau values:

0 1/2 1/2 1/2 0 1/2 1 0 0 1 1/6 1/3 1/3 1/6
A slight variation of ‘the’ Runge–Kutta method is also due to Kutta in 1901 and is called the 3/8rule. The primary advantage this method has is that almost all of the error coefficients are smaller than in the popular method, but it requires slightly more FLOPs (floatingpoint operations) per time step. Its Butcher tableau is:

0 1/3 1/3 2/3 −1/3 1 1 1 −1 1 1/8 3/8 3/8 1/8
The third form of Runge Kutta has the following Butcher tableau values:

0 1/2 1/2 1/2 1/2 + 1/sqrt(2) 11/sqrt(2) 1 0 – 1/sqrt(2) 1 + 1/sqrt(2) 1/6 1/3 – 1/(3 sqrt(2)) 1/3 + 1/(3 sqrt(2)) 1/6
References:
[1] Runge, Carl David Tolmé (1895), ‘Über die numerische Auflösung von Differential gleichungen’, Mathematische Annalen, Springer, 46 (2): 167–178, doi:10.1007/BF01446807.
[2] Kutta, Martin (1901), Beitrag zur näherungweisen Integration totaler Differentialgleichungen.
[3] Ascher, Uri M.; Petzold, Linda R. (1998), Computer Methods for Ordinary Differential Equations and DifferentialAlgebraic Equations, Philadelphia: Society for Industrial and Applied Mathematics, ISBN9780898714128.
[4] Atkinson, Kendall A. (1989), An Introduction to Numerical Analysis (2nd ed.), New York: John Wiley & Sons, ISBN9780471500230.
[5] Butcher, John C. (May 1963), Coefficients for the study of RungeKutta integration processes, 3 (2), pp. 185–201, doi:10.1017/S1446788700027932.
[6] Butcher, John C. (1975), ‘A stability property of implicit RungeKutta methods’, BIT, 15: 358–361, doi:10.1007/bf01931672.
Share Now!