微分方程建模方法

一、微分方程模型

在实际问题中经常需要寻求某个变量y随另一变量t的变化规律,这个函数关系常常不能直接求出。然而有时容易建立包含变量及导数在内的关系式,即建立变量能满足的微分方程,从而通过求解微分方程对所研究的问题进行解释说明。

二、建模步骤

建立微分方程模型一般可分成以下三步:
(1)根据实际要求确定研究的量(自变量、未知函数、必要的参数等),并确定坐标系
(2)找出这些量所满足的基本规律
(3)运用这些规律列出方程和定解条件

三、常用的建模方法

下面通过实例介绍几类常用的利用微分方程建立数学模型的方法。

1、按规律直接列方程

例:建立物体冷却过程的数学模型
将某物体置入空气中,在时刻t=0测量得它的温度为u0=150度,10min后测量得它得温度为u1=100度。要求建立此物体得温度u和时间t得关系,并计算20min后物体得温度。其中假设空气的温度保持为ue=24度。
解:牛顿冷却定律是物体温度高于周围环境得物体向周围媒质传递热量逐渐冷却时所遵循得规律:当物体表面与周围存在温度差时,单位时间从单位面积散失得热量与温度差成正比,比例系数为热传递系数。
假设该物体在时刻t得温度为u=u(t),则由牛顿冷却定律,得到
$$\frac{du}{dt}=-k(u-u_{e})$$
其中,k>0。这便是物体冷却过程得数学模型。
注意到ue=24为常熟,u-ue>0,将方程改写为
$$\frac{d(u-24)}{u-24}=-kdt$$
两边积分得到
$$\int_{150}^{u}\frac{d(u-24)}{u-24}=\int_{0}^{t}-kdt$$
化简得
$$u=24+126e^{-kt}$$
把t=10,u=u1=100代入,得k=0.0506。故,20分钟后物体的温度为69.9413度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import sympy as sp
sp.var('t, k') # 定义符号变量t,k
u = sp.var('u', cls=sp.Function) # 定义符号函数
eq = sp.diff(u(t), t) + k * (u(t) - 24) # 定义方程
uu = sp.dsolve(eq, ics={u(0): 150}) # 求微分方程的符号解
print(uu)

kk = sp.solve(uu, k) # kk返回值是列表,可能有多个解
k0 = kk[0].subs({t: 10.0, u(t): 100.0})
print(kk, '\t', k0)

u1 = uu.args[1] # 提出符号表达式
u0 = u1.subs({t: 20, k: k0}) # 代入具体值
print("20分钟后的温度为:", u0)

2、微元分析法

该方法的基本思想是通过分析研究对象的有关变量在一个很短时间内的变化规律,寻找一些微元之间的关系式。

例:有高为1m的半球形容器,水从它的底部小孔流出。小孔横截面积为1平方厘米。开始时容器内盛满了水,求水从小孔流出过程中容器里水面的高度h(水面与孔口中心的距离)随时间t变化的规律。
解:由水力学知,水从孔口流出的流量Q为”通过孔口横截面的水的体积V对时间t的变化率”,满足
$$Q=\frac{dV}{dt}=0.62S\sqrt{2gh}$$
其中,0.62为流量系数,g为重力加速度(取9.8m/s2),S为孔口横截面积(单位:m2),h为t时刻水面高度(单位:cm)。
当S=1cm2=0.0001m2
$$dV=0.000062S\sqrt{2gh}dt$$
在微小时间间隔[t,t+dt]内,水面高度由h降低到h+dh(这里dh<0),容器里水的体积改变量近似为
$$dV=-\pi r^{2}dh$$
r为t时刻的水面半径,右端置符号是尾音dh<0,这里
$$r^{2}=1^{2}-(1-h)^{2}=2h-h^{2}$$

$$0.000062S\sqrt{2gh}dt=\pi (h^{2}-2h)dh$$
考虑初始条件,得到如下微分方程
$$
\begin{cases}
\frac{dt}{dh}=\frac{10000\pi }{0.62\sqrt{2g}}(h^{\frac{3}{2}}-2h^{\frac{1}{2}}) \\
t(1)=0
\end{cases}
$$
利用分离变量法,可以求得微分方程得解为
t(h)=-15260.5042h3/2+4578.1513h5/2+10682.3530
上式表达了水从小孔流出得过程中容器内水面高度h和时间t之间得关系。

1
2
3
4
5
6
7
8
import sympy as sp
sp.var('h') # 定义符号变量
sp.var('t', cls=sp.Function) # 定义符号函数
g = 9.8
eq = t(h).diff(h) -10000*sp.pi/0.62/sp.sqrt(2*g)*(h**(3/2)-2*h**(1/2)) # 定义方程
t = sp.dsolve(eq, ics={t(1): 0}) # 求微分方程的符号解
t = sp.simplify(t)
print(t.args[1].n(9))

3、模拟近似法

该方法得基本思想是在不同的假设下模拟实际的现象,即建立模拟近似的微分方程,从数学上求解或分析解的性质,再去和实际情况作对比,观察这个模型能否模拟、近似某些实际的现象。

例:(交通管理问题)在交通十字路口,都会设置红绿灯。为了让那些正形势在交叉路口或离交叉路口太近而无法停下的车辆通过路口,红绿灯转换中间还要亮起一段时间的黄灯。那么,黄灯应亮多长时间才最合理呢?
解:记v0是法定速度,刹车速度就是从开始刹车到速度v=0时汽车驶过的距离。设W为汽车的重量,μ为摩擦系数。显然,地面对汽车的摩擦力为μW,其方向与运动方向相反。由牛顿第二定律,汽车在停车过程中,行驶的距离x与时间t的关系可由下面的微分方程表示
$$
\frac{W}{g}\cdot \frac{d^{2}x}{dt^{2}}=-μW
$$
其中,g为重力加速度。化简得
$$
\frac{d^{2}x}{dt^{2}}=-μg
$$
再考虑初始条件,建立如下得二阶微分方程模型
$$
\begin{cases}
\frac{d^{2}x}{dt^{2}}=-μg \\
x\vert_{t=0}=0,\frac{dx}{dt}\vert_{t=0}=v_{0}
\end{cases}
$$
对等式两边从0到t积分,利用初始条件得
$$
\frac{dx}{dt}=-μgt+v_{0}
$$
再积分一次,求得解为
$$x(t) = -\frac{1}{2}μgt^{2}+v_{0}t$$
令$\frac{dx}{dt}=0$,可得刹车所用时间$t_{0}=\frac{v_{0}}{μg}$,从而得到刹车距离$x(t_{0})=\frac{v_{0}^{2}}{2μg}$
下面计算黄灯状态的时间T,则
$$T= \frac{x(t_{0})+I+L}{v_{0}}+T_{0}$$
其中T0是驾驶员的反应时间,代入$x(t_{0})$得
$$T= \frac{v_{0}}{2μg}+\frac{I+L}{v_{0}}+T_{0}$$
设T0=1s,L=4.5m,I=9m。另外,取具有代表性得μ=0.7,当v_{0}=45km/h,65km/h以及80km/h,黄灯时间如下表所列

v0/(km/s) 45 65 80
T/s 4.58 5.95 7.00
1
2
3
4
5
6
7
8
9
10
from numpy import array
v0 = array([45, 65, 80])
T0 = 1
L = 4.5
I = 9
mu = 0.7
g = 9.8

T = v0/(2*mu*g)+(I+L)/v0+T0
print(T)
0%