利用MATLAB Coder生成变步长ode45数值积分C程序
Use MATLAB Coder to generate C code and MEX file for MATLAB, an example of RK45 and RKF87 for numerical integration
MATLAB虽然用着方便,但是其诟病在于太大的计算量以及很慢的计算速度,而C语言更偏向于底层,计算速度比MATLAB快很多。本博客所解决的数值积分(numerical integration)旨在求解常微分方程组的初值问题(initial value problem),即求解如下问题的
x
(
t
)
\boldsymbol{x}(t)
x(t):
x
˙
=
f
[
x
(
t
)
,
t
,
p
]
t
0
给
定
,
已
知
x
(
t
0
)
,
t
∈
[
t
0
,
t
f
]
\dot{\boldsymbol{x}}=\boldsymbol{f}[\boldsymbol{x}(t), t,p]\\ t_0给定,已知\boldsymbol{x}(t_0),t\in[t_0,t_f]
x˙=f[x(t),t,p]t0给定,已知x(t0),t∈[t0,tf]
其中等式右边为 p p p代表一个或多个参数。本篇博客的目的是构建快速计算ode结果的c程序并在MATLAB中调用。
数值积分
龙格库塔积分的单步迭代格式为:
x
i
+
1
=
x
i
+
h
i
∑
j
=
1
k
β
j
f
i
j
f
i
j
=
f
[
(
x
i
+
h
i
∑
l
=
1
k
α
j
l
f
i
l
)
,
(
t
i
+
h
i
ρ
j
)
]
\begin{aligned} x_{i+1}=&x_{i}+h_{i} \sum_{j=1}^{k} \beta_{j} f_{i j} \\ f_{i j}=&f\left[\left(x_{i}+h_{i} \sum_{l=1}^{k} \alpha_{j l} f_{i l}\right),\left(t_{i}+h_{i} \rho_{j}\right)\right] \end{aligned}
xi+1=fij=xi+hij=1∑kβjfijf[(xi+hil=1∑kαjlfil),(ti+hiρj)]
其中
h
h
h代表迭代步长。RK-4的是使用最多的,其具体的单步迭代格式为:
x
i
+
1
=
x
i
+
h
6
(
K
1
+
2
K
2
+
2
K
3
+
K
4
)
K
1
=
f
(
t
i
,
x
i
)
K
2
=
f
(
t
i
+
h
2
,
x
i
+
h
2
K
1
)
K
3
=
f
(
t
i
+
h
2
,
x
i
+
h
2
K
2
)
K
4
=
f
(
t
i
+
h
,
x
i
+
h
K
3
)
\begin{aligned} x_{i+1} &=x_{i}+\frac{h}{6}\left(K_{1}+2 K_{2}+2 K_{3}+K_{4}\right) \\ K_{1} &=f\left(t_{i}, x_{i}\right) \\ K_{2} &=f\left(t_{i}+\frac{h}{2}, x_{i}+\frac{h}{2} K_{1}\right) \\ K_{3} &=f\left(t_{i}+\frac{h}{2}, x_{i}+\frac{h}{2} K_{2}\right) \\ K_{4} &=f\left(t_{i}+h, x_{i}+h K_{3}\right) \end{aligned}
xi+1K1K2K3K4=xi+6h(K1+2K2+2K3+K4)=f(ti,xi)=f(ti+2h,xi+2hK1)=f(ti+2h,xi+2hK2)=f(ti+h,xi+hK3)
多步计算法,一个时间区间内
[
t
i
,
t
i
+
1
]
[t_i,t_{i+1}]
[ti,ti+1]内会进行多次RK4求解。RK45不同于RK4,还加了一个估计精度的步骤,以此确定下一步变步长的长度,然后进行多步求解。这个算法也叫Runge-Kutta-Fehlberg 4(5) ,也就是MATLAB自带函数[T,Y,TE,YE,IE] = ode45(odefun,tspan,y0,options)
所采用的算法,它的核心部分是以下代码:
function [x, errorEstimate] = RK45( x, t, dt ,var)
% REFERENCE1: https://github.com/recalon/rkf_calon/rkf_calon.cpp
k1 = yprime(t,x ,var);
k2 = yprime(t+dt/4, x+dt*k1/4 ,var);
k3 = yprime(t+3*dt/8, x+3*dt*k1/32+9*dt*k2/32 ,var);
k4 = yprime(t+12*dt/13, x+1932*dt*k1/2197-7200*dt*k2/2197+7296*dt*k3/2197 ,var);
k5 = yprime(t+dt, x+439*dt*k1/216-8*dt*k2+3680*dt*k3/513-845*dt*k4/4104 ,var);
k6 = yprime(t+dt/2, x-8*dt*k1/27+2*dt*k2-3544*dt*k3/2565+1859*dt*k4/4104-11*dt*k5/40 ,var);
x = x+ dt*(k1*25/216 + k3*1408/2565 + k4*2197/4104 - k5/5);
by1 = 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5;
by2 = 16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55;
errorEstimate = by2-by1;
% errorEstimate = k1/360 - k3*128/4275 - k4*2197/75240 + k5/50 + k6*2/55;
end
对应C函数可以查看GitHub: kf_calon/rkf_calon.cpp .
操作步骤
以下是操作流程图,具体步骤不再赘述。源代码与生成代码均已上传GitHub。
测试MATLAB算法
第一部分编写MATLAB代码完成之后,编写测试文件测试。这里使用的ode方程是一个二体引力场轨道动力学问题:
{
r
˙
=
v
v
˙
=
−
μ
r
3
r
\left\{\begin{array}{l} \dot{\boldsymbol{r}}=\boldsymbol{v} \\ \dot{\boldsymbol{v}}=-\frac{\mu}{r^{3}} \boldsymbol{r} \end{array}\right.
{r˙=vv˙=−r3μr
参数设定见./testM.m
,主要调用以下函数[ts,yout,tf,yend] = wrapRK45(tspan,y_0,options,var);
测试算法是否无误,以便后续转为C代码。
测试固定步长的计算结果为:
Fixed step test
时间已过 1.004216 秒。
Runge-Kutta-45: results step: 61
Accuracy Digits
2.00e-12 11.70
时间已过 0.448077 秒。
ode45 MATLAB; results step: 61
Accuracy Digits
8.57e-11 10.07
可见自己编写的RK45计算效率并不占优势。
测试变步长的计算结果为:
Variable step test
时间已过 0.364811 秒。
Runge-Kutta-45: results step: 7090
Accuracy Digits
2.12e-10 9.67
时间已过 0.636274 秒。
ode45 MATLAB; results step: 28025
Accuracy Digits
8.57e-11 10.07
可见计算时间略短于ode45,由于步长更大,实际上计算效率下降了。
MATLAB Coder转化多个.m文件
定义输入格式
这一步需要做的准备是定义C文件的调用格式,明确wrapRK45(tspan,y_0,options,var);
输入参数的类型. myIO.m文件的主要部分如下:
% time span series or t0 tf
tspan=[0:60]';
tspan =[0, 60]';
% 由于时间向量不一定完全给定,因此它的输入大小不固定
assert(all(size(tspan)<=[1e6 1]));
% initial state
y_0=[0.9;0;0;0;1.1055;0];
% options for RKF45, also could use 'odeset' in matlab code, but not in C
options.RelTol=1e-13;
options.MaxStep=0.01;
% auxilary variables in ode system
var.mu=GM;
var.Cd=2.2;
% 明确所有输入类型
[ts,yout,tf,yend] = wrapRK45(tspan,y_0,options,var);
这个文件决定了后面产生的.c文件和.mex文件的接口方式,因此非常重要。可以参考MATLABDefine Input Properties Programmatically in the MATLAB File相关章节。
转换
我所采用的是MATLAB R2015b和Visual Studio 2013. 构建mex文件的过程应该要求存在对应版本的编译器compiler,否则有可能出错。哪个版本MATLAB对应哪个编译器,可以查看官网:System Requirements and Supported Compilers
用这个,打开然后按照以下步骤
-
选中所有有依赖性的.m文件
-
指定输入类型。由于输入较多,这个过程一般是用一个脚本,包含了对所有要转化的.m文件的输入格式。
如果多个文件互相依赖,只需要给出主程序的输入参数类型即可。
-
输入类型可以手动更改。一些大小不固定的输入参数,应该在脚本中设置。必须与myIO.m文件一致,否则会报错。另外,结构体的输入字段的个数和名称都应该一致
-
点check for issues, 将会自动执行: 产生.c .h代码->构建mex动态库->测试.mex库。
-
点generate,如果报错,直接更改源文件即可。当然如果是比较大的错误,建议重头上面步骤
-
产生了所有的.c .h头文件。
-
完成。将得到以下文件列表,其中RK45_mex.mex64即为包含了所选3个.m文件的动态链接库。
-
最下边的RK45.prj是可以重复使用和调试的。注意更改.m文件后重新编译.mex时别忘了重新运行第4步check for issues重新生成mex,运行第6步generate和verify code以产生新的.c源文件.
Note: 由于变步长算法正确的实现方式有应该采用动态分配内存,也就是C语言中的malloc函数;而MATLAB转C语言不会自动添加这一句,如果后续使用,可以手动修改第6步的源文件。或者更改步骤5中的
测试
要调用.mex64里的函数而不是原来的.m源文件,应该按照mexFileName('functionName',allParametersInPut)
这样的调用格式,即[ts,yout,tf,yend] = RK45_mex('wrapRK45',tspan,y_0,options,var);
。
建立testMEX.m脚本进行测试,固定步长计算结果如下:
Fixed step test
>>时间已过 0.039140 秒。
Runge-Kutta-45 .MEX64 version; results calculation step: 61
Accuracy Digits
1.74e-12 11.76
>>时间已过 0.471013 秒。
ode45 MATLAB; results calculation step: 61
Accuracy Digits
8.57e-11 10.07
可见运算速度相较其.m文件版本提高了 1.00 / 0.04 = 25 1.00/0.04=25 1.00/0.04=25倍,精度也略高、
以下测试变步长.mex版本:
Variable step test
时间已过 0.152232 秒。
Runge-Kutta-45 .MEX64 version; results calculation step: 7090
Accuracy Digits
2.13e-10 9.67
时间已过 0.618282 秒。
ode45 MATLAB; results calculation step: 28025
Accuracy Digits
8.57e-11 10.07
可见由于我没有更改动态分配内存这一部分的源代码,其计算效率只和ode45大致一样。
如果能优化好C语言代码,再重新用mex编译一个.mex64动态库,计算效率应该会好一点。当然,MATLAB Coder转C语言只能说是给C程序编写打下了一个基础,其中的关键代码是可能复用的。这只能在.m文件较为简单的情况下起效,否则其复杂的调用框架和数据结构并不会比.m文件快多少。