最近一段时间，博主在学习利用 gpu 加速 matlab 计算，有了一些心得，拿出来和大家分享。

先声明一下版本：

matlab：R2016b

gpu：GeForce GTX 745

博主的研究课题与库仑三体系统相关，需要计算电子的经典运动轨迹，不可避免地用到 ode 求解器。遗憾的是，能够在 gpu 上运行的 matlab 重载函数并不包括 ode 函数，于是博主选择了在 matlab 中定义 gpu 内核（cuda编程的优点更多，前提是你得会编）。注意：matlab 的帮助文档明确指出，只允许用标量算术方法定义在 gpu 上执行的内核，这意味着所有的矢量、矩阵和数组都要拆成分量形式。

能够在内核中使用的内置 matlab 函数和操作，总共139个。

以最常用的 ode45 为例，它是显式龙格-库塔 (4,5) 对的一种变步长积分算法，每步有 7 个阶段。由于使用了具有一致性的 FSAL 方案，仅需估算 6 个函数值就能求得 4/5 阶精确解。

参考 matlab 的源代码，先完成一步积分，再加入步长控制和误差估计。

矩阵 A、B 和 E 构成了一个 Butcher 表格，矩阵 B 中的每一列元素的加和等于矩阵 A 中的每一列元素。neq 为常微分方程的个数，每一个常微分方程都要估算 7 个函数值，第 7 个函数值作为下一步积分中的第 1 个函数值，实际上每一步只用估算 6 个函数值。

odeFcn 是要求解的常微分方程；除去两个端点的函数值，还要在时间步骤 h 中估算 5 个函数值，因此 RK5 用来传播解。

absh * f * E 为 4/5 阶精确解的差值；y、ynew 的绝对值过小，会导致 err 过大，加入阈值项 (绝对误差 atol 除以相对误差 rtol 的商) 予以修正。

计算得到的 err 大于相对误差 rtol，缩小积分步长；反之，如果在一步积分中第一次计算得到的 err 小于相对误差 rtol，则增大积分步长。与定步长相比，变步长的计算量更小，精度也更高。

参考：微信公众号SPACEofPHD ( http://url.cn/58N4aRl  )