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

给出以下输入：

1. 一个普通的微分方程，以x和y的形式定义dy / dx的值。

   $\LARGE \frac{dy}{dx} = f(x, y)$

2. y的初始值，即y(0) 。

   $\LARGE y(0)= y_o$

任务是在给定的点x上找到未知函数y的值，即y(x) 。

例子：

Input: x₀ = 0, y₀ = 1, x = 2, h = 0.2
Output: y(x) = 0.645590

Input: x₀ = 2, y₀ = 1, x = 4, h = 0.4;
Output: y(x) = 4.122991

方法：
Runge-Kutta方法找到给定x的y近似值。使用Runge Kutta二阶方法只能求解一阶常微分方程。

下面是用来从以前的值是_否计算下一个值y _{n + 1个}公式。

所以：

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

计算y(n + 1)值的基本公式：
$\LARGE K_{1} = h*f(x_{n}, y_{n})$
$\LARGE K_{2} = h*f((x_{n} + \frac{h}{2}), (y_{n} + \frac{K_{1}*h}{2}))$
$\LARGE y_{n+1} = y_{n} + K_{2} + (h^{3})$

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

• K ₁是基于间隔开始时的斜率的增量(使用y)。
• K ₂是基于在该间隔的中点处的斜率的增加，使用(Y + H + K _1/2)。

该方法是二阶方法，这意味着局部截断误差为O(h ³ )量级，而总累积误差为O(h ⁴ )量级。

下面是上述方法的实现：

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

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

输出：

y(x) = 0.645590

参考： https : //en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods