分享我写的小程序-初值问题常微分方程（组）数值解算

从零开始

最近用python写了个初值问题微分方程数值解算程序，可以算微分方程组（测试过2x2 system,更多维的暂时还没有测试），跟scipy内置odeint函数做了结果对比，差别不大。
特点如下。

• 总计600行(包括空行和注释），包含一个非常简单RK-4th 2x2微分方程组（初值问题）解算函数，一个高精度Dormand-Prince算法函数，一个单变量微分方程实例，一个2x2微分方程组实例，一个带动画的实例。
• 自己搜索文献,wiki,和stackOverflow上的问题，实现了 Dormand-Prince算法，包括定步长和自适应步长方法，可以用于单个微分方程，也可以用于微分方程组。
• 带动画的实例，是一个2x2微分方程组实例，(形如： dx/dt = f(t,x,y), dy/dt = f(t,x,y) )运行结果会显示三张图，一张为Vector_field(向量场，matplotlib中的quiver函数实现），第二张可看作为phase_portrait(相图，用matplotlib中的streamplot函数实现），最后一张图动画显示结果x(t),y(t)随时间变化的过程。
  

如果有社友乐意试试这段小程序，注意以下两点：      

• 程序是用python写的，运行此程序需要安装python3,以及 numpy,scipy(参考对比用），matplotlib 等科学计算库（win10系统社友建议直接安装Anaconda）；
• 可以修改三个实例中的微分方程的定义和更改初值条件，个人觉得还是有点意思的，对于学习常微分方程还是有些帮助。
  

-------------------------------分割线---------------------------------------------------------------------------------------------
Dormand-Prince算法介绍：Dormand-Prince 是龙哥库塔系的中的一个比较“年轻”的成员（80年代设计出来的），目前Matlab的ode45的默认算法和simulink的默认解算器。网上已有Fortran语言版本的实现，（scipy调用的是Fortran的ODEPACK库），wiki介绍页面如下：
https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method
-------------------------------分割线--------------------------------------------------------------------------------------------
推荐一个微分方程课程系列。
https://courses.edx.org/courses/ ... 6.1x+2T2017/course/
https://courses.edx.org/courses/ ... 6.2x+2T2017/course/
https://courses.edx.org/courses/ ... 6.3x+3T2017/course/
以及另一个社友（tommy519）提供的非常好的学习资料：
http://jixietop.cn/forum.php?mod=viewthread&tid=8578

-------------------------------分割线--------------------------------------------------------------------------------------------
如果不想麻烦，最后推荐一个国外大学教授的写java程序（学习上述课程时了解到的），学习微分方程的利器，配合上述课程食用，风味最佳。
http://math.rice.edu/~dfield/dfpp.html
下载其中的 “ pplane.jar ”，电脑中装了java 平台就可以打开。
我也是从这个软件中知道的Dormand-Prince算法，（但这个软件你看不到具体实现形式）
实际上想模仿这个软件，但实力不济，写了个半成品，但还是有所收获。

疯子在雨中咆哮

之前只接触过标准容格库塔，算一阶微分方程都很普遍，二阶的就得把方程拆出俩一阶，就变方程组了。而大部分动力学方程都是二阶常微分方程组，有限差分是常规解法。而目前很多软件都拿荣格库塔做迭代，速度更快。文中提到了一些新的分支，希望楼主可以讲一讲对荣格库塔解二阶常微分方程组的看法。

第三个火枪手

一直用excel的默默点赞

从零开始

疯子在雨中咆哮 发表于 2018-6-3 01:33
之前只接触过标准容格库塔，算一阶微分方程都很普遍，二阶的就得把方程拆出俩一阶，就变方程组了。而大部分 ...

雨中咆哮大侠好，说不上看法，可以说说我学到的东西。
大侠已经说出了关键所在-就是降阶，一维二阶微分方程可转化为二维一阶方程组，会拆分基本就会算了，接下来的步骤基本跟单变量一阶微分方程一样了，把每个单步里每个节点的依次算出来，只不过积分对象不再是单个变量的微分，而是两个变量的微分组成的向量。
请看下图。

疯子在雨中咆哮

从零开始 发表于 2018-6-3 11:10
雨中咆哮大侠好，说不上看法，可以说说我学到的东西。
大侠已经说出了关键所在-就是降阶，一维二阶微分方 ...

多谢大侠回复，有限元做瞬态，拉出来就是这个式子，但a,b,c不再是一个数，分别是质量矩阵，阻尼矩阵，刚度矩阵。挼个库泰就是因为要拆，所以解大规模的有点鸡肋，以前用的少。不过应该是克服了一些问题，看到二十年有比较先驱的有限元软件拿 RK 开始解瞬态，主流的有限元软件解瞬态普遍应该从有限差分进化到挼个库塔，毕竟只要能算，就很稳健，收敛快，还可以用相对有限差分更大的时间步长。不过，方法阶数越高反而还容易抽，比较娇气，看到一些是在拿 2th-RK 在算。RK算常微分比有限差分格调要高，值得琢磨。

回火马氏体

6，顶！！！！！

dxx347

牛X,羡慕中。。。