最近有用python解微分方程的需求，然后找了网上很多的资料和帖子，然后结合个人的想法做了一些研究。

本篇博客主要讲解了二阶常系数微分方程的解析解和数值解的解法以及一阶常系数微分方程组的通解解法，并对每一个步骤做了详细的说明（网上很多帖子讲的并不是很详细），并做了一些探究性的尝试

我们这里用到的库是sympy；

用sympy.dslove常系数齐次微分方程代码思路如下：

例题：解微分方程 f ( x ) ′ ′ + ω 2 f ( x ) = 0 f(x)''+\omega^2f(x)=0 f(x)′′+ω2f(x)=0

现在对照上面的三步很简单的思路来看一个例子：

输出的结果：

这个结果可以化成对应的三角函数形式。

我们拿同济高数书上第346页上的一道例题来检验一下：求方程 y ( 4 ) − 2 y ′ ′ ′ + 5 y ′ ′ = 0 y^{(4)}-2y'''+5y''=0 y(4)−2y′′′+5y′′=0的通解。

代码：

输出：

对照书上的答案，是正确的。

我们还是拿同济高数书来检验一道例题：求方程通解 y ′ ′ − 5 y ′ + 6 y = x e 2 x y''-5y'+6y=xe^{2x} y′′−5y′+6y=xe2x

代码：

输出结果：

谜一般的输出，只有 C 1 ∗ e x p ( 2 ∗ x ) + C 2 ∗ e x p ( 3 ∗ x ) C1*exp(2*x) + C2*exp(3*x) C1∗exp(2∗x)+C2∗exp(3∗x)这一部分是跟答案相同的，其他长长的数字就很莫名其妙了。

之所以会有前面那堆奇怪的数字可能跟算法有关系，这里就不再深究了，只能说用dslove解非齐次方程并不是很好使。

这里我们用到求数值解的库叫做scipy；

scipy.odeint解数值解的思路是这样：

下面还是看一个例子（scipy官网文档中的例子）：解微分方程 y ′ ′ + b ∗ y ′ + c ∗ s i n y = 0 y''+b*y'+c*siny = 0 y′′+b∗y′+c∗siny=0，其中b，c是参数，对t求导；

按照上面的思路，首先化成两个一阶方程： y ′ ( t ) = z ( t ) (1) y'(t)=z(t) \tag 1 y′(t)=z(t)(1)

z ′ ( t ) = − b ∗ z ( t ) − c ∗ s i n ( y ( t ) ) (2) z'(t)=-b*z(t)-c*sin(y(t)) \tag 2 z′(t)=−b∗z(t)−c∗sin(y(t))(2)

上面的第二步中说的x和y是相对方程而言，拿方程（2）来讲就是将z‘(t)当作y’，z(t)当作y，t当作x，其他的当作参数，化为 y ′ = f ( x , y ) y'=f(x,y) y′=f(x,y)的形式，以便这个odeint解方程使用；

接下来看代码和注释，注释只解释基本方法：

接下来我参考了一下说明文档，对代码中一些关键部分做一个详细的说明。

得到的方程数值解可视化结果如下：

对照官网的答案可知是正确的。

我们还是用dsolve来，整个过程跟解一个方程是差不多的，所以直接上例子：

我们还是用高数书上的例题来验证(357页)：

求微分方程组： { x ′ ( t ) = 3 ∗ x ( t ) − 2 ∗ y ( t ) y ′ ( t ) = 2 ∗ z ( t ) − y ( t ) \begin{cases} x'(t)=3*x(t)- 2*y(t) \\ y'(t)=2*z(t)-y(t) \end{cases} {x′(t)=3∗x(t)−2∗y(t)y′(t)=2∗z(t)−y(t)​ 代码：

只是多加了一个函数关系和一个方程，然后解方程的时候需要把两个方程和两个因变量写成一个数组传进对应的位置，依葫芦画瓢即可，最后得到的结果是这样的：

跟书上答案有略微差别但是仔细看一下就会发现是等价的。

直接看高数书358页例题2：

发现这段代码会报错：

它说没有实现hhh，翻看了一下ode.py文件，查看他能做到的求解微分方程的种类：

确实没有二阶常系数微分方程组这个类别，我的sympy是1.6.2版本的，不知道新版会不会增加这个功能。

总结：python解微分方程比较方便，可以满足绝大部分基础的需求（简单的数学建模方程基本都ok），如果需要更加高级的求解的话，建议转matlab。。。