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

54 篇文章 0 订阅
28 篇文章 0 订阅

📜常微分方程-用例

📜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+mcosana4/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+mcos7n74/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) V 0=(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)= X2X30g
一般问题,求解 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 t0t1
  • 时间步 d t = t i + 1 − t i d t=t_{i+1}-t_i dt=ti+1ti

💦欧拉法:
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()

🔗参阅一:计算思维

🔗 参阅二:亚图跨际

  • 9
    点赞
  • 24
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值