![cover](https://img-blog.csdnimg.cn/f60dec47487b4ae7a1f653f3b19577e5.jpeg#pic_center)
MATLAB中ode45()和Runge-Kutta算法(4阶)的比较
写这篇博客目的是自己在求解微分方程的时候,考虑到 `ode45()` 可能求解速度比较慢,想用一种快速一点的微分方程求解算法,又想到 `ode45()` 用的就是 Runge-Kutta (RK)算法,所以想是不是自己可以自己编写一个 RK 算法用来代替 `ode45()`。因为 MATLAB 的 `ode45()` 算法中可能存在较多的判断条件,这也许是让 `ode45()` 速度较慢的一个原因
引言
写这篇博客目的是自己在求解微分方程的时候,考虑到 ode45()
可能求解速度比较慢,想用一种快速一点的微分方程求解算法,又想到 ode45()
用的就是 Runge-Kutta (RK)算法,所以想是不是自己可以自己编写一个 RK 算法用来代替 ode45()
。因为 MATLAB 的 ode45()
算法中可能存在较多的判断条件,这也许是让 ode45()
速度较慢的一个原因,而自己编写的 RK 算法省略了很多不必要的判断,也许会快一些。
ode45()
ode45()
算法是 MATLAB 中专门用于求解常微分方程(Ordinary differential equations,ODE)的函数,它在求解微分方程时的步长是变步长的,使用也是 RK 算法。
45
代表了 ode45()
使用 4 阶- 5阶 RK 算法。其中,4 阶方法用于提供 ODE 的候选解,5 阶方法用于控制误差,整体截断误差为
Δ
x
5
\Delta x^5
Δx5。
Runge-Kutta 算法
RK 算法是德国数学家 Runge 和 Kutta 在 1900 年前后提出的一种高精度求解微分方程的数值算法。该算法间接使用 Taylor 级数展开构造高精度数值方法,利用函数 f f f 在若干点的值的线性组合代替 f f f 的导数,再利用 Taylor 级数展开方法确定其中系数。
下面详细说明 RK 算法如何求解一阶微分方程。
对于如下问题:
{ y ˙ = f ( t , y ) , y ( t 0 ) = y 0 , \left \{\begin{matrix} \begin{aligned} &\boldsymbol{\dot{y}} = f(t,\boldsymbol{y}),\\ &\boldsymbol{y}(t_0) = \boldsymbol{y}_0, \end{aligned} \end{matrix} \right. {y˙=f(t,y),y(t0)=y0,
式中, t 0 t_0 t0 为初始时间, y 0 \boldsymbol{y}_0 y0 为初始状态, f ( t , y ) f(t,\boldsymbol{y}) f(t,y) 为关于时间 t t t 和状态 t \boldsymbol{t} t 的函数。
那么 4 阶 RK 算法为:
{ k 1 = f ( t ( k ) , y ( k ) ) , k 2 = f ( t ( k ) + h / 2 , y ( k ) + k 1 h / 2 ) , k 3 = f ( t ( k ) + h / 2 , y ( k ) + k 2 h / 2 ) , k 4 = f ( t ( k ) + h , y ( k ) + k 3 h ) , y ( k + 1 ) = y ( k ) + h ( k 1 + 2 k 2 + 2 k 3 + k 4 ) / 6 , y ( 0 ) = y 0 . \left \{\begin{matrix} \begin{aligned} &k_1 = f(t(k) ,\boldsymbol{y}(k)), \\ &k_2 = f(t(k) + h/2,\boldsymbol{y}(k) + k_1 h/2), \\ &k_3 = f(t(k) + h/2,\boldsymbol{y}(k) + k_2 h/2), \\ &k_4 = f(t(k) + h ,\boldsymbol{y}(k) + k_3 h), \\ &\boldsymbol{y}(k+1) = \boldsymbol{y}(k) + h(k_1+2k_2+2k_3+k_4)/6, \\ &\boldsymbol{y}(0) = \boldsymbol{y}_0. \end{aligned} \end{matrix} \right. ⎩ ⎨ ⎧k1=f(t(k),y(k)),k2=f(t(k)+h/2,y(k)+k1h/2),k3=f(t(k)+h/2,y(k)+k2h/2),k4=f(t(k)+h,y(k)+k3h),y(k+1)=y(k)+h(k1+2k2+2k3+k4)/6,y(0)=y0.
上式即为 y ( k ) \boldsymbol{y}(k) y(k) 向 y ( k + 1 ) \boldsymbol{y}(k+1) y(k+1) 递推的形式,根据初始条件可以求出 y ( 1 ) \boldsymbol{y}(1) y(1), y ( 2 ) \boldsymbol{y}(2) y(2), y ( 3 ) \boldsymbol{y}(3) y(3), y ( 4 ) \boldsymbol{y}(4) y(4), ⋯ \cdots ⋯ y ( N ) \boldsymbol{y}(N) y(N),该离散序列即为微分方程的数值解。
RK 算法程序
在这里给出 RK 算法的 MATLAB 程序。
%--- RK 算法 ---%
% ufunc - 微分方程组的函数名
% y0 - 初始值
% h - 步长
% a - 初始时间
% b - 末端时间
function [x,y]=runge_kutta1(ufunc,y0,h,a,b)
% 步数
n = floor((b-a)/h);
% 时间起点
x(1) = a;
% 初始值
y(:,1)=y0;
% 算法开始
for iter=1:n
% 时间变量
x(iter+1) = x(iter)+h;
% 开始迭代
k1 = ufunc(x(iter), y(:,iter));
k2 = ufunc(x(iter) + h/2, y(:,iter) + h*k1 / 2);
k3 = ufunc(x(iter) + h/2, y(:,iter) + h*k2 / 2);
k4 = ufunc(x(iter) + h, y(:,iter) + h*k3 );
% 得到结果
y(:,iter+1) = y(:,iter) + h * (k1 + 2*k2 + 2*k3 + k4) / 6;
end
end
仿真
用一个一阶微分方程数值求解例子来对比 ode45()
和 自编 RK 算法的效果。
对于问题
{ y ˙ 1 = y 2 y 3 , y ˙ 2 = − y 1 + y 3 , y ˙ 3 = − 0.51 y 1 y 2 , y = ( y 10 , y 20 , y 30 ) = ( 1 , 1 , 3 ) . \left \{\begin{matrix} \begin{aligned} \dot{y}_1 &= y_2 y_3, \\ \dot{y}_2 &= -y_1 + y_3, \\ \dot{y}_3 &= -0.51 y_1 y_2, \\ \boldsymbol{y} &= (y_{10},y_{20},y_{30}) = (1,1,3). \end{aligned} \end{matrix} \right. ⎩ ⎨ ⎧y˙1y˙2y˙3y=y2y3,=−y1+y3,=−0.51y1y2,=(y10,y20,y30)=(1,1,3).
在时间区间 [ 0 , t f ] [0,t_f] [0,tf] 上的解,其中 t f t_f tf 未定。
仿真代码
给出 MATLAB 解决上述问题的仿真代码。
clc;clear;close all;
%% 主函数
% 时间变量
t0 = 0;
tf = 1000;
% 初值
y10 = 1;
y20 = 1;
y30 = 3;
y = [y10,y20,y30];
% RK 算法步长
h = 0.25;
% 对该微分方程组用ode45和自编的龙格库塔函数进行比s较,调用如下:
% ode45()函数
tic;
[T,F] = ode45(@fun,[t0 tf],y);
time_record_ode = toc;
toc;
% 自编RK函数
tic;
[T1,F1]=runge_kutta1(@fun,y,h,t0,tf);
time_record_rk = toc;
toc;
%% 画图
figure('color',[1 1 1],'position',[600,100,500*1.5,416*1.5]);
subplot(121);
plot(T,F,'LineWidth',1.5);
title(['ode45','($t_f$=',num2str(tf),')'],'Interpreter','Latex');
set(gca,'FontSize',15,'FontName','Times New Roman','LineWidth',1.5);
subplot(122);
plot(T1,F1,'LineWidth',1.5);
title(['RK4','($t_f$=',num2str(tf),'),($h$=',num2str(h),')'],'Interpreter','Latex');
set(gca,'FontSize',15,'FontName','Times New Roman','LineWidth',1.5);
% 保存数据
str = ['tf_',num2str(tf),'_h_',num2str(h),'.mat'];
str_fig = ['tf_',num2str(tf),'_h_',num2str(h),'.jpg'];
save(str);
saveas(gcf,str_fig)
%% 子函数部分
% 微分方程
function dy = fun(t,y)
dy = zeros(3,1);%初始化列向量
dy(1) = y(2) * y(3);
dy(2) = -y(1) + y(3);
dy(3) = -0.51 * y(1) * y(2);
end
%--- RK 算法 ---%
% ufunc - 微分方程组的函数名
% y0 - 初始值
% h - 步长
% a - 初始时间
% b - 末端时间
function [x,y]=runge_kutta1(ufunc,y0,h,a,b)
% 步数
n = floor((b-a)/h);
% 时间起点
x = zeros(n,1);
x(1) = a;
% 初始值
y(:,1)=y0;
% 算法开始
for iter=1:n
% 时间变量
x(iter+1) = x(iter)+h;
% 开始迭代
k1 = ufunc(x(iter), y(:,iter));
k2 = ufunc(x(iter) + h/2, y(:,iter) + h*k1 / 2);
k3 = ufunc(x(iter) + h/2, y(:,iter) + h*k2 / 2);
k4 = ufunc(x(iter) + h, y(:,iter) + h*k3 );
% 得到结果
y(:,iter+1) = y(:,iter) + h * (k1 + 2*k2 + 2*k3 + k4) / 6;
end
end
算法主要对比的有 2 项指标,计算精度和计算效率。
不同终端时间下的算法对比
分别对比
t
f
=
15
,
1
0
2
,
1
0
3
,
1
0
4
,
1
0
5
t_f=15,10^2,10^3,10^4,10^5
tf=15,102,103,104,105 时,ode45()
和 RK 算法的计算精度和计算效率。
从计算精度来看,计算时间比较小的时候, RK 算法和 ode45()
的区别不是很大,但是计算时间增大之后,到了
t
f
=
1
0
3
,
1
0
4
,
1
0
5
t_f=10^3,10^4,10^5
tf=103,104,105 的时候,RK 算法的精度就逐渐跟不上 ode45()
,此时 RK 算法累积的误差更大一些了。
t f = 15 t_f=15 tf=15 | t f = 1 0 2 t_f=10^2 tf=102 | t f = 1 0 3 t_f=10^3 tf=103 | t f = 1 0 4 t_f=10^4 tf=104 | t f = 1 0 5 t_f=10^5 tf=105 | |
---|---|---|---|---|---|
ode45() | 0.000863 | 0.002703 | 0.020696 | 0.165518 | 1.204862 |
RK4() | 0.000448 | 0.001650 | 0.015817 | 0.136802 | 1.336549 |
从计算效率来看,RK 算法的计算效率要优于 ode45()
,不过计算时间增大到
t
f
=
1
0
5
t_f=10^5
tf=105 的时候,ode45()
的计算效率超过了 RK 算法,猜想是计算时间太长,自编的 RK 算法无法针对这么长计算时间的微分方程做特定优化,所以导致计算时间慢于 ode45()
。
不同步长下的算法对比
分别对比
t
f
=
15
t_f=15
tf=15 和
t
f
=
1
0
3
t_f=10^3
tf=103 ,Rk 算法步长为
h
=
0.01
,
0.10
,
0.20
,
0.25
,
0.35
,
0.50
h=0.01,0.10,0.20,0.25,0.35,0.50
h=0.01,0.10,0.20,0.25,0.35,0.50 时,ode45()
和 RK 算法的计算精度和计算效率。
从计算精度来看,计算时间比较小的时候, RK 算法的步长越小,算的结果越和 ode45()
的结果接近,步长增加到
h
=
0.35
,
0.5
h=0.35,0.5
h=0.35,0.5 的时候,得到的曲线已经呈现不光滑的趋势了。
计算时间增大之后,,RK 算法精度跟不上 ode45()
,这个时候的步长选取就值得考虑了,可以发现,不论是步长取得过小(
h
=
0.01
,
0.10
,
0.20
h=0.01,0.10,0.20
h=0.01,0.10,0.20),还是步长取得过大(
h
=
0.35
,
0.50
h=0.35,0.50
h=0.35,0.50),得到的结果都和 ode45()
的结果不相符合,在
h
=
0.25
h=0.25
h=0.25 时,Rk 算法求解微分方程的效果会好一些,和 ode45()
的结果有差距,但相较其他步长的结果来说,要更好一些。
t f = 15 t_f=15 tf=15 | h = 0.01 h=0.01 h=0.01 | h = 0.10 h=0.10 h=0.10 | h = 0.20 h=0.20 h=0.20 | h = 0.25 h=0.25 h=0.25 | h = 0.35 h=0.35 h=0.35 | h = 0.50 h=0.50 h=0.50 |
---|---|---|---|---|---|---|
ode45() | 0.000795 | 0.000781 | 0.000881 | 0.000863 | 0.000992 | 0.000870 |
RK4() | 0.005705 | 0.000901 | 0.000493 | 0.000448 | 0.000380 | 0.000251 |
t f = 1 0 3 t_f=10^3 tf=103 | h = 0.01 h=0.01 h=0.01 | h = 0.10 h=0.10 h=0.10 | h = 0.20 h=0.20 h=0.20 | h = 0.25 h=0.25 h=0.25 | h = 0.35 h=0.35 h=0.35 | h = 0.50 h=0.50 h=0.50 |
ode45() | 0.021820 | 0.020883 | 0.020322 | 0.020696 | 0.020698 | 0.020758 |
RK4() | 0.343084 | 0.035581 | 0.017536 | 0.015817 | 0.011288 | 0.007863 |
从计算效率来看,RK 算法的步长越小,计算效率越低,因为步长小的话,每次更新迭代的值也小了,需要更多步才能迭代到理想值;步长大的时候,计算效率较高,但是计算精度较低。
结论
因此,可以得出结论,
- 在计算时间短的微分方程求解问题中,使用自编 RK 算法的计算精度与使用
ode45()
计算精度相近,计算效率更高; - 在计算时间长的微分方程求解问题中,建议使用
ode45()
会得到更好的结果; - RK 算法中,计算步长和计算时间区间是需要协调好的重要变量,把这 2 个变量设置好,计算精度和计算效率才会更加优越。
欢迎通过邮件联系我:lordofdapanji@foxmail.com
![Logo](https://devpress.csdnimg.cn/79de2bf0b7994defa4242ef90d5513fa.jpg)
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)