🍠Python机器人动力学和细胞酶常微分方程

关键词

Python | 物理 | 化学 | 数值 | 算法 | 模型 | Julia | 非线性 | C++ | 数学 | 分形 | 热力学 | 静电学 | 波动 | 方程 | 资本 | 成本 | 收益 | 粒子 | C | 傅里叶 | 推理 | 分类 | 常微分 | 光学 | 机器人 | 动力学 | 细胞 | 酶

Python机器人动力学和细胞酶常微分方程 | Notion亚图跨际 on Notion

🏈指点迷津 | Brief

📜常微分方程-用例

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

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

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

如果你知道函数的导数是什么，那么如何求函数本身呢？你需要找到反导数，即你需要积分。例如，如果给你

$$\frac{ d x}{ d t}(t)=\cos t$$

那么函数 x(t) 是什么？由于\cos t 的反导数是\sin t，那么x(t) 必定是\sin t。只是我们忘记了一个重要的一点：如果我们只知道导数，那么总是存在一个我们无法确定的任意常数。因此，我们从上面的方程可以确定的是

$$x(t)=\sin t+C$$

对于某个任意常数 C。您可以验证 $x(t)$确实满足方程 $\frac{ d x}{ d t}=\cos t$。

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

不过，让我们从更简单的开始吧。最简单的 ODE 是什么？令 x(t) 为满足 ODE 的 t 函数：

$$\frac{ d x}{ d t}=0 \qquad(1)$$

我们可以问一些简单的问题。什么是$x(t)$？ $x(t)$ 是由该方程唯一确定的吗？如果不是，还需要说明什么吗？

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

让我们把事情变得更复杂一点。考虑方程

$$\frac{ d x}{ d t}=m \sin t+n t^3 \qquad(2)$$

其中 $m$和 $n$ 只是一些实数。方程 (2) 并不比方程 (1) 复杂多少，因为右侧不依赖于 $x$。它仅取决于$t$。我们只是简单地用 $t$来指定导数。解就是反导数或积分。

这次让我们以稍微不同的方式进行积分。我们将使用从时间 $t=a$ 到时间 $t=b$的定积分。使用微积分基本定理，$\frac{ d x}{ d t}$ 从 $a$ 到 $b$ 的积分必须为

$$\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}$$

我们可以用不同的方式编写解。我们可以用任意时间 $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{ d x}{ d t}$的方程不依赖于函数 $x(t)$，而仅依赖于变量$t$。另一方面，一旦方程同时依赖于 $\frac{ d x}{ d t}$和 $x(t)$，我们就需要做更多的工作来求解函数 $x(t)$。

这是一个包含 $x(t)$的 ODE：

$$\frac{ d x}{ d t}=a x(t)+b$$

💦Python自由落体示例

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

• 重力场 $\vec{g}=(0,-g)$

• 质量 $m$

• 初始位置$P_0=(0,0)$

• 初始速度 $\vec{V}_0=\left(v_{x 0}, v_{y 0}\right)$

问题数学表述：

$$\left\{\begin{array}{l} \ddot{x}=0 \\ \ddot{y}=-g \end{array}\right.$$

封闭解：

$$\left\{\begin{array}{l} x(t)=v_{x 0} t \\ y(t)=-g \frac{t^2}{2}+v_{y 0} t \end{array}\right.$$

 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、\ldots$

• 时间步 $d t=t_{i+1}-t_i$

💦欧拉法：

$$X_{i+1}=X_i+f\left(X, t_i\right) d t$$

 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()