Carl Runge
Martin Wilhelm Kutta
1 龙格-库塔四阶方法
数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
对于一阶精度的欧拉公式有：
yi+1=yi+h*K1 K1=f(xi,yi)
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式：
yi+1=yi+h*( K1+ K2)/2
K1=f(xi,yi)
K2=f(xi+h,yi+h*K1)
依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：
yi+1=yi+h*( K1+ 2*K2 +2*K3+ K4)/6
K1=f(xi,yi)
K2=f(xi+h/2,yi+h*K1/2)
K3=f(xi+h/2,yi+h*K2/2)
K4=f(xi+h,yi+h*K3)
通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头” 。
2 源代码
using System;
using System.Collections;
using System.Collections.Generic;namespace Legalsoft.Truffer.Algorithm
{/// <summary>/// 给定微分方程的一阶偏导方程/// </summary>/// <param name="x"></param>/// <param name="y"></param>/// <returns></returns>public delegate double SDE_Equation(double x, double y);/// <summary>/// 解微分方程的龙格-库塔四阶方法/// </summary>public static partial class Algorithm_Gallery{public static SDE_Equation dydx = null;/// <summary>/// Finds value of y for a given x/// using step size h/// and initial value y0 at x0./// </summary>/// <param name="x0">初值</param>/// <param name="y0">初值</param>/// <param name="x">求值点</param>/// <param name="h">步长</param>/// <returns></returns>public static double Runge_Kutta_4th_Order(double x0, double y0, double x, double h){int n = (int)((x - x0) / h);double y = y0;for (int i = 1; i <= n; i++){double k1 = h * (dydx(x0, y));double k2 = h * (dydx(x0 + 0.5 * h, y + 0.5 * k1));double k3 = h * (dydx(x0 + 0.5 * h, y + 0.5 * k2));double k4 = h * (dydx(x0 + h, y + k3));y = y + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);x0 = x0 + h;}return y;}}
}
3 源程序
using System;
using System.Collections;
using System.Collections.Generic;
namespace Legalsoft.Truffer.Algorithm
{
/// <summary>
/// 给定微分方程的一阶偏导方程
/// </summary>
/// <param name="x"></param>
/// <param name="y"></param>
/// <returns></returns>
public delegate double SDE_Equation(double x, double y);
/// <summary>
/// 解微分方程的龙格-库塔四阶方法
/// </summary>
public static partial class Algorithm_Gallery
{
public static SDE_Equation dydx = null;
/// <summary>
/// Finds value of y for a given x
/// using step size h
/// and initial value y0 at x0.
/// </summary>
/// <param name="x0">初值</param>
/// <param name="y0">初值</param>
/// <param name="x">求值点</param>
/// <param name="h">步长</param>
/// <returns></returns>
public static double Runge_Kutta_4th_Order(double x0, double y0, double x, double h)
{
int n = (int)((x - x0) / h);
double y = y0;
for (int i = 1; i <= n; i++)
{
double k1 = h * (dydx(x0, y));
double k2 = h * (dydx(x0 + 0.5 * h, y + 0.5 * k1));
double k3 = h * (dydx(x0 + 0.5 * h, y + 0.5 * k2));
double k4 = h * (dydx(x0 + h, y + k3));
y = y + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
x0 = x0 + h;
}
return y;
}
}
}
