高等数学--微分方程的求解（sympy库）

微分方程的意义

函数是客观事物内部联系的量化，利用函数关系可以研究客观事物的内在规律，因此，如何寻找函数关系在实践中有着极其重要的意义！
但是在很多问题当中，往往是若干“因”导致若干“果”，想要通过观察直接找出变量之间的函数关系比较困难，那么我们该怎么办呢?
我们可以通过问题中反映的情况，列出含有要找的函数及其导数关系式，这样的关系式就是微分方程
微分方程建立之后，需要对其进行研究进而找到未知函数来，这一过程就是解微分方程

引自《高等数学 第七版 同济大学出版社》

sympy求解微分方程的方式

symbols()+Eq()

求解下面微分方程
$$y^{\prime \prime }-2y^{\prime \prime }+y=0$$

import sympy as sp

x=sp.symbols('x')
y=sp.symbols('y',cls=sp.Function)

eq1 = sp.Eq(y(x).diff(x, x) - 2*y(x).diff(x) + y(x), 0)
res1=sp.dsolve(eq1,y(x))
print(res1)
>>> Eq(y(x), (C1 + C2*x)*exp(x))
1
2
3
4
5
6
7
8
9

Function()

求解下面微分方程
$$y^{\prime \prime }-2y^{\prime }+y=0,其中y(0)=0，y^{\prime }(0)=1$$

from sympy.abc import x
import sympy as sp
y=sp.Function('y')
eq=sp.diff(y(x),x,2)-2*sp.diff(y(x),x)+y(x)

con={y(0):0,sp.diff(y(x),x).subs(x,0):1}
res=sp.dsolve(eq,ics=con)
res
>>> Eq(y(x), x*exp(x))
1
2
3
4
5
6
7
8
9

solveset()

求解下面微分方程
$$y^{\prime }+y+y^2=0，其中y(0)=1$$

from sympy import *
f = symbols('f', cls=Function)
x = symbols('x')
eq = Eq(f(x).diff(x,1)+f(x)+f(x)**2, 0)
dsolve(eq, f(x))
>>>  Eq(f(x), -C1/(C1 – exp(x)))
C1 = symbols('C1')
eq1 = -C1/(C1 - exp(x))
eq2 = eq1.subs(x, 0)
solveset(eq2 - 1, C1)
>>> {1/2}

1
2
3
4
5
6
7
8
9
10
11
12

求解实际问题

知道了如何使用sympy库求解各种各样的微分方程，下面我们来解决一点现实中遇到的问题

1. 降落伞问题

某降落伞在下降过程中，受到的空气阻力与速度成正比，设降落伞的初始速度为0，求降落伞下降速度与时间的函数关系

假设变量

降落伞的速度：v(t)
重力：P，大小为mg，方向与v一致，
阻力：R，大小为kv，方向与v相反
降落伞所受到的合力为：F

建立模型

根据题目中给出的比那变量之间的函数关系可得
$$F=mg-kv$$
根据牛顿第二定律
$$\begin{gathered} F=ma， \\ t时刻降落伞的加速度a=\frac{d{s}^2}{d^{2}t}=\frac{dv}{dt} \end{gathered}$$

那么列出微分方程：
$$m\frac{dv}{dt}=mg-kv，其中v_{t=0}=0$$

求解微分方程模型

$$\frac{dv}{mg-kv}=\frac{dt}{m}$$
两端积分得到
$$\int \frac{1}{mg-kv}dv=\int \frac{1}{m}dt$$

由于mg-kv>0，则
$$-\frac{1}{k}ln(mg-kv)=\frac{t}{m}+C_1$$

即
$$mg-kv=e^{-\frac{k}{m}t-k{C}_1}$$

于是
$$v=\frac{mg}{k}(1-e^{-\frac{k}{m}t})$$

sympy求解

from sympy import *

m,g,k,t=symbols('m g k t')
v=Function('v')
eq=Eq(m*diff(v(t),t),m*g-k*v(t))
con={v(0):0}
res=dsolve(eq,ics=con)
res
>>> Eq(v(t), g*m/k - g*m*exp(-k*t/m)/k)
1
2
3
4
5
6
7
8
9

从上面我们可以看出，我们在求解微分方程时经过这么多步骤推导出的$v与t$的函数关系式，被sympy的几行代码就解决了！不得不叹服sympy库的便利

2. 放射性元素衰变问题

镭元素的衰变规律如下，衰变速度与其现存量R成正比现在已知1600年后，镭的衰变变为原来质量的$\frac{1}{2}$，求镭元素现存量R与时间t之间的函数关系

假设变量

现存量R
衰变常数$\lambda$
衰变速度v

构建与求解模型

根据题目中的已知关系

$$v=\frac{dR}{dt}=-\lambda R$$
那么
$$\frac{dR}{R}=-\lambda dt$$
两边积分得到
$$\int \frac{1}{R}dR=\int -\lambda dt$$
即
$$\begin{gathered} lnR=-\lambda t+lnC \\ R=C{e}^{-\lambda t} \end{gathered}$$
由于
$$\begin{gathered} t=0时，R=R_0 \\ t=1600,R=\frac{1}{2}{R}_0 \end{gathered}$$
代入上式可得
$$\lambda =\frac{ln2}{1600}=0.000433$$

于是
$$R=R_0{e}^{-0.000433t}$$

sympy求解

from sympy import *
t, k ,R0= symbols('t,k,R0')
R = Function('R')
R_ = Derivative(R(t),t)
Eqn = Eq(R_,-k*R(t))
res = dsolve(Eqn,ics={R(0):R0})
print(res)
>>> Eq(R(t), R0*exp(-k*t))
kk=solve(res,k)
k0=kk[0].subs({t:1600,R(t):0.5*R0})# 求式子中的参数k（元素衰变系数）
print(kk,k0)
>>> [log(R0/R(t))/t] 0.000433216987849966
k1=res.args[1]
k11=k1.subs({t:99999,k:k0})
print(k11)
>>> 1.53395780934154e-19*R0   # 99999年之后元素的剩余量
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16