【Matlab】使用龙格库塔方法求积分

在实际工作中，有一些需要求积分的场合，突然想到可否使用龙格库塔的方式求积分，然后就查找了相关的资料并使用了一个简单的函数验证了一下。

基本原理：

在各种龙格－库塔法当中有一个方法十分常用，以至于经常被称为“RK4”或者就是“龙格－库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。 [1]

令初值问题表述如下：

则，对于该问题的RK4由如下方程给出：

其中

这样，下一个值（yn+1）由现在的值（yn）加上时间间隔（h）和一个估算的斜率的乘积所决定。该斜率是以下斜率的加权平均：

• k1是时间段开始时的斜率；
• k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；
• k3也是中点的斜率，但是这次采用斜率k2决定y值；
• k4是时间段终点的斜率，其y值用k3决定。

当四个斜率取平均时，中点的斜率有更大的权值：

RK4法是四阶方法，也就是说每步的误差是h阶，而总积累误差为h阶。

注意上述公式对于标量或者向量函数（y可以是向量）都适用。

若已知函数f：

可按照上述式子，根据本次函数值以及函数关系 f ，预测下次的函数值 y(n + 1)，进而求得积分

若未知函数 f 求过往积分：

若未知函数 f 求过往值的积分时，在上述的龙格库塔原理中需要由表达式预测出下一步的各个斜率，因此实际，可以将公式中要计算的的 y(n + 1) 改为 y(n) 即：将要预测的下一步函数值修改为本次函数值； 将 y(n) 改为 y(n - 1) 即：使用上一步的函数值来得出预测的本次函数值，使用本次预测函数值来进行积分运算。

全部代码如下：

%*************************************************************************%
                %龙格库塔积分法
%**********************************************************************%

clc;
clear;

%定义一个数组存储时间
Ts = 0.001;    % 离散周期0.01秒
FinalT = 1;  %终止时间
t = Ts : Ts : FinalT;   % 1k
y = zeros(1,(FinalT / Ts));

%计算 y = t^2 的图像
for i = 1 : 1 : (FinalT / Ts)
    y(i) = t(i) * t(i);
end

sum = 0;  %声明一个变量存储积分
yRK = zeros(1,(FinalT / Ts));  %使用一个数组存储龙格库塔计算出来的图像数据
h = Ts;   %h代表时间间隔

%按照四阶龙格库塔法求解
for i = 1 : 1 : (FinalT / Ts)   
    if i > 1
        k1 = ( (y(i) - y(i - 1)) / h);  %计算初始斜率k1
        k2 = ( y(i-1) + k1*h/2 - y(i-1) ) / (h/2);   %计算斜率为k1时对应的中点斜率
        k3 = ( y(i-1) + k2*h/2 - y(i-1) ) / (h/2);   %计算斜率为k2时对应的中点斜率
        k4 = ( y(i-1) + k3*h - y(i-1) ) / h;
        yRK(i) = yRK(i-1) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
        sum = sum + yRK(i) * h;   %计算积分
    else
        k1 = 0;
        k2 = 0;
        k3 = 0;
        k4 = 0;
        yRK(1) = y(1);
%        sum = 0;  %在初始时，积分为0
    end
    
%    sum = y()
end

subplot(2,1,1),plot(t,y);   %函数图像
grid on;
title("y = t^2 图像");

subplot(2,1,2),plot(t,yRK);   %函数图像
grid on;
title("龙格库塔 图像");

积分结果与图像如下：