显式积分,隐式积分和弹簧质点系统(详细公式推导和太极源码)

本人初学者,零基础入门(大二数学基础),因此本教程还算比较舒适,但是也免不了有错误,还请批评指正。


数值积分

数值积分,是用于求定积分的近似值的一种方法。在数学分析中,有很多计算给定函数的定积分是不可行的,而数值积分是利用黎曼积分等数学定义,用数值逼近的方法近似计算给定的定积分值。借助计算机和编程,数值积分可以快速而有效地计算复杂的积分。


欧拉方法

欧拉方法是一种数值积分方法,又称为欧拉折线法,是用折线来逼近曲线的一种方法。

例如 d y d x = f ( x , y ) \frac{dy}{dx}=f(x,y) dxdy=f(x,y),可以转化为 y n + 1 − y n = f ( x n , y n ) h y_{n+1}-y_n=f(x_n,y_n)h yn+1yn=f(xn,yn)h,其中h则为折线的步长。

由泰勒公式 y ( x + h ) = y ( x ) + y ′ ( x ) h + o ( h ) y(x+h)=y(x)+y^{'}(x)h+o(h) y(x+h)=y(x)+y(x)h+o(h),可以看出欧拉公式实际上是泰勒公式的离散形式。很显然,h越小,欧拉方法的结果越精确,h越大,结果越不精确。

但是h较大时,除了不精确之外,还会导致使用欧拉方法得到的值不收敛,即不稳定。

欧拉方法分为显式积分和隐式积分两种形式,其中显式积分条件稳定,隐式积分无条件稳定:
在这里插入图片描述


弹簧质点系统

在清楚了数值积分的解决方法之后,我们使用它来解决一个最简单的物理模拟问题——弹簧质点系统。

弹簧质点系统中主要有弹力和阻力。
在这里插入图片描述
在这里插入图片描述
总结弹力和阻力的计算公式如下:

{ f ( x a ) s = − k s x a − x b ∣ ∣ x a − x b ∣ ∣ ( ∣ ∣ x a − x b ∣ ∣ − l ) f ( x a ) d = − k d x a − x b ∣ ∣ x a − x b ∣ ∣ ( x a ′ − x b ′ ) ⋅ x a − x b ∣ ∣ x a − x b ∣ ∣ \left\{ \begin{array}{c} \LARGE f(x_a)_s=-k_s\frac{x_a-x_b}{||x_a-x_b||}(||x_a-x_b||-l) \\ \\ \LARGE f(x_a)_d=-k_d\frac{x_a-x_b}{||x_a-x_b||}(x_a^{'}-x_b^{'})·\frac{x_a-x_b}{||x_a-x_b||} \end{array} \right. f(xa)s=ksxaxbxaxb(xaxbl)f(xa)d=kdxaxbxaxb(xaxb)xaxbxaxb

显式方法

{ f t = ∑ b = 0 n f ( x a ) s + ∑ b = 0 n f ( x a ) d v t + d t = v t + d t ∗ f t m x t + d t = x t + d t ∗ v t \left\{ \begin{array}{c} \LARGE f_t=\sum_{b=0}^n f(x_a)_s+\sum_{b=0}^nf(x_a)_d \\ \\ \LARGE v_{t+dt}=v_t+dt*\frac{f_t}{m} \\ \\ \LARGE x_{t+dt}=x_t+dt*v_t \end{array} \right. ft=b=0nf(xa)s+b=0nf(xa)dvt+dt=vt+dtmftxt+dt=xt+dtvt

显式方法直接实现即可。

半隐式方法

{ v t + d t = v t + d t ∗ f t m x t + d t = x t + d t ∗ v t + d t \left\{ \begin{array}{c} \LARGE v_{t+dt}=v_t+dt*\frac{f_{t}}{m} \\ \\ \LARGE x_{t+dt}=x_t+dt*v_{t+dt} \end{array} \right. vt+dt=vt+dtmftxt+dt=xt+dtvt+dt

可以直接实现,求速度时使用显式,求位移时使用隐式。

隐式方法

{ v t + d t = v t + d t ∗ f t + d t m x t + d t = x t + d t ∗ v t + d t \left\{ \begin{array}{c} \LARGE v_{t+dt}=v_t+dt*\frac{f_{t+dt}}{m} \\ \\ \LARGE x_{t+dt}=x_t+dt*v_{t+dt} \end{array} \right. vt+dt=vt+dtmft+dtxt+dt=xt+dtvt+dt

隐式方法需要进行特殊的推导。

由泰勒公式的一阶展开,可以得到f的近似:

f t + d t = f t + ∂ f ∂ x △ x + ∂ f ∂ v △ v \LARGE f_{t+dt}=f_t+\frac{∂f}{∂x}△x+\frac{∂f}{∂v}△v ft+dt=ft+xfx+vfv

代入上式,可以求得△v的表示形式为:

△ v = v t + d t − v t = d t ∗ f t + d t m \LARGE △v=v_{t+dt}-v_t=dt*\frac{f_{t+dt}}{m} v=vt+dtvt=dtmft+dt

△ v = d t m ∗ ( f t + ∂ f ∂ x △ x + ∂ f ∂ v △ v ) \LARGE △v=\frac{dt}{m}*(f_t+\frac{∂f}{∂x}△x+\frac{∂f}{∂v}△v) v=mdt(ft+xfx+vfv)

将△x也表示为△v的式子:

△ x = x t + d t − x t = d t ∗ v t + d t = d t ∗ ( v t + △ v ) \LARGE △x=x_{t+dt}-x_t=dt*v_{t+dt}=dt*(v_t+△v) x=xt+dtxt=dtvt+dt=dt(vt+v)

代入上式,消去△x:

△ v = d t m ∗ ( f t + d t ∗ ∂ f ∂ x ( v t + △ v ) + ∂ f ∂ v △ v ) \LARGE △v=\frac{dt}{m}*(f_t+dt*\frac{∂f}{∂x}(v_t+△v)+\frac{∂f}{∂v}△v) v=mdt(ft+dtxf(vt+v)+vfv)

展开括号:

△ v = d t m ∗ f t + d t 2 m ∂ f ∂ x v t + d t 2 m ∂ f ∂ x △ v + d t m ∂ f ∂ v △ v \LARGE △v=\frac{dt}{m}*f_t+\frac{dt^2}{m}\frac{∂f}{∂x}v_t+\frac{dt^2}{m}\frac{∂f}{∂x}△v+\frac{dt}{m}\frac{∂f}{∂v}△v v=mdtft+mdt2xfvt+mdt2xfv+mdtvfv

移项:

△ v − d t 2 m ∂ f ∂ x △ v − d t m ∂ f ∂ v △ v = d t m ∗ f t + d t 2 m ∂ f ∂ x v t \LARGE △v-\frac{dt^2}{m}\frac{∂f}{∂x}△v-\frac{dt}{m}\frac{∂f}{∂v}△v=\frac{dt}{m}*f_t+\frac{dt^2}{m}\frac{∂f}{∂x}v_t vmdt2xfvmdtvfv=mdtft+mdt2xfvt

整理:

( 1 − d t 2 m ∂ f ∂ x − d t m ∂ f ∂ v ) △ v = d t m ( f t + d t ∗ ∂ f ∂ x v t ) \LARGE (1-\frac{dt^2}{m}\frac{∂f}{∂x}-\frac{dt}{m}\frac{∂f}{∂v})△v=\frac{dt}{m}(f_t+dt*\frac{∂f}{∂x}v_t) (1mdt2xfmdtvf)v=mdt(ft+dtxfvt)

别忘了我们操作的是一个矩阵:

( I − d t 2 M ∂ f ∂ x − d t M ∂ f ∂ v ) △ v = d t M ( f t + d t ∗ ∂ f ∂ x v t ) \LARGE (I-\frac{dt^2}{M}\frac{∂f}{∂x}-\frac{dt}{M}\frac{∂f}{∂v})△v=\frac{dt}{M}(f_t+dt*\frac{∂f}{∂x}v_t) (IMdt2xfMdtvf)v=Mdt(ft+dtxfvt)

两边乘以M:

( M − d t 2 ∂ f ∂ x − d t ∂ f ∂ v ) △ v = d t ( f t + d t ∗ ∂ f ∂ x v t ) \LARGE (M-dt^2\frac{∂f}{∂x}-dt\frac{∂f}{∂v})△v=dt(f_t+dt*\frac{∂f}{∂x}v_t) (Mdt2xfdtvf)v=dt(ft+dtxfvt)

解隐式积分需要求解上述的线性方程组,并解出△v。

其中 ∂ f ∂ x \frac{∂f}{∂x} xf ∂ f ∂ v \frac{∂f}{∂v} vf的计算方法在参考文献中写的很清楚,这里直接给出结果:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述


雅克比迭代

在隐式方法中,最终需要求解的线性方程组往往是一个巨大的稀疏矩阵,因此很难通过矩阵求逆的方式求解,这里介绍最简单的迭代求解线性方程组的方法——雅克比迭代。

在这里插入图片描述
在这里插入图片描述


实例代码

弹簧质点系统的显式方法

弹簧质点系统的隐式方法

代码由python语言以及taichi框架编写而成。

单击鼠标添加质点,相邻质点自动添加弹簧,其中红色弹簧表示伸长中,绿色弹簧表示收缩中。
在这里插入图片描述


公式和求导部分参考
隐式方法推导部分参考
雅克比迭代部分参考

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 终极编程指南 设计师:CSDN官方博客 返回首页