matlab实现RK45（Runge-Kutta45、ode45）求解器算法

lidenghui

RK45求解器，又称为Dormand-Prince求解器。这是比较精确的求解器，可以快速地求解微分方程，但是，需要消耗一些内存。在matlab simulink中默认条件下，系统自动选择RK45求解器。用户可以根据实际问题，选择合适的求解器。

Dopri54是Dormand / Prince龙格-库塔的一种方法，Dopri54由龙格-库塔简单法构成，阶为5和4。因此，五阶龙格-库塔法是利用一步向前+四阶龙格-库塔法估计误差。

本文分享一个简单例子来从m代码实现RK45求解器，matlab也可以用自带的ode45函数来求解微分方程:Matlab通过ode系列函数求解微分方程

假定y'=y，y(0) = 1，很明显结果的理论解为y(t)=e^t，

matlab代码

clcclose allcleary0 = 1;[t,y] = dopri54c('fun', 0, 1, y0, 0.0001);
figureplot(t,y)function u = fun(x, y)u=y;function[t_1, y_1]=dopri54c(funcion,t0,t1,y0,tol)%Dormand and Prince 54 code.%INPUT%funcion - 函数句柄%t0 - 开始时间.%t1 - 结束时间.%y0 - 初始值.%tol - 局部误差%OUTPUT%y - 输出t=t0;y=y0;h=tol^(1/5)/4;step=0;nrej=0;fcall=1;a4=[44/45 -56/15 32/9]';a5=[19372/6561 -25360/2187 64448/6561 -212/729]';a6=[9017/3168 -355/33 46732/5247 49/176 -5103/18656]';a7=[35/384 0 500/1113 125/192 -2187/6784 11/84]';e=[71/57600 -1/40 -71/16695 71/1920 -17253/339200 22/525]';k1=feval(funcion,t,y);t_1 = t;y_1 = y;
whilet < t1ift+h > t1h=t1-t; endk2=feval(funcion,t+h/5,y+h*k1/5);k3=feval(funcion,t+3*h/10,y+h*(3*k1+9*k2)/40);k4=feval(funcion,t+4*h/5,y+h*(a4(1)*k1+a4(2)*k2+a4(3)*k3));k5=feval(funcion,t+8*h/9,y+h*(a5(1)*k1+a5(2)*k2+a5(3)*k3+a5(4)*k4));k6=feval(funcion,t+h,y+h*(a6(1)*k1+a6(2)*k2+a6(3)*k3+a6(4)*k4+a6(5)*k5));yt=y+h*(a7(1)*k1+a7(3)*k3+a7(4)*k4+a7(5)*k5+a7(6)*k6);k2=feval(funcion,t+h,yt);est=norm(h*(e(1)*k1+e(2)*k2+e(3)*k3+e(4)*k4+e(5)*k5+e(6)*k6),inf);fcall=fcall+6;
ifest < tolt=t+h;k1=k2;step=step+1;y=yt;t_1 = [t_1; t];y_1 = [y_1; y];elsenrej=nrej+1;endh=.9*min((tol/(est+eps))^(1/5),10)*h;
end





结果符合预期

suzuki

关于RK45求解器（也被称为Dormand-Prince求解器）在MATLAB中的实现，ode45函数已经内置了此算法，用户无需自行编写。但为了满足您的需求，我可以提供一个简单的示例代码，展示如何使用MATLAB的ode45函数求解微分方程。

1. matlab<br>% 定义微分方程 dy/dt = f(t, y) 的函数<br>function dy = myfun(t, y)<br> % 在这里填写你的微分方程式<br> dy = ...; % 以矩阵形式返回dy，对应于y的每个分量<br>end<br><br>% 初始条件<br>y0 = [initialvalue1; initialvalue2; ...]; % 根据问题设定初始值<br>tspan = [starttime, endtime]; % 定义时间区间<br><br>% 使用ode45求解微分方程<br>[t, y] = ode45(@myfun, tspan, y0); % @myfun是之前定义的函数句柄<br><br>% 绘图展示结果<br>plot(t, y(:, 1), 'b-', tspan(1), y0(1), 'ro'); % 绘制第一个分量的结果，可根据需要绘制其他分量<br>xlabel('Time'); ylabel('Solution'); title('Solution of Differential Equation');<br>legend('Solution', 'Initial Point'); grid on;

复制代码

在上面的代码中，myfun函数需要您根据具体问题定义微分方程式。ode45函数会自动选择RK45算法（如果适用的话）来求解微分方程。返回的t和y分别代表时间向量和对应的解向量。在绘图部分可以根据需求绘制解的变化曲线。需要注意的是，对于复杂的系统可能需要更高级的设置和调试过程。