Runge-Kutta四阶方法求解微分方程

鉴于以下输入，

• 一个普通的微分方程，以x和y的形式定义dy / dx的值。
• y的初始值，即y(0)

因此，我们在下面给出。
$\frac{\mathrm{dx} }{\mathrm{d} y} = f(x, y),y(0)= y_o$
任务是在给定点x处找到未知函数y的值。
Runge-Kutta方法找到给定x的y的近似值。使用Runge Kutta四阶方法只能求解一阶常微分方程。
下面是用来从以前的值是_否计算下一个值y _{n + 1个}公式。 n的值为0、1、2、3，….(x – x0)/ h。这里h是台阶高度， x _{n + 1} = x ₀ + h
。步长越小意味着精度越高。
$K_1 = hf(x_n, y_n)\\ K_2 = hf(x_n+\frac{h}{2}, y_n+\frac{k_1}{2})\\ K_3 = hf(x_n+\frac{h}{2}, y_n+\frac{k_2}{2})\\ K_4 = hf(x_n+h, y_n+k_3)\\ y_{n+1} = y_n + k_1/6 + k_2/3 + k_3/3 + k_4/6 + O(h^{5})$

该公式基本上使用当前y _n加上四个增量的加权平均值来计算下一个值y _{n + 1} 。

• k ₁是基于间隔开始时的斜率的增量，使用y
• K ₂是使用Y + HK2_分之1基于在该间隔的中点处的斜率的增加，。
• k ₃还是使用y + hk _2/2基于中点处的斜率的增量。
• k ₄是使用y + hk ₃基于间隔结束时的斜率的增量。

该方法是四阶方法，这意味着局部截断误差为O(h ⁵ )量级，而总累积误差为O(h ⁴ )量级。
来源：https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

以下是上述公式的实现。

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 method
import 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

Python

# 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 method
x0 = 0
y = 1
x = 2
h = 0.2
print '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 method
using 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.

Javascript

 
// 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 = parseInt((x - x0) / h, 10);
 
    let k1, k2, k3, k4, k5;
 
    // 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);
        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 / 6) * (k1 + 2 * k2 +
                            2 * k3 + k4);;
 
        // Update next value of x
        x0 = x0 + h;
    }
    return y.toFixed(6);
}
 
// Driver code
let x0 = 0, y = 1, x = 2, h = 0.2;
 
document.write("The value of y at x is : " +
               rungeKutta(x0, y, x, h));   
                
// This code is contributed by divyesh072019
 

输出：

The value of y at x is : 1.103639

上述解决方案的时间复杂度为O(n)，其中n为(x-x0)/ h。
一些有用的资源，用于详细的示例和更多的解释。
http://w3.gazi.edu.tr/~balbasi/mws_gen_ode_txt_runge4th.pdf
https://www.youtube.com/watch?v=kUcc8vAgoQ0