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

Python | MATLAB | 物理 | 数学 | 方程 | 线性 | 非线性 | 动态 | 运动 | 力 | 算法 | 数值 | 微分方程 | 偏微分方程 | 模型 | 雨滴 | 坐标系 | 自由体草图 | 定律 | 空气阻力 | 简化模型

Python和MATLAB粘性力接触力动态模型半隐式欧拉算法 | 亚图跨际亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

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

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

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

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

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

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

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

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

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

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

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

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

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

🍇Python半隐式欧拉方法

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

半隐式欧拉方法可以应用于一对以下形式的微分方程：

$$\begin{aligned} & \frac{d x}{d t}=f(t, v) \\ & \frac{d v}{d t}=g(t, x) \end{aligned}$$

其中$f$ 和$g$ 是给定函数。这里，$x$和$v$可以是标量或向量。如果哈密顿量具有以下形式，则哈密顿力学中的运动方程采用这种形式

$$H=T(t, v)+V(t, x)$$

微分方程需在初始条件下求解

$$x\left(t_0\right)=x_0, \quad 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{aligned} 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{aligned}$$

如果在 $v$ 的方程中您使用了刚刚计算的值 $x_{i+1}$会怎样？像这样：

$$\begin{aligned} & x\left(t_0\right)=x_0, \quad x_{i+1}=x_i+v_i \Delta t \\ & v\left(t_0\right)=v_0, \quad v_{i+1}=v_i-\omega^2 x_{i+1} \Delta t \\ & \end{aligned}$$

请注意第二个方程右侧的 $x_{i+1}$：这是更新后的值，给出时间间隔结束时的加速度。这种修改后的方案称为欧拉-克罗默方法。

代码实现：

 def euler_cromer(state, rhs, dt):
   
     mid_state = state + rhs(state)*dt 
     mid_derivs = rhs(mid_state)       
     
     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');

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