四阶龙格库塔（Runge-Kutta）求解微分方程-多种编程语言

✅作者简介：热爱科研的算法开发者，Python、Matlab项目可交流、沟通、学习。

🍎个人主页：算法工程师的学习日志

前期是分享了matlab下面实现四阶龙格库塔（Runge-Kutta）求解微分方程，这期分享一下C++、C、Java、Python下面的四阶龙格库塔（Runge-Kutta）求解微分方程。

前文传送门：​​matlab代码实现四阶龙格库塔求解微分方程​​

C++方法

#include <bits/stdc++.h>
using namespace std;


//  "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 Code
int 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;
}

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, 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));
  }
}

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 method
x0 = 0
y = 1
x = 2
h = 0.2
print ('The value of y at x is:', rungeKutta(x0, y, x, h))

The value of y at x is: 1.1036393232374955