作者：18级潘愽浩

导言：

今天让我们来了解一下如何用scipy中的odeint函数解一阶常微分方程，不过由于scipy的功能实现是在numpy库的基础之上，我还会简单的为大家介绍numpy库。话不多说，让我们进入正题！

我们先了解一下python在科学计算领域中四个非常受欢迎的库，numpy、scipy、matplotlib，pandas。numpy是一个高性能的多维数组（矩阵）的计算库，SciPy是构建在numpy的基础之上的，它提供了许多的操作numpy的数组的函数。SciPy是一款方便、易于使用、专为科学和工程设计的python工具包，它包括了统计、优化、整合以及线性代数模块、傅里叶变换、信号和图像图例，常微分方差的求解等。Matplotlib 是 Python 的绘图库。它可与 NumPy 一起使用,提供了一种有效的 MatLab 开源替代方案。pandas 是基于NumPy 的一种工具，该工具是为了解决数据分析任务而创建的。Pandas 纳入了大量库和一些标准的数据模型，提供了高效地操作大型数据集所需的工具。这四个库可实现matlab的绝大部分功能。

开讲之前，我们首先要明白：对于一阶常微分方程，odeint函数只能解dy/dx = f(x,y)形式的常微分方程，只能解dy/dx = f(x,y)形式的常微分方程，只能解dy/dx = f(x,y)形式的常微分方程！！！

记得自行cmd安装所需库（pip install scipy；pip install numpy；pip install matplotlib）

odeint函数介绍：

odeint()是scipy库中一个数值求解微分方程的函数。它的完整形式为odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)（太过冗长可不用看）。而实际操作中，我们一般只需要给它传三个参数：第一个是微分方程的函数，第二个是微分方程初值，第三个是微分的自变量。我们就以这三个参数为主，来实现常微分方程的数值求解。

一．定义微分方程的函数：

一般方法为直接定一个函数，将因变量y和自变量x传给它作为参数，返回值就是上文所强调的等式右边的f(x,y)的表达式。

比如说我们要求dy/dx = x的话，就定义函数：

1 2	def fun(y,x): return x

要求dy/dx = (x-1)y^2 + (1-2x)*y + x的话，就定义函数：

1 2	def fun(y,x): retuen (x-1)y2 + (1-2x)*y + x

二．给自变量赋值：

学过matlab的应该都知道，matlab作f（x，y）图一般是定义一组x的值，然后依次求出y的值，最后将这许许多多的（x，y）坐标绘出并连接起来。例如，我们要画一个在-5<=x<=5的定义域上y = x^3的图像，就应该是

1
2
3

x = [-5:0.01:5];
y = x.^3;
plot(x,y)

而在python中我们要画常微分方程的图形也是同样的道理，定义一个一行n列的矩阵，矩阵中的值x0,x1,x2…xn均在定义域中，然后将每个自变量对应的因变量y0,y1,y2…yn的值求出，这样就可以将常微分方程的坐标图画出来啦！

在python中，我们要想定义一个矩阵（array），就需要用到numpy库了。定义一个解常微分方程所需的矩阵则有以下几种形式：

1. 创建列表，再把其转化成array形式。

1 2	a = [0,1,2,3,4,5,6,7,8,9] numpy.array(a)

但是这样我们最终的作图只能描下十个点，直线连接之后的坐标图不够圆滑，所以我们需要盘它，把它的步长调小点：

1 2	a = range(0,10,0.01) numpy.array(a)

2步长和0.01步长对比（以dy/dx = x为例）：

当然，这多一部创造列表的方法还是有些麻烦，我们可以用numpy直接创造所需矩阵：

2. Numpy直接创造array：

可以用numpy.arange(0,10,0.01)，也可以用numpy.linspace(0,10,100)（意思是在0到10范围内等距离取100个值并存放到1行n阶矩阵里）

这样，我们就完成了自变量X的赋值工作。

三．利用x的值根据函数关系求y的值：

这一步比较简单，我们只需要使用我们的主角odeint函数。一行代码即可搞定：

y = odeint(fun,0,x)

在这个函数中，第一个参数fun是我们之前定义好的微分方程表达式的右边f(x,y)的表达式。

至于第二个参数，我们以dy/dx = x（0<=x<=10）（以下简称fun1）和dy/dx = x（5<=x<=10）（以下简称fun2）为例子作讲解：

对于dy/dx = x，我们可求得最终表达式为y = (x^2)/2 + c.

而第二个参数0就是解的初值.比如对于fun1，自变量最小为0，令y(x0) = 0,x0 = 0,所以c的值就是0。对于fun2，自变量最小为5，令y(x0) = 0,x0 = 5,所以c的值为-25/2.也就是说，我们所作出的坐标图永远是从y = 0开始的。

第三个参数就是自变量，这个就没有什么好说的了。

四．根据所求得的x0,x1,x2…,y0,y1,y2…作坐标图：

做图我们需要用到matplotlib库，代码如下:

matplotlib.pyplot.plot(x,y) matplotlib.pyplot.grid() matplotlib.pyplot.show()

这样，我们就可以得到在0<=x<=10的范围内，常数c = 0，dy/dx = x的坐标图：

总代码如下：

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#定义函数
def diff(y,x):
    return x
#定义自变量x的取值范围
x = np.arange(0,10,0.01)
#根据x求因变量的值
y = odeint(diff,0,x)
#作图
plt.plot(x,y)
plt.grid()
plt.show()

我们再以常微分课本74页最后一题为例

代码为：

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def diff(y,x):
    return (x-1)*y**2 + (1-2*x)*y + x
x = np.arange(0,10,0.01)/
y = odeint(diff,0,x)

plt.plot(x,y)
plt.grid()
plt.show()

可以看到，仅定义的函数返回值不一样，其他的全部一样。欢迎各位大佬验证答案的正确性。

当然，专业解微分方程还是要用maple，不过多了解一些其他的也肯定有其好处。

好了，关于一阶常微分方程的数值求解就到这里啦，因为我也是刚开始看，肯定会有写的不对的地方，欢迎各位大佬帮忙更正。下一次我们学习二阶或者更高阶的求解问题（默默给自己挖个坑，填不填看情况吧/笑哭）。

导言：

odeint函数介绍：

一． 定义微分方程的函数：

二． 给自变量赋值：