Python和MATLAB粘性力接触力动态模型半隐式欧拉算法

🎯要点

🎯运动力模型计算制作过程：🖊相机捕捉网球运动图，制定运动数学模型，数值微分运动方程 | 🖊计算运动，欧拉算法离散积分运动，欧拉-克罗默算法微分运动方程

🎯粘性力模型计算制作过程：🖊绘制雨滴、环境和坐标系下落自由体草图，使用牛顿第二定律制定运动数学模型，定义无空气阻力和恒速简化数学模型，数值模型解析

🎯固体-固体接触力模型计算制作过程：🖊实验绘制悬块弹簧形变量图，使用牛顿第二定律制定平衡模型，定义无空气阻力简化数学模型，欧拉-克罗默算法数值求解，制作三维动态可视图

🎯移动物体间复合力模型计算制作过程：🖊实验绘制电梯，物重和人自由体图，模型使用牛顿第三定律确定力对和牛顿第二定律确定图中运动体，解析无运动和相对运动简化力模型

🎯龙卷风、弹跳球和河流上船只三维运动力模型计算制作过程 | 🎯加速实验车、旋转杆、圆周运动，风中的珠子和地震期间的振荡中约束运动力模型计算制作过程

📜物理学和数学模型用例 | 本文

📜物理数学波形方程：Python数值和符号算法计算及3D视图物理数学波形方程

📜物理数学热力学静电学和波动方程：Python和C++数学物理计算分形热力学静电学和波动方程

📜物理数学拉格朗日和哈密顿动力学：Python计算物理粒子及拉格朗日和哈密顿动力学

📜物理数学气体动能和粒子速度：MATLAB和Python数值和符号计算可视化物理学气体动能和粒子速度

📜物理统计推理模型：Python射频电磁肿瘤热疗数学模型和电磁爆炸性变化统计推理模型

📜物理数学流体力学：C++风流和MATLAB | Python | CUDA 库埃特流泊肃叶流薄膜流体

📜物理数学气体运动模型：C++ | Python气泡表面张力和预期形态及上升速度数值模型

🍇Python半隐式欧拉方法

在数学中，半隐式欧拉方法也称为辛欧拉、半显式欧拉、欧拉-克罗默和牛顿-斯托默-韦莱，是欧拉方法的一种改进，用于求解哈密顿方程，哈密顿方程是经典力学中出现的常微分方程组。半隐式欧拉方法是一种辛积分器，因此比标准欧拉方法能得到更好的结果。

半隐式欧拉方法可以应用于一对以下形式的微分方程：
$$\begin{array}{rl}& \frac{dx}{dt}=f(t,v)\\ & \frac{dv}{dt}=g(t,x)\end{array}$$
其中 $f$ 和 $g$ 是给定函数。这里，$x$和$v$可以是标量或向量。如果哈密顿量具有以下形式，则哈密顿力学中的运动方程采用这种形式
$$H=T(t,v)+V(t,x)$$
微分方程需在初始条件下求解
$$x\left(t_0\right)=x_0,v\left(t_0\right)=v_0$$

欧拉方法对于振荡系统存在一个根本问题。再看一下欧拉方法的近似，得到下一个时间间隔的位置：
$$x\left(t_i+\Delta t\right)\approx x\left(t_i\right)+v\left(t_i\right)\Delta t$$
它使用时间间隔开始时的速度值来将解逐步推向未来。

由于欧拉方法通过线性近似将解投影到未来，并假设区间开始时的导数值，因此它对于振荡函数来说不是很好。改进欧拉方法的一个聪明的想法是使用第二个方程的导数的更新值。

纯欧拉方法适用：
$$\begin{array}{rl}x\left(t_0\right)=x_0,& x_{i+1}=x_i+v_i\Delta t\\ v\left(t_0\right)=v_0,& v_{i+1}=v_i-{\omega }^2{x}_i\Delta t\end{array}$$
如果在 $v$ 的方程中您使用了刚刚计算的值 $x_{i+1}$ 会怎样？像这样：
$$\begin{array}{rl}& x\left(t_0\right)=x_0,x_{i+1}=x_i+v_i\Delta t\\ & v\left(t_0\right)=v_0,v_{i+1}=v_i-{\omega }^2{x}_{i+1}\Delta t\\ & \end{array}$$
请注意第二个方程右侧的 $x_{i+1}$：这是更新后的值，给出时间间隔结束时的加速度。这种修改后的方案称为欧拉-克罗默方法。

代码实现：

def euler_cromer(state, rhs, dt):
  
    mid_state = state + rhs(state)*dt # Euler step
    mid_derivs = rhs(mid_state)       # updated derivatives
    
    next_state = np.array([mid_state[0], state[1] + mid_derivs[1]*dt])
    
    return next_state

模拟数据

w = 2
period = 2*np.pi/w
dt = period/200  
T = 800*period  
N = round(T/dt)

print('The number of time steps is {}.'.format( N ))
print('The time increment is {}'.format( dt ))

t = np.linspace(0, T, N)

x0 = 2    
v0 = 0    

num_sol = np.zeros([N,2])
num_sol[0,0] = x0
num_sol[0,1] = v0

for i in range(N-1):
    num_sol[i+1] = euler_cromer(num_sol[i], springmass, dt)

The number of time steps is 160000. The time increment is 0.015707963267948967

首先，得到解析解。然后，您选择绘制振荡运动的前几个周期：数值和解析。

x_an = x0*np.cos(w * t)

iend = 800 
fig = plt.figure(figsize=(6,4))
plt.plot(t[:iend], num_sol[:iend, 0], linewidth=2, linestyle='--', label='Numerical solution')
plt.plot(t[:iend], x_an[:iend], linewidth=1, linestyle='-', label='Analytical solution')
plt.xlabel('Time [s]')
plt.ylabel('$x$ [m]')
plt.title('Spring-mass system, with Euler-Cromer method.\n');

该图显示，欧拉-克罗默不存在振幅增大的问题。从这个意义上讲，你应该对此感到满意。但是，如果你绘制一段较长模拟的末尾，你就会发现它确实开始偏离解析解。

istart = 400

fig = plt.figure(figsize=(6,4))

plt.plot(t[-istart:], num_sol[-istart:, 0], linewidth=2, linestyle='--', label='Numerical solution')
plt.plot(t[-istart:], x_an[-istart:], linewidth=1, linestyle='-', label='Analytical solution')
plt.xlabel('Time [s]')
plt.ylabel('$x$ [m]')
plt.title('Spring-mass system, with Euler-Cromer method. \n');

观察一段很长的运行中的最后几次振荡，即使时间增量很小，也会发现轻微的相位差。因此，尽管欧拉-克罗默方法解决了欧拉方法的一个大问题，但它仍然存在一些错误。它仍然是一阶方法！

👉参阅：计算思维 | 亚图跨际