利用MATLAB Coder自动生成RK45数值积分C程序和MEX动态库

最新推荐文章于 2022-11-12 17:55:02 发布

倪偲001

最新推荐文章于 2022-11-12 17:55:02 发布

阅读量1k

点赞数

分类专栏： coding 文章标签： matlab 算法

本文链接：https://blog.csdn.net/NICAI001/article/details/118033655

版权

coding 专栏收录该内容

13 篇文章 3 订阅

订阅专栏

利用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)，即求解如下问题的 $\boldsymbol{x}(t)$ ：
$\dot{\boldsymbol{x}}=\boldsymbol{f}[\boldsymbol{x}(t), t,p]\\ t_0给定,已知\boldsymbol{x}(t_0),t\in[t_0,t_f]$

其中等式右边为 $p$ 代表一个或多个参数。本篇博客的目的是构建快速计算ode结果的c程序并在MATLAB中调用。

数值积分

龙格库塔积分的单步迭代格式为：
$\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}$
其中 $h$ 代表迭代步长。RK-4的是使用最多的，其具体的单步迭代格式为：
$\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}$
多步计算法，一个时间区间内 $t_i,t_{i+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方程是一个二体引力场轨道动力学问题:
$\left\{\begin{array}{l} \dot{\boldsymbol{r}}=\boldsymbol{v} \\ \dot{\boldsymbol{v}}=-\frac{\mu}{r^{3}} \boldsymbol{r} \end{array}\right.$
参数设定见./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$ 倍，精度也略高、

以下测试变步长.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