Runge-Kutta 4th Order Method to Solve Differential Equation

Last Updated : 17 Jan, 2023

Given the following inputs,

• An ordinary differential equation that defines value of dy/dx in the form x and y.
• Initial value of y, i.e., y(0)

Thus we are given below.

The task is to find the value of the unknown function y at a given point x.
The Runge-Kutta method finds the approximate value of y for a given x. Only first-order ordinary differential equations can be solved by using the Runge Kutta 4th order method.
Below is the formula used to compute next value yn+1 from previous value yn. The value of n are 0, 1, 2, 3, ….(x – x0)/h. Here h is step height and xn+1 = x0 + h
. Lower step size means more accuracy.

The formula basically computes next value yn+1 using current yn plus weighted average of four increments.

• k1 is the increment based on the slope at the beginning of the interval, using y
• k2 is the increment based on the slope at the midpoint of the interval, using y + hk1/2.
• k3 is again the increment based on the slope at the midpoint, using  y + hk2/2.
• k4 is the increment based on the slope at the end of the interval, using y + hk3.

The method is a fourth-order method, meaning that the local truncation error is on the order of O(h5), while the total accumulated error is order O(h4).
Source: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

Below is the implementation for the above formula.

C++

 // C++ program of the above approach#include using namespace std; // A sample differential equation "dy/dx = (x - y)/2"float dydx(float x, float y){    return((x - y)/2);}  // Finds value of y for a given x using step size h// and initial value y0 at x0.float rungeKutta(float x0, float y0, float x, float h){    // Count number of iterations using step size or    // step height h    int n = (int)((x - x0) / h);      float k1, k2, k3, k4, k5;      // Iterate for number of iterations    float y = y0;    for (int i=1; i<=n; i++)    {        // Apply Runge Kutta Formulas to find        // next value of y        k1 = h*dydx(x0, y);        k2 = h*dydx(x0 + 0.5*h, y + 0.5*k1);        k3 = h*dydx(x0 + 0.5*h, y + 0.5*k2);        k4 = h*dydx(x0 + h, y + k3);          // Update next value of y        y = y + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4);;          // Update next value of x        x0 = x0 + h;    }      return y;} // Driver Codeint main(){    float x0 = 0, y = 1, x = 2, h = 0.2;    cout << "The value of y at x is : " <<            rungeKutta(x0, y, x, h);     return 0;} // This code is contributed by code_hunt.

C

 // C program to implement Runge Kutta method #include  // A sample differential equation "dy/dx = (x - y)/2" float dydx(float x, float y) {     return((x - y)/2); }  // Finds value of y for a given x using step size h // and initial value y0 at x0. float rungeKutta(float x0, float y0, float x, float h) {     // Count number of iterations using step size or     // step height h     int n = (int)((x - x0) / h);      float k1, k2, k3, k4, k5;      // Iterate for number of iterations     float y = y0;     for (int i=1; i<=n; i++)     {         // Apply Runge Kutta Formulas to find         // next value of y         k1 = h*dydx(x0, y);         k2 = h*dydx(x0 + 0.5*h, y + 0.5*k1);         k3 = h*dydx(x0 + 0.5*h, y + 0.5*k2);         k4 = h*dydx(x0 + h, y + k3);          // Update next value of y         y = y + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4);;          // Update next value of x         x0 = x0 + h;     }      return y; }  // Driver method int main() {     float x0 = 0, y = 1, x = 2, h = 0.2;     printf("\nThe value of y at x is : %f",             rungeKutta(x0, y, x, h));     return 0; }

Java

 // Java program to implement Runge Kutta methodimport java.io.*;class differential{    double dydx(double x, double y)    {        return ((x - y) / 2);    }         // Finds value of y for a given x using step size h    // and initial value y0 at x0.    double rungeKutta(double x0, double y0, double x, double h)    {        differential d1 = new differential();        // Count number of iterations using step size or        // step height h        int n = (int)((x - x0) / h);         double k1, k2, k3, k4, k5;         // Iterate for number of iterations        double y = y0;        for (int i = 1; i <= n; i++)         {            // Apply Runge Kutta Formulas to find            // next value of y            k1 = h * (d1.dydx(x0, y));            k2 = h * (d1.dydx(x0 + 0.5 * h, y + 0.5 * k1));            k3 = h * (d1.dydx(x0 + 0.5 * h, y + 0.5 * k2));            k4 = h * (d1.dydx(x0 + h, y + k3));             // Update next value of y            y = y + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);                         // Update next value of x            x0 = x0 + h;        }        return y;    }         public static void main(String args[])    {        differential d2 = new differential();        double x0 = 0, y = 1, x = 2, h = 0.2;                 System.out.println("\nThe value of y at x is : "                    + d2.rungeKutta(x0, y, x, h));    }} // This code is contributed by Prateek Bhindwar

Python3

 # Python program to implement Runge Kutta method# A sample differential equation "dy / dx = (x - y)/2"def dydx(x, y):    return ((x - y)/2) # Finds value of y for a given x using step size h# and initial value y0 at x0.def rungeKutta(x0, y0, x, h):    # Count number of iterations using step size or    # step height h    n = (int)((x - x0)/h)     # Iterate for number of iterations    y = y0    for i in range(1, n + 1):        "Apply Runge Kutta Formulas to find next value of y"        k1 = h * dydx(x0, y)        k2 = h * dydx(x0 + 0.5 * h, y + 0.5 * k1)        k3 = h * dydx(x0 + 0.5 * h, y + 0.5 * k2)        k4 = h * dydx(x0 + h, y + k3)         # Update next value of y        y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4)         # Update next value of x        x0 = x0 + h    return y # Driver methodx0 = 0y = 1x = 2h = 0.2print ('The value of y at x is:', rungeKutta(x0, y, x, h)) # This code is contributed by Prateek Bhindwar

C#

 // C# program to implement Runge// Kutta methodusing System; class GFG {         static double dydx(double x, double y)    {        return ((x - y) / 2);    }         // Finds value of y for a given x    // using step size h and initial    // value y0 at x0.    static double rungeKutta(double x0,                double y0, double x, double h)    {             // Count number of iterations using        // step size or step height h        int n = (int)((x - x0) / h);         double k1, k2, k3, k4;         // Iterate for number of iterations        double y = y0;                 for (int i = 1; i <= n; i++)         {                         // Apply Runge Kutta Formulas            // to find next value of y            k1 = h * (dydx(x0, y));                         k2 = h * (dydx(x0 + 0.5 * h,                             y + 0.5 * k1));                                          k3 = h * (dydx(x0 + 0.5 * h,                            y + 0.5 * k2));                                         k4 = h * (dydx(x0 + h, y + k3));             // Update next value of y            y = y + (1.0 / 6.0) * (k1 + 2                        * k2 + 2 * k3 + k4);                         // Update next value of x            x0 = x0 + h;        }                 return y;    }         // Driver code    public static void Main()    {                 double x0 = 0, y = 1, x = 2, h = 0.2;                 Console.WriteLine("\nThe value of y"                             + " at x is : "                 + rungeKutta(x0, y, x, h));    }} // This code is contributed by Sam007.

PHP

 

Javascript

 

Output
The value of y at x is : 1.10364

Time Complexity: O(n), where n is (x-x0)/h.
Auxiliary Space: O(1) as constant space for variables is being used

Some useful resources for detailed examples and more explanation.
http://w3.gazi.edu.tr/~balbasi/mws_gen_ode_txt_runge4th.pdf