龙格库塔(runge-kutta，RK)法求解微分方程_龙格库塔法解微分方程组

求解微分方程的意思就是，已知导函数，求原函数。

先声明一点，欧拉法、中值法、龙哥库塔法求解微分方程，得出的结果不是表达式，而是一系列离散点。

一、欧拉法递推

问题描述：已知y'=f(x,y)，求y(x)。

例题：已知y'=y，y(0)=1，求y(x)。

解：这个题目高数知识很容易求解，答案是y=e^x。但是这里我们不使用解析法，而使用迭代递推法去求y(x)。

把y(x)在x0处泰勒展开到一阶，有y(x0+h)=y(x0)+y'(x0)h，有了这个式子，我们就能根据y(0)递推出y(0+h)的值，进而推出y(0+h+h)、y(0+h+h+h)的值····。根据这一原理就能很容易的写出递推程序了：

clear;
euler(0.1);%减小步长可以提高运算精度
 
function euler(h)
%欧拉法解微分方程:y'=y     真实解为y=e^t
%形参：自变量的步长h
hold on;
idx = 1: 1: 10/h;
 
x(1)=0;
y(1)=1;
 
for i = idx
    x(i+1) = x(i) + h;
    dy = y(i);    
    y(i+1) = y(i) + h * dy;  %y1 = y0 + h * y0'
end
 
plot(x, exp(x), 'r');
plot(x, y, 'b');
legend('真实解','递推解');
end

运行结果如下：

我们可以发现，递推的精度有点差，不过下面我们还有更好的递推方法。

二、中值法递推

该方法是对欧拉法的改进，原理就是拉格朗日中值定理，y(x0+h)=y(x0)+h*y'(a)，其中a∈[x0, x0+h]，虽然a不是区间的中点，但是这里我们认为用[x0, x0+h]区间中点处x0+h/2的导数来递推，比单纯用左端点的导数做递推，精度更高。该方法所用的递推公式为：

中值法有两种，第二种是这样的：用区间左右端点处的导数的平均值来递推，本质上就是左右端点各预测一个增量△y，这两个增量权重各占50%。

该方法的递推公式:

在接下来的章节中，我们将会看到，这两种中值方法，本质上都是二阶龙格库塔方法的特例。

下面是第一种中值法的程序：

clear;
 
center(0.1);
 
function center(h)
%中点法解微分方程:y'=y     真实解为y=e^t
%形参：自变量的步长h
hold on;
 
idx = 1: 1: 10/h;
 
x(1)=0;
y(1)=1;
 
for i = idx    
    x(i+1) = x(i) + h;
    dy = y(i);    
    y_center = y(i) + h/2 * dy;%预测y(xi+h/2)函数值    
    dy = y_center;%中点的导数代替左端点的导数进行迭代
    y(i+1) = y(i) + h * dy;
end
 
plot(x, exp(x), 'r');
plot(x, y, 'b');
legend('真实解','递推解');
end

可以发现，几乎重合，这个方法跟前面的欧拉法相比，步长相等的情况下，递推精度大幅提高，红蓝曲线几乎重合！

三、龙格库塔法（RungeKutta）求解

上面的方法2，看起来几乎重合，实际上放大以后还是能看到误差的，我们还有精度更高的方法。

https://wenku.baidu.com/view/d69e8f1f77c66137ee06eff9aef8941ea76e4b2c.html

推导过程会用到一元函数泰勒公式、复合函数求导公式。

泰勒公式：

复合函数求导公式（全导数公式）：

我们如果用泰勒二阶截断式来做递推，效果肯定比欧拉法要好，因为欧拉法的本质就是用的泰勒一阶截断式。

还是用前面的例题：已知y'=f(x,y)，y(0)=1，求y(x)。

用泰勒二阶递推公式能提高精度是显而易见的，但是由泰勒公式可见，其缺点是，通过提高阶次来运算时，需要用到二阶导、三阶导.....，有些高阶导并不是那么好求的，龙格库塔递推法就是对这一问题进行优化，用低阶导的泰勒展开来避免掉求高阶导。