常用算法-几种数字积分法

1、几种常用的数字积分方法(微分方程的数字解)25数字积分法1欧拉法(折线法)设一阶微分方程dy=y=f(t,y)y(t)二ydx00由图可知，过(t0,y0)点的斜率为y二f(t,y)如果/离t很近，即At很小，曲线y(t)可用切线来近似，其1切线方程y二y+f(t,y)(t-1)其微分方程在t=t1时，可近似表示为y(t)=y=y+f(t,y)(t-1)重复上述近似过程，当t=t时，2y(t)=y=y+f(t,y)(t-1)则有2一2般近1似公1式121y(t)二y二y+f(t,y)(t-1)如果令-h，称为计算步矩，则n+1nny(t)=y=y+f(t,y)-h(1)这就是欧拉法数字积分的递

2、推计算公式。由公式可看出，只要我们给出方程的初值(t0,y0)以及相应的步距，逐步进行递推就可获得微分方程的近似数字解。欧拉法的计算是十分简单的，其计算误差正比于h2，由此，要获得高精度解，必须减小步距，但这使得计算次数增加，又由于计算机的字长有限，h减小得过小，将引入舍入误差，所以此方法的精度提高有限，实际应用中较少采用。2梯形法(预报一一校正法)欧拉法精度低，却给我们一些启发，对微分方程y=f(t，y)可改写成y(t)=y+Itf(t,y)dt当t=时则0y(t)二y+hf(t,y(t)dt0t从此式可以看出0，要求得y(t)的值，等式右边中含有未知函数，所以不能得到y(t)的值，但如果我

3、们用已知的函数值y(t)来代替y(t)，用不变取代变化的函数，即Itif(t,y(t)dt沁Itif(t,y(t)dtt0t000实0际上右边是一个矩0形面积1tif(t,y(t)dt二f(t,y(t)(t-t)t0010则y二y+hf(t,y)递推公式为y=y+hf(t,y)用此矩形的面积的算法，其计算误差是显然的(欧拉法)，为了提高精度，我们可以用梯形面积来取代矩形的面积，即Itif(t,y(t)dt=丄(f+f)ht2010贝y二y+i(f+f)h02010递推形式为y=y+lh(f+f)n+1n2nnn+1或y二y+hf(t,y)+f(t,y)n+1n2nnnn+1n+1应用上式求积分

4、，产生了新的问题，即在计算y时，要用y，而y不知，则f(ty)是未知的，要获得y，通常可用迭代方法，即在t与J之间迭代多次，使其计算的y逐步收敛于y(t)，即n+1n+11(5)00yo二y+hf(t,y)n+1nnnny1二y+hf(t,y)+f(t,yo)n+1n2nnnn+1n+11yk二y+hf(t,y)+f(t,yk-1)n+1n2nnnn+1n+1如果序列尹极限存在，则当k”时，ykTy(t),要保证上述极限存在，只要选取h小到一定程度，就能得到满足。当选取一定的满足了收敛条件，但在计算上要迭代多次才认为求得了准确值，迭代次数越多，计算精度越高，但计算工作量加大，所以一般只迭代一次

5、即可，则算法写成(2)yo二y+hf(t,y)n+11nnnny二y+-hf(t,y)+f(t,yo)n+1n2nnnn+1n+1上式为预报校正公式。应用梯形近似进行校正求得的值，实际上此方法是将欧拉法与梯形法的结合的一种算法，计算量比欧拉法增加了一倍。3龙格库塔法(runge-kutta)将梯形法进一步扩展，可以得到经常使用的一种算法考虑如下的微分方程：TOC o 1-5 h z学二y二f(t,y)y(t)二y(3)dxoo设h为步长，即取t二t+h，贝U1oy二y(t+h)当h取值相对小时，可应用泰劳级数将y在y处展开，而保留h2项，即10yy+f(t,y)h+丄(竺+竺今)h21ooo2

6、ataydt(4)t=t0y=打在计算时，为避免求微分，我们设y=y+(ak+ak)-h101122k二f(t,y)1ook=f(t+b-h,y+b-k-h)o1o21当h很小，在(t,y)处泰勒展开，有1Qtt二t2Qyy二y代入(5)式得dfk=f(t,y)+(b200将k，k12dfdfy=y+(a+a)f(t,y)h+(ab+abf)h2101200213t22Qy(6)a+a=11b21ab=212ab=222方程中有四个未知数，则解有无穷多组，可取一个解，将(6)与(4)比较，可得选取a=a，则有：121,TOC o 1-5 h za=a=b=1，b=112212代入(5)，可得二

7、阶龙格库塔法的计算公式y=y+(k+k)h0212(rk=f(t,y)(7)100k=f(t+h,y+kh)001由于在台劳级数中只保留了h2以下的项，所以称为二阶龙格一库塔法，此法的截断误差正比于h3。比较(2)式，这组公式完全一样，其计算工作量完全相同，从而也证明了梯形法的截断误差正比于h3。如果我们要求更高的计算精度，可保留级数的h4及以下的项，其此时的截断误差正比于h5，其公式就是在仿真中用得较广泛的四阶龙格库塔法，它有多种写法，其中一种为：y=y+h(k+2k+2k+k)1061234k=f(t,y)100k=f(t+1h,y+壬hk)202021k=f(t+1h,y+1hk)020

8、22k=f(t+h,y+hk)003推导方法与二阶方法相同，(8)但比较麻烦，这里就不再推导了。采用公式(7)和(8)，在计算时只用到上一次的值y，而与更前的值无关，具有这种计算法称为一步法。在计算中，步长h是可以变化的，其变化范围是可以根据精度要求而定。前面的算法是以一阶一元微分方程而得到的，但一般系统均是高于一阶的，根据控制原理可知，高阶系统的模型可化成n元一次微分方程组的形式，这种模型结构为线性微分方程组：oiy=Ay+Buy叫=02Ly个nXn阶方阵,nB是一个nXm阶方阵,向量，n是系统的阶数，m为系统的输入个数，u为m阶输入向量，其中y=f(u(t),y)=ay+ayHbay+buHbbuiiiiii22inniiiimmi=1,2,n对这种方程求解，为使程序设计紧凑，我们引入一个向量和两个矩阵：H=0,1h,1h,h=h1,h2,h3,h4220k11k12k13k14K=k0,k1,k2,k3,k4=0k21k22k23k240kn1kn2kn3kn4+h)0t(1uMhn+2+:0t(1uh+2+: