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

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

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

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

### 🎯要点

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

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

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

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

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

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

📜[物理数学波形方程：Python数值和符号算法计算及3D视图物理数学波形方程](python-shu-zhi-he-fu-hao-suan-fa-ji-suan-ji-3d-shi-tu-wu-li-shu-xue-bo-xing-fang-cheng)

📜[物理数学热力学静电学和波动方程：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)

📜[物理数学拉格朗日和哈密顿动力学：Python计算物理粒子及拉格朗日和哈密顿动力学](python-ji-suan-wu-li-li-zi-ji-la-ge-lang-ri-he-ha-mi-dun-dong-li-xue)

📜[物理数学气体动能和粒子速度：MATLAB和Python数值和符号计算可视化物理学气体动能和粒子速度](../jiao-rong-guan-tong-intersection/matlab-he-python-shu-zhi-he-fu-hao-ji-suan-ke-shi-hua-wu-li-xue-qi-ti-dong-neng-he-li-zi-su-du)

📜[物理统计推理模型：Python射频电磁肿瘤热疗数学模型和电磁爆炸性变化统计推理模型](../suan-fa-mo-xing-algorithm/python-she-pin-dian-ci-zhong-liu-re-liao-shu-xue-mo-xing-he-dian-ci-bao-zha-xing-bian-hua-tong-ji-tu)

📜[物理数学流体力学：C++风流和MATLAB | Python | CUDA 库埃特流泊肃叶流薄膜流体](c++-feng-liu-he-matlab-python-cuda-ku-ai-te-liu-bo-su-ye-liu-bao-mo-liu-ti)

📜[物理数学气体运动模型：C++ | Python气泡表面张力和预期形态及上升速度数值模型](../suan-fa-mo-xing-algorithm/c++-python-qi-pao-biao-mian-zhang-li-he-yu-qi-xing-tai-ji-shang-sheng-su-du-shu-zhi-mo-xing)

{% @mermaid/diagram content="pie title 语言分比
 "Python":90
 "MATLAB":80
 "C++":30
 "C":5" %}

{% @mermaid/diagram content="pie title 内容分比
 "物理、数学":98
 "波动":50
 "热力学、静电学":45
 "动力学":47
 "气体、液体流体力学":60
 "算法":40" %}

### 🍇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}$$：这是更新后的值，给出时间间隔结束时的加速度。这种修改后的方案称为欧拉-克罗默方法。

代码实现：

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

模拟数据

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

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

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

```python
 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');
```

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

```python
 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');
```

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