Runge-Kutta 2nd order method to solve Differential equations //

Given the following inputs: 

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

The task is to find the value of unknown function y at a given point x, i.e. y(x).

Example:

Input: x0 = 0, y0 = 1, x = 2, h = 0.2 
Output: y(x) = 0.645590
Input: x0 = 2, y0 = 1, x = 4, h = 0.4; 
Output: y(x) = 4.122991   

Approach:

The Runge-Kutta method finds an approximate value of y for a given x. Only first-order ordinary differential equations can be solved by using the Runge-Kutta 2nd-order method.

Below is the formula used to compute the next value y_n+1 from the previous value y_n. Therefore: 

yn+1 = value of y at (x = n + 1)
yn = value of y at (x = n)
where
  0 ≤ n ≤ (x - x0)/h
  h is step height
  xn+1 = x0 + h

The essential formula to compute the value of y(n+1): 

K1 = h * f(x,y)
K2 = h * f(x/2, y/2) or K1/2
yn+1 = yn + K2 + (h3)

The formula basically computes the next value y_n+1 using current y_n plus the weighted average of two increments:  

• K₁ is the increment based on the slope at the beginning of the interval, using y.
• K₂ is the increment based on the slope at the midpoint of the interval, using (y + h*K₁/2).

The method is a second-order method, meaning that the local truncation error is in the order of O(h³), while the total accumulated error is order of O(h⁴).

Below is the implementation of the above approach: 

C++

// C++ program to implement Runge
// Kutta method

#include <bits/stdc++.h>

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;

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

        // Update next value of y
        y = y + (1.0 / 6.0) * (k1 + 2 * k2);

        // Update next value of x
        x0 = x0 + h;
    }

    return y;
}

// Driver Code
int main()
{
    float x0 = 0, y = 1, x = 2, h = 0.2;

    cout << fixed << setprecision(6)
         << "y(x) = " << rungeKutta(x0, y, x, h);
    return 0;
}

// This code is contributed by shivani

Java

// Java program to implement Runge
// Kutta method

class GFG {

    // A sample differential equation
    // "dy/dx = (x - y)/2"
    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;

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

            // Update next value of y
            y = y + (1.0 / 6.0) * (k1 + 2 * k2);

            // Update next value of x
            x0 = x0 + h;
        }

        return y;
    }

    // Driver Code
    public static void main(String[] args)
    {
        double x0 = 0, y = 1, x = 2, h = 0.2;

        System.out.println(rungeKutta(x0, y, x, h));
    }
}

// This code is contributed by Yash_R

Python3

# Python3 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 = round((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)

        # Update next value of y
        y = y + (1.0 / 6.0) * (k1 + 2 * k2)

        # Update next value of x
        x0 = x0 + h

    return y


# Driver Code
if __name__ == "__main__":

    x0 = 0
    y = 1
    x = 2
    h = 0.2

    print("y(x) =", rungeKutta(x0, y, x, h))

# This code is contributed by Yash_R

C

// C program to implement Runge
// Kutta method

#include <stdio.h>

// 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;

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

        // Update next value of y
        y = y + (1.0 / 6.0) * (k1 + 2 * k2);

        // Update next value of x
        x0 = x0 + h;
    }

    return y;
}

// Driver Code
int main()
{
    float x0 = 0, y = 1, x = 2, h = 0.2;

    printf("y(x) = %f", rungeKutta(x0, y, x, h));
    return 0;
}

C#

// C# program to implement Runge
// Kutta method
using System;

class GFG {

    // A sample differential equation
    // "dy/dx = (x - y)/2"
    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;

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

            // Update next value of y
            y = y + (1.0 / 6.0) * (k1 + 2 * k2);

            // Update next value of x
            x0 = x0 + h;
        }

        return y;
    }

    // Driver Code
    public static void Main(string[] args)
    {
        double x0 = 0, y = 1, x = 2, h = 0.2;

        Console.WriteLine(rungeKutta(x0, y, x, h));
    }
}

// This code is contributed by Yash_R

Javascript

<script>

// JavaScript program to implement Runge
// Kutta method

    // A sample differential equation
    // "dy/dx = (x - y)/2"
    function 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.
    function rungeKutta(x0, y0, x, h)
    {
        // Count number of iterations
        // using step size or
        // step height h
        let n = ((x - x0) / h);

        let k1, k2;

        // Iterate for number of iterations
        let y = y0;
        for (let 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);

            // Update next value of y
            y = y + (1.0 / 6.0) * (k1 + 2 * k2);

            // Update next value of x
            x0 = x0 + h;
        }

        return y;
    }

// Driver Code

    let x0 = 0, y = 1,
              x = 2, h = 0.2;

    document.write(rungeKutta(x0, y, x, h).toFixed(6));


</script>

Output:

y(x) = 0.645590

Time Complexity: O(n)
Auxiliary Space: O(1)