Runge-Kutta解微分方程

经典单步显式四阶Runge-Kutta（RK4）

1、计算过程：

$$\{\begin{array}{rl}{x}_{n+1}& =x_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(t_n,x_n)\\ k_2& =f(t_n+\frac{h}{2},x_n+\frac{h}{2}{k}_1)\\ k_3& =f(t_n+\frac{h}{2},x_n+\frac{h}{2}{k}_2)\\ k_4& =f(t_n+h,x_n+h{k}_3)\end{array}$$

2、matlab代码

function [x,Y] = rk(f,dt,x,Y)

if nargin < 4 
    Y = {};
end

if ~isempty(Y)
    if  strcmp(class(f),'function_handle')
        k1 = x;
        dk1 = feval(f,k1,Y{1}) * dt;
        k2 = k1 + 0.5 * dk1;
        dk2 = feval(f,k2,Y{2}) * dt;
        k3 = k2 + 0.5 * dk2;
        dk3 = feval(f,k3,Y{3}) * dt;
        k4 = k3 + dk3;
        dk4 = feval(f,k4,Y{4}) * dt;
    else
        k1 = x;
        dk1 = f(k1,Y{1}) * dt;
        k2 = k1 + 0.5 * dk1;
        dk2 = f(k2,Y{2}) * dt;
        k3 = k2 + 0.5 * dk2;
        dk3 = f(k3,Y{3}) * dt;
        k4 = k3 + dk3;
        dk4 = f(k4,Y{4}) * dt;
    end
else
    if  strcmp(class(f),'function_handle')
        k1 = x;
        dk1 = feval(f,k1) * dt;
        k2 = k1 + 0.5 * dk1;
        dk2 = feval(f,k2) * dt;
        k3 = k2 + 0.5 * dk2;
        dk3 = feval(f,k3) * dt;
        k4 = k3 + dk3;
        dk4 = feval(f,k4) * dt;
    else
        k1 = x;
        dk1 = f(k1) * dt;
        k2 = k1 + 0.5 * dk1;
        dk2 = f(k2) * dt;
        k3 = k2 + 0.5 * dk2;
        dk3 = f(k3) * dt;
        k4 = k3 + dk3;
        dk4 = f(k4) * dt;
    end
end

Y = {k1 k2 k3 k4};
x = x + 1/6 * (dk1 + 2*dk2 + 2*dk3 + dk4); 

3、求解算例

单变量：

微分方程：$\frac{df}{dx}=sin(x+3)+x^2$

x = 0.2; dt = 0.1;
f = @(x)(sin(x+3)+x.^2);  % f = inline('sin(x+3)+x.^2','x')
X = [];
for i = 1:100
    x = rk(f,dt,x);
    X = [X x];
end

多变量：

微分方程：
$$\{\begin{array}{rl}\frac{df}{dx}& =sin(x+y)+x^2\\ \frac{df}{dy}& =e^{-x+y}+sin(x+y)\end{array}$$

 XY = [];
 xy = [1;2];
 fxy = @(xy)([sin(xy(1)+xy(2))+xy(1).^2;exp(-xy(1)+xy(2))+sin(xy(1)+xy(2))]);
 for i=1:50
    xy = rk(fxy,0.1,xy);
    XY = [XY xy];
 end

参考

https://zhuanlan.zhihu.com/p/34261490