📜常微分方程-用例

📜Python物理量和化学量数值计算 | 📜Julia和Python蛛网图轨道图庞加莱截面曲面确定性非线性系统 | 📜Python和C++数学物理计算分形热力学静电学和波动方程 | 📜C++计算资本市场收益及成本分配数学方程 | 📜Python计算物理粒子及拉格朗日和哈密顿动力学 | 📜C代码快速傅里叶变换-分类和推理-常微分和偏微分方程 | 📜Python和C++计算物理光学波形化学结构数学方程 | 📜机器人动力学仿真 | 📜无细胞系统酶模型分析。

✒️Python自由落体运动封闭式解

常微分方程 (ODE) 是涉及函数的一些常导数（而不是偏导数）的方程。通常，我们的目标是求解 ODE，即确定哪个或哪些函数满足方程。

如果你知道函数的导数是什么，那么如何求函数本身呢？你需要找到反导数，即你需要积分。例如，如果给你
$$\frac{dx}{dt}(t)=\cos t$$
那么函数 $x(t)$ 是什么？由于$\cos t$ 的反导数是$\sin t$，那么$x(t)$ 必定是$\sin t$。只是我们忘记了一个重要的一点：如果我们只知道导数，那么总是存在一个我们无法确定的任意常数。因此，我们从上面的方程可以确定的是
$$x(t)=\sin t+C$$
对于某个任意常数 $C$。您可以验证 $x(t)$ 确实满足方程 $\frac{dx}{dt}=\cos t$。

一般来说，求解 ODE 比简单积分更复杂。即便如此，基本原则始终是积分，因为我们需要从导数到函数。通常，困难的部分是确定我们需要进行哪些积分。

不过，让我们从更简单的开始吧。最简单的 ODE 是什么？令 $x(t)$ 为满足 ODE 的 $t$ 函数：
$$\frac{dx}{dt}=0(1)$$
我们可以问一些简单的问题。什么是 $x(t)$ ？ $x(t)$ 是由该方程唯一确定的吗？如果不是，还需要说明什么吗？

上述等式仅仅意味着$x(t)$是一个常数函数，$x(t)=C$。它当然不是唯一确定的，因为如果我们只有 $x$ 的导数方程，则无法指定常数 $C$。为了唯一地确定 $x(t)$，必须根据函数 $x(t)$ 本身提供一些附加数据。

让我们把事情变得更复杂一点。考虑方程
$$\frac{dx}{dt}=m\sin t+n{t}^3(2)$$
其中 $m$ 和 $n$ 只是一些实数。方程 (2) 并不比方程 (1) 复杂多少，因为右侧不依赖于 $x$。它仅取决于$t$。我们只是简单地用 $t$ 来指定导数。解就是反导数或积分。

这次让我们以稍微不同的方式进行积分。我们将使用从时间 $t=a$ 到时间 $t=b$ 的定积分。使用微积分基本定理，$\frac{dx}{dt}$ 从 $a$ 到 $b$ 的积分必须为
$$\begin{array}{rl}x(b)-x(a)& ={\int }_a^b\frac{dx}{dt}dt\\ & ={\int }_a^b\left(m\sin t+n{t}^3\right)dt\\ & =-m\cos b+n{b}^4/4-\left(-m\cos a+n{a}^4/4\right).\end{array}$$
我们可以用不同的方式编写解。我们可以用任意时间 $t$ 替换 $b$，
$$x(t)=-m\cos t+n{t}^4/4+m\cos a-n{a}^4/4+x(a).$$
这种形式非常明显地表明解 $x(t)$ 如何依赖于初始条件 $x\left(t_0\right)=x_0$。如果$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+C$$
到目前为止，我们看到的 ODE 示例可以简单地通过积分来求解。它们如此简单的原因是 $\frac{dx}{dt}$ 的方程不依赖于函数 $x(t)$，而仅依赖于变量 $t$。另一方面，一旦方程同时依赖于 $\frac{dx}{dt}$ 和 $x(t)$，我们就需要做更多的工作来求解函数 $x(t)$。

这是一个包含 $x(t)$ 的 ODE：
$$\frac{dx}{dt}=ax(t)+b$$
💦Python自由落体示例

自由落体中的点质量 $P$。

• 重力场 $\vec{g}=(0,-g)$
• 质量 $m$
• 初始位置$P_0=(0,0)$
• 初始速度 ${\vec{V}}_0=\left(v_{x0},v_{y0}\right)$

问题数学表述：
$$\{\begin{array}{l}\ddot{x}=0\\ \ddot{y}=-g\end{array}$$
封闭解：
$$\{\begin{array}{l}x(t)=v_{x0}t\\ y(t)=-g\frac{t^2}{2}+v_{y0}t\end{array}$$

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=\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]$$
结果为：
$$\dot{X}=\left[\begin{array}{l}\dot{x}\\ \dot{y}\\ \ddot{x}\\ \ddot{y}\end{array}\right]$$
然后，初始二阶方程可以重新表述为：
$$\dot{X}=f(X,t)=\left[\begin{array}{c}{X}_2\\ X_3\\ 0\\ -g\end{array}\right]$$
一般问题，求解$\dot{Y}=f(Y,t)$。

一般公式：
$$\dot{X}=f(X,t)$$

• 近似解：需要误差估计
• 离散时间：$t_0、t_1、\dots$
• 时间步 $dt=t_{i+1}-t_i$

💦欧拉法：
$$X_{i+1}=X_i+f\left(X,t_i\right)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] = X0for i in range(nt-1):X[i+1] = X[i] + func(X[i], t[i]) * dtreturn 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()

🔗参阅一：计算思维

🔗 参阅二：亚图跨际