2.机电系统CAD算法基础 2.1.数值积分算法 2.1.1基本数值积分方法 连续系统的动态特性一般可由一个微分方程或一组一阶微分方程加以描述。因 此对系统进行计算机仿真,就需要研究如何利用计算机对微分方程进行求解,目前 采用的方法是应用数值积分法对微分方程求数值解,其中欧拉法、梯形法、龙格 库塔法等数值积分法是几种基本的方法 1.欧拉( Euler)法 欧拉法是一种最简单的数值积分方法,但其精度差,实际应用很少。然而由于 其算法简单且具有明显的几何意义,有利于初学者学习应用数值积分法求解微分方 程,因此在讨论数值积分法时通常先介绍欧拉法。 假设一阶微分方程为 dy dt 初始条件为 y(to)=yo 欧拉法的基本思想就是将(2-1)式中的积分曲线解用直线段所组成的折线加以 近似。根据导数的定义可知,存在 dt △ 故 hny(+△)-y() (2-2) A→0 fIt,y(] 成立,当M足够小时,上式可近似表示为 y(t+△r)-y() ≈f,y( (23) △t 若令M=h,t=lo则得到
24 2．机电系统 CAD 算法基础 2．1．数值积分算法 2．1．1 基本数值积分方法 连续系统的动态特性一般可由一个微分方程或一组一阶微分方程加以描述。因 此对系统进行计算机仿真，就需要研究如何利用计算机对微分方程进行求解，目前 采用的方法是应用数值积分法对微分方程求数值解，其中欧拉法、梯形法、龙格一 库塔法等数值积分法是几种基本的方法。 1．欧拉（Euler）法 欧拉法是一种最简单的数值积分方法，但其精度差，实际应用很少。然而由于 其算法简单且具有明显的几何意义，有利于初学者学习应用数值积分法求解微分方 程，因此在讨论数值积分法时通常先介绍欧拉法。 假设一阶微分方程为 f (t y) dt dy = , （2-1） 初始条件为 ( ) 0 0 y t = y 欧拉法的基本思想就是将（2-1）式中的积分曲线解用直线段所组成的折线加以 近似。根据导数的定义可知，存在 ( ) ( ) t y t t y t dt dy t + − = →0 lim 故 ( ) ( ) t y t t y t t + − →0 lim = f[t, y(t)] （2-2） 成立，当 t 足够小时，上式可近似表示为 ( ) ( ) f [t, y(t)] t y t t y t + − （2-3） 若令 t = h ， 0 t = t 则得到
y(n+h)≈y(a)+h/[n,yn 其中h称为计算步长。 因为y(tn)=y,tn+h=t1,所以有y(t1)≈yo+h/(a 由此可得到所欲求解在t=1处的近似值y1为 y1=y()+h/[n,y( (2-5) 这样就可由已知点(n2y)去确定它在1y平面上的另一点(t1y1)若设P0点是初始 给定点,且微分方程的精确解y(0经过P点,对应1=1+h的y()的精确解为 y(t1),即P点,而由(2-5)式求得的y1值所对应的点为P1点,其坐标为(1,y1) 它是从P点作斜率为f(n,y)的线段与直线【=1的交点。 若从P{(12y1)点出发,可求解在t=12=n+2h上的近似值y2为 y2=y1 (2-6) 即从P1几点作斜率等于f(1,y1)的线段与直线t=t2相交而得到P2点。根据类似的 思想由式(2-6)不难推出 Vk+I=y Kr) P〔! 32 3(r;) 图2-1欧拉数值积分方法
25 ( ) ( ) ( ) 0 0 0 0 y t + h y t + hf t , y t （2-4） 其中 h 称为计算步长。 因为 ( ) 0 0 y t = y ， 0 1 t + h = t ，所以有 ( ) ( ) 1 0 0 0 y t y + hf t , y 由此可得到所欲求解在 1 t = t 处的近似值 1 y 为 ( ) ( ) 1 0 0 0 y = y t + hf t , y t （2-5） 这样就可由已知点 ( ) 0 0 t , y 去确定它在 t-y 平面上的另一点 ( ) 1 1 t , y 。若设 P0 点是初始 给定点，且微分方程的精确解 y(t) 经过 P0 点，对应 t 1 = t 0 + h 的 y(t) 的精确解为 ( ) 1 y t ，即 P 点，而由（2-5）式求得的 1 y 值所对应的点为 P1 点，其坐标为 ( ) 1 1 t , y ， 它是从 P0 点作斜率为 ( ) 0 0 f t , y 的线段与直线 1 t = t 的交点。 若从 ( ) 1 1 1 P t , y 点出发，可求解在 t = t 2 = t 0 + 2h 上的近似值 2 y 为 ( ) 2 1 1 1 y = y + hf t , y （2-6） 即从 P1 几点作斜率等于 ( ) 1 1 f t , y 的线段与直线 2 t = t 相交而得到 P2 点。根据类似的 思想由式（2-6）不难推出 ( ) k k k k y = y + hf t , y +1 （2-7） 图 2-1 欧拉数值积分方法
利用上式不断对微分方程求近似解,直至t=1点P将P,P1,P,…,P诸点连 成一条折线,即得到所求微分方程精确解y()的近似曲线解,如图21所示。(27) 式称为欧拉法递推公式 该方法简单,计算量小,由前一点即可推出后一点值,属于单步法。因其从初 值即可开始进行递推运算,无需其他条件要求,因此能够自启动。该算法精度低 适当减小计算步长h有助于提高计算精度。 2.梯形法 1).欧拉法误差原因的分析 根据欧拉法的分析可知,该方法简单但较粗 ,精度差。考虑对于微分方程 山 存在精确解 ()=y+f(,y 当t=t1时 y)=y+/(y (2-8) 由前述推导可知 y(t1)≈yn+(1-0)/(n,y) (2-9) 对照(2-8)式与(2-9)式得到,在欧拉法中存在 广(y知≈(nx-4) 用∫(ta,yn)代替变量,由矩形近似曲边梯形,因此造成欧拉法误差较大
26 利用上式不断对微分方程求近似解，直至 N t = t 点 PN 。将 P P P PN , , , , 0 1 2 诸点连 成一条折线，即得到所求微分方程精确解 y(t) 的近似曲线解，如图 2-1 所示。（2-7） 式称为欧拉 法递推公式。 该方法简单，计算量小，由前一点即可推出后一点值，属于单步法。因其从初 值即可开始进行递推运算，无需其他条件要求，因此能够自启动。该算法精度低， 适当减小计算步长 h 有助于提高计算精度。 2．梯形法 1）．欧拉法误差原因的分析 根据欧拉法的分析可知，该方法简单但较粗 糙，精度差。考虑对于微分方程 f (t y) dt dy = , ( ) 0 0 y t = y 存在精确解 ( ) ( ) = + t t y t y f t y dt 0 0 , 当 1 t = t 时 ( ) ( ) = + 1 0 1 0 t t y t y f t, y dt （2-8） 由前述推导可知 ( ) ( ) ( ) 1 0 1 0 0 0 y t y + t −t f t , y （2-9） 对照（2-8）式与（2-9）式得到，在欧拉法中存在 ( ) ( )( ) 0 0 1 0 1 0 f t y dt f t y t t t t − , , 用 ( ) 0 0 f t , y 代替变量，由矩形近似曲边梯形，因此造成欧拉法误差较大
2).改进的欧拉法一一梯形法 由前述分析可以看出,进一步改进欧拉法的方法就是利用梯形面积代替矩形面 积,以逼近曲边梯形面积,即 广(,y)≈ h 2 推而广之,得到梯形法公式为 Yl=y +l(,y)+frui (2-10) 由上式可以看出,在计算yx+1时,在等式右边存在有yxn1,因此无法求出。在 这种情况下,一般采用欧拉法进行一次选代,求出yx1的预估值,再将其代人梯形 公式计算出yx+1真值,从而构成预估校正梯形公式为 yo=y +hf(,y) (2-11) 离+1 (2-12) (2-11)式称为预报式,(2-12)式称为校正式。 3.龙格一库塔( Runge-kutt)法 考虑如下一阶微分方程 dy dt f(,y) (a)=y 设1=0+h,这里h为计算步长,在1时刻的解为y1=y{(n)+h,可在1附 近展成泰勒级数,取展开式的前三项得到 VI=yo+ h (2-13) hf(o, yo) h2( 2 at
27 2）．改进的欧拉法——梯形法 由前述分析可以看出，进一步改进欧拉法的方法就是利用梯形面积代替矩形面 积，以逼近曲边梯形面积，即 ( ) ( ) 0 1 2 1 0 f f h f t y dt t t + , 推而广之，得到梯形法公式为 ( ) ( ) 1 1 1 2 + = + + + + f t y f t y h y y , , （2-10） 由上式可以看出，在计算 +1 y 时，在等式右边存在有 +1 y ，因此无法求出。在 这种情况下，一般采用欧拉法进行一次选代，求出 +1 y 的预估值，再将其代人梯形 公式计算出 +1 y 真值，从而构成预估校正梯形公式为 ( ) y = y + hf t , y + 0 1 （2-11） ( ) ( ) 0 1 1 1 2 + = + + + + f t y f t y h y y , , （2-12） （2-11）式称为预报式，（2-12）式称为校正式。 3．龙格一库塔（Runge- kutta）法 考虑如下一阶微分方程 f (t y) dt dy = , ( ) 0 0 y t = y 设 t 1 = t 0 + h ，这里 h 为计算步长，在 1 t 时刻的解为 y1 = y(t 0 )+ h ，可在 0 t 附 近展成泰勒级数，取展开式的前三项得到 ( ) 0 0 0 0 2 2 1 2 0 0 0 2 2 2 1 0 y y t t t t t t t f f t h f y hf t y dt d y h dt dy y y h = = = = + = + + = + + , （2-13）
假设上式的解可表示为 y1=y+(a1K1+a2K2 K,=flo, yo r(,+b, h,y,+b, k,h) (2-16) 将K2在(a,y)附近展开成泰勒级数,并取前三项得到 K2=f(,y)+42+bK,9 将(2-15)式(2-17)式代入(2-14)式中,经整理得 y,=yo +a,hf(o, yo)+arh f(o, yo )+b, h2+b,hi 将(218)式与(213)式相比较得到{2a2b1=1 b2=1 若取h==1,则得到a1=a2=2 将其代入(2-14)式、(2-15)式、(2-16) 式,得 +K K y 将(2-19)式写成递推形式 y +-(K, +K K1=f(t4,y) (2-20) 2=f(t4+h,y4+K1b) 由于(220)式只取了级数展开式的前三项,略去了h2以上的高阶无穷小,其
28 假设上式的解可表示为 y1 = y0 + (a1K1 + a2K2 )h (2-14) ( ) 1 0 0 K = f t , y (2-15) K f (t b h y b K h) 2 = 0 + 1 0 + 2 1 , (2-16) 将 K2 在 ( ) 0 0 t , y 附近展开成泰勒级数，并取前三项得到 ( ) y f b K h t f K f t y b h + 2 0 0 + 1 2 1 , (2-17) 将(2-15)式( 2-17)式代入(2-14)式中，经整理得 ( ) ( ) + = + + + y f b hK t f y1 y0 a1hf t 0 y0 a2 h f t 0 y0 b1h 2 1 , , (2-18) 将（2-18）式与（2-13）式相比较得到 = = + = 2 1 2 1 1 2 2 2 1 1 2 a b a b a a 若取 b1 = b2 = 1 ，则得到 2 1 a1 = a2 = 。将其代入（2-14）式、（2-15）式、（2-16） 式，得 ( ) ( ) ( ) = + + = = + + K f t h y K h K f t y K K h y y 2 0 0 1 1 0 0 1 0 1 2 2 , , （2-19） 将（2-19）式写成递推形式 ( ) ( ) ( ) = + + = + = + + K f t h y K h K f t y K K h y y k k k k k k 2 1 1 1 1 2 2 , , （2-20） 由于（2-20）式只取了级数展开式的前三项，略去了 2 h 以上的高阶无穷小，其