🍠Python机器人动力学和细胞酶常微分方程
📜常微分方程-用例
📜Python物理量和化学量数值计算 | 📜Julia和Python蛛网图轨道图庞加莱截面曲面确定性非线性系统 | 📜Python和C++数学物理计算分形热力学静电学和波动方程 | 📜C++计算资本市场收益及成本分配数学方程 | 📜Python计算物理粒子及拉格朗日和哈密顿动力学 | 📜C代码快速傅里叶变换-分类和推理-常微分和偏微分方程 | 📜Python和C++计算物理光学波形化学结构数学方程 | 📜机器人动力学仿真 | 📜无细胞系统酶模型分析。
✒️Python自由落体运动封闭式解
常微分方程 (ODE) 是涉及函数的一些常导数(而不是偏导数)的方程。通常,我们的目标是求解 ODE,即确定哪个或哪些函数满足方程。
如果你知道函数的导数是什么,那么如何求函数本身呢?你需要找到反导数,即你需要积分。例如,如果给你
那么函数 x(t) 是什么?由于\cos t 的反导数是\sin t,那么x(t) 必定是\sin t。只是我们忘记了一个重要的一点:如果我们只知道导数,那么总是存在一个我们无法确定的任意常数。因此,我们从上面的方程可以确定的是
对于某个任意常数 C。您可以验证 确实满足方程 。
一般来说,求解 ODE 比简单积分更复杂。即便如此,基本原则始终是积分,因为我们需要从导数到函数。通常,困难的部分是确定我们需要进行哪些积分。
不过,让我们从更简单的开始吧。最简单的 ODE 是什么?令 x(t) 为满足 ODE 的 t 函数:
我们可以问一些简单的问题。什么是? 是由该方程唯一确定的吗?如果不是,还需要说明什么吗?
上述等式仅仅意味着是一个常数函数,。它当然不是唯一确定的,因为如果我们只有 的导数方程,则无法指定常数 C。为了唯一地确定,必须根据函数 本身提供一些附加数据。
让我们把事情变得更复杂一点。考虑方程
其中 和 只是一些实数。方程 (2) 并不比方程 (1) 复杂多少,因为右侧不依赖于 。它仅取决于。我们只是简单地用 来指定导数。解就是反导数或积分。
这次让我们以稍微不同的方式进行积分。我们将使用从时间 到时间 的定积分。使用微积分基本定理, 从 到 的积分必须为
我们可以用不同的方式编写解。我们可以用任意时间 替换,
这种形式非常明显地表明解 x(t) 如何依赖于初始条件 。如果,那么
另一方面,如果我们不关心常数的形式,我们可以将通解写为
到目前为止,我们看到的 ODE 示例可以简单地通过积分来求解。它们如此简单的原因是 的方程不依赖于函数 ,而仅依赖于变量。另一方面,一旦方程同时依赖于 和 ,我们就需要做更多的工作来求解函数 。
这是一个包含 的 ODE:
💦Python自由落体示例
自由落体中的点质量 。
重力场
质量
初始位置
初始速度
问题数学表述:
封闭解:
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 都可以重新表述为一阶系统方程。我们假设
结果为:
然后,初始二阶方程可以重新表述为:
一般问题,求解。
一般公式:
近似解:需要误差估计
离散时间:
时间步
💦欧拉法:
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()
Last updated
Was this helpful?