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

<details>

<summary>关键词</summary>

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

</details>

{% embed url="https://viadean.notion.site/Python-322a1366e3ff4ce08853956a0f090d63" %}

{% content-ref url="../zhi-dian-mi-jin-brief" %}
[zhi-dian-mi-jin-brief](../zhi-dian-mi-jin-brief)
{% endcontent-ref %}

### 📜常微分方程-用例

📜[Python物理量和化学量数值计算](../jiao-rong-guan-tong-intersection/python-wu-li-liang-he-hua-xue-liang-shu-zhi-ji-suan) | 📜[Julia和Python蛛网图轨道图庞加莱截面曲面确定性非线性系统](../gong-cheng-engineering/julia-he-python-zhu-wang-tu-gui-dao-tu-pang-jia-lai-jie-mian-qu-mian-que-ding-xing-fei-xian-xing-xi) | 📜[Python和C++数学物理计算分形热力学静电学和波动方程](../gong-cheng-engineering/python-he-c++-shu-xue-wu-li-ji-suan-fen-xing-re-li-xue-jing-dian-xue-he-bo-dong-fang-cheng) | 📜[C++计算资本市场收益及成本分配数学方程](../jiao-rong-guan-tong-intersection/c++-ji-suan-zi-ben-shi-chang-shou-yi-ji-cheng-ben-fen-pei-shu-xue-fang-cheng) | 📜[Python计算物理粒子及拉格朗日和哈密顿动力学](../wu-li-xue-physics/python-ji-suan-wu-li-li-zi-ji-la-ge-lang-ri-he-ha-mi-dun-dong-li-xue) | 📜[C代码快速傅里叶变换-分类和推理-常微分和偏微分方程](c-dai-ma-kuai-su-fu-li-ye-bian-huan-fen-lei-he-tui-li-chang-wei-fen-he-pian-wei-fen-fang-cheng) | 📜[Python和C++计算物理光学波形化学结构数学方程](../wu-li-xue-physics/python-he-c++-ji-suan-wu-li-guang-xue-bo-xing-hua-xue-jie-gou-shu-xue-fang-cheng) | 📜机器人动力学仿真 | 📜无细胞系统酶模型分析。

### ✒️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.
$$

```python
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib as mpl
```

```python
 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
$$

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