📜常微分方程-用例
📜Python物理量和化学量数值计算 | 📜Julia和Python蛛网图轨道图庞加莱截面曲面确定性非线性系统 | 📜Python和C++数学物理计算分形热力学静电学和波动方程 | 📜C++计算资本市场收益及成本分配数学方程 | 📜Python计算物理粒子及拉格朗日和哈密顿动力学 | 📜C代码快速傅里叶变换-分类和推理-常微分和偏微分方程 | 📜Python和C++计算物理光学波形化学结构数学方程 | 📜机器人动力学仿真 | 📜无细胞系统酶模型分析。
✒️Python自由落体运动封闭式解
常微分方程 (ODE) 是涉及函数的一些常导数(而不是偏导数)的方程。通常,我们的目标是求解 ODE,即确定哪个或哪些函数满足方程。
如果你知道函数的导数是什么,那么如何求函数本身呢?你需要找到反导数,即你需要积分。例如,如果给你
d
x
d
t
(
t
)
=
cos
t
\frac{ d x}{ d t}(t)=\cos t
dtdx(t)=cost
那么函数
x
(
t
)
x(t)
x(t) 是什么?由于
cos
t
\cos t
cost 的反导数是
sin
t
\sin t
sint,那么
x
(
t
)
x(t)
x(t) 必定是
sin
t
\sin t
sint。只是我们忘记了一个重要的一点:如果我们只知道导数,那么总是存在一个我们无法确定的任意常数。因此,我们从上面的方程可以确定的是
x
(
t
)
=
sin
t
+
C
x(t)=\sin t+C
x(t)=sint+C
对于某个任意常数
C
C
C。您可以验证
x
(
t
)
x(t)
x(t) 确实满足方程
d
x
d
t
=
cos
t
\frac{ d x}{ d t}=\cos t
dtdx=cost。
一般来说,求解 ODE 比简单积分更复杂。即便如此,基本原则始终是积分,因为我们需要从导数到函数。通常,困难的部分是确定我们需要进行哪些积分。
不过,让我们从更简单的开始吧。最简单的 ODE 是什么?令
x
(
t
)
x(t)
x(t) 为满足 ODE 的
t
t
t 函数:
d
x
d
t
=
0
(
1
)
\frac{ d x}{ d t}=0 \qquad(1)
dtdx=0(1)
我们可以问一些简单的问题。什么是
x
(
t
)
x(t)
x(t) ?
x
(
t
)
x(t)
x(t) 是由该方程唯一确定的吗?如果不是,还需要说明什么吗?
上述等式仅仅意味着 x ( t ) x(t) x(t)是一个常数函数, x ( t ) = C x(t)=C x(t)=C。它当然不是唯一确定的,因为如果我们只有 x x x 的导数方程,则无法指定常数 C C C。为了唯一地确定 x ( t ) x(t) x(t),必须根据函数 x ( t ) x(t) x(t) 本身提供一些附加数据。
让我们把事情变得更复杂一点。考虑方程
d
x
d
t
=
m
sin
t
+
n
t
3
(
2
)
\frac{ d x}{ d t}=m \sin t+n t^3 \qquad(2)
dtdx=msint+nt3(2)
其中
m
m
m 和
n
n
n 只是一些实数。方程 (2) 并不比方程 (1) 复杂多少,因为右侧不依赖于
x
x
x。它仅取决于
t
t
t。我们只是简单地用
t
t
t 来指定导数。解就是反导数或积分。
这次让我们以稍微不同的方式进行积分。我们将使用从时间
t
=
a
t=a
t=a 到时间
t
=
b
t=b
t=b 的定积分。使用微积分基本定理,
d
x
d
t
\frac{ d x}{ d t}
dtdx 从
a
a
a 到
b
b
b 的积分必须为
x
(
b
)
−
x
(
a
)
=
∫
a
b
d
x
d
t
d
t
=
∫
a
b
(
m
sin
t
+
n
t
3
)
d
t
=
−
m
cos
b
+
n
b
4
/
4
−
(
−
m
cos
a
+
n
a
4
/
4
)
.
\begin{aligned} x(b)-x(a) & =\int_a^b \frac{ d x}{ d t} d t \\ & =\int_a^b\left(m \sin t+n t^3\right) d t \\ & =-m \cos b+n b^4 / 4-\left(-m \cos a+n a^4 / 4\right) . \end{aligned}
x(b)−x(a)=∫abdtdxdt=∫ab(msint+nt3)dt=−mcosb+nb4/4−(−mcosa+na4/4).
我们可以用不同的方式编写解。我们可以用任意时间
t
t
t 替换
b
b
b,
x
(
t
)
=
−
m
cos
t
+
n
t
4
/
4
+
m
cos
a
−
n
a
4
/
4
+
x
(
a
)
.
x(t)=-m \cos t+n t^4 / 4+m \cos a-n a^4 / 4+x(a) .
x(t)=−mcost+nt4/4+mcosa−na4/4+x(a).
这种形式非常明显地表明解
x
(
t
)
x(t)
x(t) 如何依赖于初始条件
x
(
t
0
)
=
x
0
x\left(t_0\right)=x_0
x(t0)=x0。如果
x
(
7
)
=
5
x(7)=5
x(7)=5,那么
x
(
t
)
=
−
m
cos
t
+
n
t
4
/
4
+
m
cos
7
−
n
7
4
/
4
+
5
x(t)=-m \cos t+n t^4 / 4+m \cos 7-n 7^4 / 4+5
x(t)=−mcost+nt4/4+mcos7−n74/4+5
另一方面,如果我们不关心常数的形式,我们可以将通解写为
x
(
t
)
=
−
m
cos
t
+
n
t
4
/
4
+
C
x(t)=-m \cos t+n t^4 / 4+C
x(t)=−mcost+nt4/4+C
到目前为止,我们看到的 ODE 示例可以简单地通过积分来求解。它们如此简单的原因是
d
x
d
t
\frac{ d x}{ d t}
dtdx 的方程不依赖于函数
x
(
t
)
x(t)
x(t),而仅依赖于变量
t
t
t。另一方面,一旦方程同时依赖于
d
x
d
t
\frac{ d x}{ d t}
dtdx 和
x
(
t
)
x(t)
x(t),我们就需要做更多的工作来求解函数
x
(
t
)
x(t)
x(t)。
这是一个包含
x
(
t
)
x(t)
x(t) 的 ODE:
d
x
d
t
=
a
x
(
t
)
+
b
\frac{ d x}{ d t}=a x(t)+b
dtdx=ax(t)+b
💦Python自由落体示例
自由落体中的点质量 P P P。
- 重力场 g ⃗ = ( 0 , − g ) \vec{g}=(0,-g) g=(0,−g)
- 质量 m m m
- 初始位置 P 0 = ( 0 , 0 ) P_0=(0,0) P0=(0,0)
- 初始速度 V ⃗ 0 = ( v x 0 , v y 0 ) \vec{V}_0=\left(v_{x 0}, v_{y 0}\right) V0=(vx0,vy0)
问题数学表述:
{
x
¨
=
0
y
¨
=
−
g
\left\{\begin{array}{l} \ddot{x}=0 \\ \ddot{y}=-g \end{array}\right.
{x¨=0y¨=−g
封闭解:
{
x
(
t
)
=
v
x
0
t
y
(
t
)
=
−
g
t
2
2
+
v
y
0
t
\left\{\begin{array}{l} x(t)=v_{x 0} t \\ y(t)=-g \frac{t^2}{2}+v_{y 0} t \end{array}\right.
{x(t)=vx0ty(t)=−g2t2+vy0t
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
tmax = .2
t = np.linspace(0., tmax, 1000)
x0, y0 = 0., 0.
vx0, vy0 = 1., 1.
g = 10.
x = vx0 * t
y = -g * t**2/2. + vy0 * t
fig = plt.figure()
ax.set_aspect("equal")
plt.plot(x, y, label = "Exact solution")
plt.grid()
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
任何 ODE 都可以重新表述为一阶系统方程。我们假设
X
=
[
X
0
X
1
X
2
X
3
]
=
[
x
y
x
˙
y
˙
]
X=\left[\begin{array}{l} X_0 \\ X_1 \\ X_2 \\ X_3 \end{array}\right]=\left[\begin{array}{l} x \\ y \\ \dot{x} \\ \dot{y} \end{array}\right]
X=
X0X1X2X3
=
xyx˙y˙
结果为:
X
˙
=
[
x
˙
y
˙
x
¨
y
¨
]
\dot{X}=\left[\begin{array}{l} \dot{x} \\ \dot{y} \\ \ddot{x} \\ \ddot{y} \end{array}\right]
X˙=
x˙y˙x¨y¨
然后,初始二阶方程可以重新表述为:
X
˙
=
f
(
X
,
t
)
=
[
X
2
X
3
0
−
g
]
\dot{X}=f(X, t)=\left[\begin{array}{c} X_2 \\ X_3 \\ 0 \\ -g \end{array}\right]
X˙=f(X,t)=
X2X30−g
一般问题,求解
Y
˙
=
f
(
Y
,
t
)
\dot{Y}=f(Y, t)
Y˙=f(Y,t)。
一般公式:
X
˙
=
f
(
X
,
t
)
\dot{X}=f(X, t)
X˙=f(X,t)
- 近似解:需要误差估计
- 离散时间: t 0 、 t 1 、 … t _0、t _1、\ldots t0、t1、…
- 时间步 d t = t i + 1 − t i d t=t_{i+1}-t_i dt=ti+1−ti
💦欧拉法:
X
i
+
1
=
X
i
+
f
(
X
,
t
i
)
d
t
X_{i+1}=X_i+f\left(X, t_i\right) d t
Xi+1=Xi+f(X,ti)dt
dt = 0.02
X0 = np.array([0., 0., vx0, vy0])
nt = int(tmax/dt)
ti = np.linspace(0., nt * dt, nt)
def derivate(X, t):
return np.array([X[2], X[3], 0., -g])
def Euler(func, X0, t):
dt = t[1] - t[0]
nt = len(t)
X = np.zeros([nt, len(X0)])
X[0] = X0
for i in range(nt-1):
X[i+1] = X[i] + func(X[i], t[i]) * dt
return X
%time X_euler = Euler(derivate, X0, ti)
x_euler, y_euler = X_euler[:,0], X_euler[:,1]
plt.figure()
plt.plot(x, y, label = "Exact solution")
plt.plot(x_euler, y_euler, "or", label = "Euler")
plt.grid()
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()