《常用算法几种数字积分法》由会员分享，可在线阅读，更多相关《常用算法几种数字积分法（7页珍藏版）》请在金锄头文库上搜索。
1、1584.n 阶代数方程（n 大于等于 5）（搜索法求一个根）对于 n 5 的代数方程，没有求根公式可寻，要求方程的根，通常是采用根搜索方法，求得方程的某一个根，然后将此根从原方程中劈去，使方程降一阶或二阶，继续求根。寻求方程在根平面上的某一个根，其搜索方法有很多种，但大部分方法对重根或密集根得方程搜索将会失败，下面介绍一种搜索较保险并能得到一定精度的搜索方法，即牛顿下山与撒网格结合的搜索法。此方法可求得任何形式代数方程的根，设代数方程 0azza)z(f n11n0 并且不失一般性， 0设方程的一个试验根为 0yjx当在此试验根附近存在方程的一个根，则有 )dy(j)dx(zyjxz 00
2、代入方程得 )z(vu)(f如果 时，)(s222则 为方程得一个根。关键是怎样求得试验根的增量 dx, dy 值，使得 z 向方程yjxz的某一个根趋近。这就是利用牛顿下山法，其方法是将 f(z) 按泰劳级数：20000 )z(!2)fz)(f()f展开后取一阶项得： df0根据 随 z 变化的下降性来判断下山迭代是否成功，如果 )(f )z(s)0则所选择的 dz 是成功的，将 继续迭代，直到 =0 或 或 txmax 时，搜索过程结束。计算步长选择原则是，保证在一个步长区间内的根不能多于 1 个互异根（重根除外），否则将造成丢根现象，这就要求 h 取得较小，但增加计算时间，所以应综合考虑
3、，对周期函数，h 应小于周期的 1/2。控制精度 的大小直接影响求根精度，应根据需要设定，一般可取 10-5 10-10 ， 可取 10.1 。矩阵的逆矩阵 A 是 nn 阶矩阵，若 det(A)不为 0，则 A 的逆存在。这里介绍常用的高斯约当列主元消去法求矩阵的逆，其原理是 A-1A=I 相当于EmEm-1E1A I , 显然有 A-1=C= EmEm-1E1I对 A 作初等变换的同时，对单位矩阵作同样的初等变换，当 A 化为单位矩阵时，原单位矩阵变成了 A 的逆矩阵，即 I：AA -1：I在初等变换时，为防止主对角线上元素较小，造成计算机溢出，则选一列中最大的元素作为主元，并移到对角线上
4、，将主元变为 1，再将该列的其它元素变成 0。在初等变换的同时，将主元相乘可得 A 的行列式值 det，注意每交换一次行元素，其行列式值要变号。例： , 求 A-1，构造(设 det 初值为 1)30165a11a21(第一列), 主元 b =a21=10,交换 1,2 行,有行交换,det=-detb=-100主元 归化，第一行除以 10651消去第列其它元素，2 行1 行503b=主元=a 22=-9, 没有行交换，det=detb=9095.1第二行对角线元素归 1 化,消去第二列非对角元素1083则 A-1=913 8/953/由上过程，可以年到初等变换过程是：对于第 I 次变换，找主
5、元（对角线及以下的最大元素）移主元行到对角线行a ij 归 1 化消去第 I 列除对角线外的其它元素 ，即aki=aki-aji*akj几种常用的数字积分方法（微分方程的数字解）25 数字积分法1 欧拉法（折线法）161设一阶微分方程 )y,t(fdx0t(由图可知，过（t 0, y0）点的斜率为),t(fy0如果 离 很近，即 很小，曲线 y(t)可用切线来近似，其切线方程1t)t(,00其微分方程在 t=t1 时，可近似表示为 )ty,fy)t( 011重复上述近似过程，当 时，2t )t(y,tfy( 12112则有一般近似公式 )(,)(1 nnnntf如果令 ，称为计算步矩，则n1h
6、t(1)ny,tfy)t( 这就是欧拉法数字积分的递推计算公式。由公式可看出，只要我们给出方程的初值（t 0, y0）以及相应的步距，逐步进行递推就可获得微分方程的近似数字解。欧拉法的计算是十分简单的，其计算误差正比于 ，由此，要获得高精度解，必须2h减小步距，但这使得计算次数增加，又由于计算机的字长有限，h 减小得过小，将引入舍入误差，所以此方法的精度提高有限，实际应用中较少采用。2 梯形法（预报校正法）欧拉法精度低，却给我们一些启发，对微分方程 ),(ytfy可改写成 d)y,(ft0当 时，则1t10t t,)(t0 t1 ty(t)hyy(t0)y(t1)图 2-5-1t0 t1 tf
7、(t)ff(t0)图 2-5-2h162从此式可以看出，要求得 的值，等式右边中含有未知函数，所以不能得到)t(y1的值，但如果我们用已知的函数值 来代替 ，用不变取代变化的函数，)t(y1 )t(0)t(y即 1010 t 0t dtt,fdt)(y,f实际上右边是一个矩形面积 )t()t(y,tft)(,tf 0100t10 则 ,fhy01递推公式为 ),t(fnn1n用此矩形的面积的算法，其计算误差是显然的（欧拉法），为了提高精度，我们可以用梯形面积来取代矩形的面积，即 01021t h)f(dt)(y,f10 则 0211hf递推形式为 )f(1nnn 或 y,t),tfy1n应用上
8、式求积分，产生了新的问题，即在计算 时，要用 ，而 不知，1n1ny1n则 是未知的，要获得 ，通常可用迭代方法，即在 与 之间迭),t(f1n 1n t代多次，使其计算的 逐步收敛于 ，即1n)t()y,t(fhy01n )y,t(f201nnn ),t(f)y,t(fhy1knnnk1n 如果序列 极限存在，则当 时， ，要保证上述极限存在，只k)t(y1n要选取 h 小到一定程度，就能得到满足。当选取一定的满足了收敛条件，但在计算上要迭代多次才认为求得了准确值，迭代次数越多，计算精度越高，但计算工作量加大，所以一般只迭代一次即可，则算法写成(2)y,t(f),t(fh21y01nnn1n
9、0上式为预报校正公式。应用梯形近似进行校正求得的值，实际上此方法是将欧拉法与梯形法的结合的一种算法，计算量比欧拉法增加了一倍。1633 龙格库塔法(runge-kutta)将梯形法进一步扩展，可以得到经常使用的一种算法。考虑如下的微分方程：(3)y,t(fdx0t(设 h 为步长，即取 ，则 h1)t(01当 h 取值相对小时，可应用泰勒级数将 在 处展开，而保留 项，即1y02h(4)2001 th)dtft(2h)y,t(fy 在计算时，为避免求微分，我们设(5)hkby,ht(fk)t(fa1201021当 h 很小，在 处泰勒展开，有 )hkyftf(),t(f 1020102 将 ，
10、 代入（5）式得1k2(6)2212010 )fbatfa(h)y,t(f)a(y 将（6）与（4）比较，可得 21ba21方程中有四个未知数，则解有无穷多组，可取一个解，选取 ，则有：21a, , 21a1b2代入(5)，可得二阶龙格库塔法的计算公式(7)hky,t(fky102120由于在泰勒级数中只保留了 以下的项，所以称为二阶龙格库塔法，此法的截断2误差正比于 。比较(2) 式，这组公式完全一样，其计算工作量完全相同，从而也证3h明了梯形法的截断误差正比于 。3164如果我们要求更高的计算精度，可保留级数的 及以下的项，其此时的截断误差正4h比于 ，其公式就是在仿真中用得较广泛的四阶龙
11、格库塔法，它有多种写法，其5h中一种为：(8)hky,t(fk,t)(f)k2k61y304212103204301推导方法与二阶方法相同，但比较麻烦，这里就不再推导了。采用公式（7）和（8），在计算时只用到上一次的值 ，而与更前的值无关，具有0y这种计算法称为一步法。在计算中，步长 h 是可以变化的，其变化范围是可以根据精度要求而定。前面的算法是以一阶一元微分方程而得到的，但一般系统均是高于一阶的，根据控制原理可知，高阶系统的模型可化成 n 元一次微分方程组的形式，这种模型结构为线性微分方程组：BuAy n0201y)t(A 是一个 nn 阶方阵，B 是一个 nm 阶方阵，y 为 n 阶向量，n 是系统的阶数，m 为系统的输入个数，u 为 m 阶输入向量，其中 n,21iubaa)y,t(f mi1ini2i1iii 对这种方程求解，为使程序设计紧凑，我们引入一个向量和两个矩阵：H=0, , , h=h1, h2, h3, h4h21K=k0, k1, k2, k3, k4= 5n4k32n1k0V=v1, v2, v3, v4= 4m0210m210m01 )ht(u)t(u)ht(u)t( thttt 由（8）式可得计算矩阵 k 的递推公式： 4,3jvBhy(Akj1j0j hF613614321001
