2023年哈工大计算方法上机实验报告.pdf

《2023年哈工大计算方法上机实验报告.pdf》由会员分享，可在线阅读，更多相关《2023年哈工大计算方法上机实验报告.pdf（41页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、实验报告题目：非线性方程求解摘要:非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改善的Newton法。前言：(目的和意义)掌握二分法与Newton法的基本原理和应用。数学原理:对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法:二分法和Newton法。对于二分法，其数学实质就是说对于给定的待求解的方程f(x),其在 a,夕上连续，如0,且f(x)在 a,6内仅有一个实根x*,取 区 间 中 点c,若，则c恰为其根，否则根据F 2加匕)是否成立判断根在区间 a,c 和c力中的哪一个，从而得出新区间，仍称为 a,句。反

2、复运营计算，直至满足精度为止。这就是二分法的计算思想。N e wto n法通常预先要给出一个猜测初值x o,然后根据其迭代公式产生逼近解一 的迭代数列 h，这 就 是Newton法的思想。当m接 近x*时收敛不久，但是当xo选择不好时，也许会发散，因此初值的选取很重要。此外，若将该迭代公式改善为/(4)其 中 为规定的方程的根的重数，这就是改善的Newton法，当求解已知重数的方程的根时,在同种条件下其收敛速度要比Ne w t o n 法快的多。程序设计：本实验采用Ma t/a 匕的M 文献编写。其中待求解的方程写成力me制 n 的方式，如下func t ion y-f(x);y=-xx s

3、i n(x)；写成如上形式即可，下面给出主程序。二分法源程序：clear%给定求解区间b=1.5；a=0;%误差R=l;攵=0;%迭代次数初值whi 1 e(R 5 e-6);c=(a+b)/2;iffl2(a)f l 2(c)0;a=c;el s eb=c；endR=b-a;%求出误差k=k+J；endx=c%给出解N ewton法及改善的N e wton法源程序：c J ear%输入函数f=inpu t(,请输入需规定解函数,s)%求解f(x)的导数%改善常数或重根数m i u=2;%初始值x 0 x 0=in p u t(i npu t i n itial value x0);4=0;%

4、迭代次数m a x =700;%最大迭代次数/?=0。/($“6 5 /,80。r)；求 解 (*0),以拟定初值*0 时否就是解while(a bs(R)1 e-8)xl=x0-mi u*eval(subs(f；xO ,x)/eval(subs(df,x 0,x,);R=x 1-xO；x0=x 1；k=k+l;if(eval(s u bs(f,x 0,x)m a x；%假如迭代次数大于给定值，认为迭代不收敛，重新输入初值s s=i nput Cmayh e result i s e nor,c hoose a new x 0,y/n?2)；if s trcmp(ss,ry)x O=inp u

5、 t(i np u t i nit i al value x 0 z);k=0;elseb reaken dendendk;%给出迭代次数x=x O;%给出解结果分析和讨论：1.用二分法计算方程s i n无-=0在 1,2 内的根。(=5*10下同)计算结果为x=1.2 3 ;f(x )=-3 .3 1 1 e-0 0 7;k=l 8 ；由f (x)知结果满足规定，但迭代次数比较多，方法收敛速度比较慢。2 .用二分法计算方程/一%1 =0在 1,1.5 内的根。计算结果为x =1.8 0;f (x)=2 .8 15 e-O 0 6；k=17;由f(x)知结果满足规定，但迭代次数还是比较多。3

6、.用N e w ton法求解下列方程a)xex-1=0 x o=O.5;计算结果为x=0.5 6 7;f (x)=2.3 1 3 e-0 1 6；k=4 ；由f (x)知结果满足规定，并且又迭代次数只有4次看出收敛速度不久。b)x，x 1 0 x o=1 ；c)(x-l)2(2 x 1)=0 x o=o.4 5,X O=0 .6 5;当x o=O.4 5时，计算结果为x=0 .8 3；f(x)=-8.5 8 4 e-0 1 4 ；k =4;由f(x)知结果满足规定，并且又迭代次数只有4次看出收敛速度不久，事实上该方程的确有真解X =0 .5 o当x()=0.6 5时，计算结果为x=0.0 0

7、;f(x)=0 ;k=9 ；由f(x)知结果满足规定，事实上该方程的确有真解x=0 .5,但迭代次数增多，事实上当取x o 0.6 8时，xul,就变成了方程的另一个解,这说明N e w t o n法收敛与初值很有关系,有的时候甚至也许不收敛。4.用改善的N e w t o n法求解，有2重根，取=2(X-1)2(2X-1)=0 x()=0 .5 5;并与 3.中的 c )比较结果。当X o=O.5 5时，程序死循环，无法计算，也就是说不收敛。改=1.5时，结果收敛为x=0.8 6；f(x )=4.3 8 5 7 e-0 0 7;k=l 6 ;显然这个结果不是很好，并且也不是收敛至方程的2重根

8、上。当x o=O.8 5时,结果收敛为x=1.8 9;f(x)=2 .7 3 7 e -0 2 3；k=4；这次达成了预期的结果，这说明初值的选取很重要，直接关系到方法的收敛性，事实上直接用N e w t o n法，在给定同样的条件和精度规定下，可得其迭代次数k=1 5,这说明改善后 的Newt o n法法速度的确比较快。结论：对于二分法,只要可以保证在给定的区间内有根，使可以收敛的，当时收敛的速度和给定的区间有关，二且总体上来说速度比较慢。N ew ton法，收敛速度要比二分法快，但是最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也也许不同样，也就是结果也许不时预期需要得结果。改 善

9、 的Newt o n法求解重根问题时,假如初值不妥，也许会不收敛，这一点非常重要，当然初值合适，相同情况下其速度要比N ew ton法快得多。实验报告二题目：G auss列主元消去法摘要:求解线性方程组的方法很多，重要分为直接法和间接法。本实验运用直接法的Guass消去法，并采用选主元的方法对方程组进行求解。前言：（目的和意义）1.学习G auss消去法的原理。2.了解列主元的意义。3.拟定什么时候系数阵要选主元数学原理：由 于 一 般 线 性 方 程 在 使 用G a u s s消去法求解时，从求解的过程中可以看到，若4 =0,则必须进行行互换,才干使消去过程进行下去。有 的 时 候 即 使

10、*0,但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不对的。因此有必要进行列主元技术，以最大也许的消除这种现象。这一技术要寻找行，使得I 魔f 1=max|4i-，）|tk 1并将第，行和 第k行的元素进行互换，以使得当前的或-D的数值比0要大的多。这种列主元的消去法的重要环节如下：1.消元过程对 仁1,2,，27-1,进行如下环节。1）选主元，记Iark|=maxkli k 1 1若14rttl很小，这说明方程的系数矩阵严重病态，给出警告,提醒结果也许不对。2）互换增广阵A 的 r,衣两行的元素。%-%（J=k，+1）3）计算消元aij=aij aikak

11、Jakk（i=k+1 j=k+1,.,+l）2.回代过程对左=，-1,1,进行如下计算Xk=1 -由凡。j=k-l至此,完毕了整个方程组的求解。程序设计：本 实 验 采 用 匕 的 M 文献编写。G au ss消去法源程序：cl e ara=J npu t（输入系数阵:n）b=input（z 输入列阵b:n）n=1 eng t h（b）；A=a bx-zero s（n,l）;%函数主体fo r k=l:n 1;%是否进行主元选取i f a b s(A(kfk)ah s(t)p=r;els e p=k;endend%互换元素i f P=k；f or q=k：n+1；s=A(k,q)；A(k,q)

12、=A(p,q);A(p,q)=s；endendend%判断系数矩阵是否奇异或病态非常严重i f a b s(A(kyk)yipu s HongdispC矩阵奇异，解也许不对的)e nd%计算消元，得三角阵fo r r=k+1:n;m=A(r,k)/A(kf k)；for q=k：H+I;A(rfq)=A C rfq)-A(k,q)*m；ende n dend%求解xx(n)-A(n,n+1)/A(nfn);for k=n-l:-l:l;s-0；for r=k+l:n;s=s+A(kfr)x(r);e ndt=(A(k,n-vl)s)x(k)=(A(kf n+1)-s)/A(k,k)en d结果

13、分析和讨论:例:求解方程53272651xyZ2234 O 其中 为一小数，当 =10-5,1()70,10-14,10-2 0 时分别10采用列主元和不列主元的G a u s s 消去法求解,并比较结果。记的心为求出的解代入方程后的最大误差，按规定，计算结果如下:当 =10一 5时,不选主元和选主元的计算结果如下，其中前一列为不选主元结果,后一列为选主元结果，下同。0.91 0.5 12.72 2.632.51 2.21Emax=9.3024e-010,0此时，由于 不是很小,机器误差就不是很大，由蜃见可以看出不选主元的计算结果精度还可以，因此此时可以考虑不选主元以减少计算量。当 =10年时

14、,不选主元和选主元的计算结果如下1.7 7 0.4 81.0 7 2.7 43.31 2.09Emax=2.66 8e-005,0此时由项山可以看出不选主元的计算精度就不好了，误差开始增大。当 =10*时,不选主元和选主元的计算结果如下1.23 1.0 01.6 62.0 03.11 0 00Emax=0.03,0此 时 由&1s可以看出，不选主元的结果应当可以说是不对的了,这是由机器误差引起的。当 =10回 时,不选主元和选主元的计算结果如下N a N 1NaN 2N a N 3Eg=NaN,0不选主元时，程序报错:W a rn i n g：Di v ide by zero.o这是由于机器计

15、算的最小精 度 为IO 所以此时的=IO-2。就 认 为 是o,故出现了错误现象。而选主元时则没有这种现象，并且由曷山可以看出选主元时的结果应当是精确解。结论：采用G a u ss消去法时，假如在消元时对角线上的元素始终较大（假 如 大 于10-5）,那么本方法不需要进行列主元计算，计算结果一般就可以达成规定，否则必须进行列主元这一步，以减少机器误差带来的影响，使方法得出的结果对的。实验报告三题目：Rung现象产生和克服摘要：由于高次多项式插值不收敛，会产生Rung e现象,本实验在给出具体的实例后，采用分段线性插值和三次样条插值的方法有效的克服了这一现象，并且还取的很好的插值效果。前言：（目

16、的和意义）1.深刻结识多项式插值的缺陷。2.明确插值的不收敛性如何克服。3.明确精度与节点和插值方法的关系。数学原理：在给定+1个节点和相应的函数值以后构造次的Lagrange插值多项式，实验结果表白（见后面的图）这种多项式并不是随着次数的升高对函数的逼近越来越好，这种现象就是Rung现象。解决Rung现象的方法通常有分段线性插值、三次样条插值等方法。分段线性插值：设 在 区 间 句 上，给定n+1个插值节点a=xox/xn=b和相应的函数值yo,y/,y”，求作一个插值函数。（x）,具有如下性质：1）。（勺）=x J=0,1,，2)0(x)在每个区间 x,刈上是线性连续函数。则插值函数0(幻

17、称为区间 a,b上相应个数据点的分段线性插值函数。三次样条插值：给定区间 a,8 一个分划/：a=xoxi V XN=b若函数S满足下列条件：1)S在每个区间上，芍 上是不高于3 次的多项式。2)S及其2 阶导数在出，句上连续。则称S 使关于分划/的三次样条函数。程序设计：本实验采用M a-M 的 M 文献编写。其中待插值的方程写成了 Actio A的方式，如下funct i on y=f(x);y=1/(1+25*x*x);写成如上形式即可，下面给出主程序L agrange插值源程序：n=input(将区间分为的等份数输入:5D ;S=/-1+2/A*/0：刀；%给定的定点,R f 为给定的

18、函数x=-1：0.01:1;f=0；fo r q=l:n+1;/=/;%求插值基函数for k=l:n+l;if k=q;I-1.*(x-s(k)./(s(q)s(k)；e 1 se1=1;enden df=f +Rf(s(q)*l；%求插值函数endplot(x,J%作出插值函数曲线grid onh o ld o n分段线性插值源程序clearn=inpu t 将区间分为的等份数输入:n);s=-l+2/n。;%给定的定点人 为给定的函数m=0;hh=O.001;f or x-1：h h：1;ff=o；fo r k=l:n/;%求插值基函数s w 1 tc h kcase 1i f xs(n

19、);l=(x-s(n)./(s(n+l)-s(n)；else1=0；e ndoth e rwi s eif x=s(k-l)&x=s(k)&x=s(k+l);l=(x s(k+1)./(s(k)s(k+1);el s eI-0；ende nde n dff=jf+Rf(s(k)*1；%求插值函数值en dm=m+1;f (m)=f f;end%作出曲线x=1:h h：1；pl o t(x,f,r)；g r id onhold o n三次样条插值源程序：(采用第一边界条件)c 1 ea rinput(，将区间分为的等份数输入:n);%插值区间a=-1;b=1;hh=0.0 0 1;%画图的步长s

20、=a+(b-a)/nO:n；%给定的定点,R f 为给定的函数%第一边界条件 R f(-1),Rf (1)v=5 0 00*1/(1+25*a*编 入 3-50/(1+25*a 八4；f o r k=1：；%取出节点间距h(k)=s(k+1)-s(k);en dfa r 攵=/;-/;%求出系数向量1 a m u da,miula(k)=h(k+l)/(h(kl)+h(k);m i u(k)=1 1 a(k);end%赋值系数矩阵Afork-1：n-1;for p=l：n-1;s w i tc h pca s e kA(k,p)=2;c a se k-1A(k,p)=mi u(p+1)；c a

21、s e k+1A(k,p)=la(p-1)；otherwis eA(k,p)=0;e n den dend%求出d 阵fo r k-1:n-1;s wit ch kc ase 1d(k)=6*f2c Cs(k)s Ck+1)s(k+2)J)mi u(k)v;case n 1d(k)=6*f2 c (s(k)s Ck+1)s(k+2)-la(k)*v;o th e rwised(k)=6*f2c(s(k)s(k+l)s(k+2);e ndend%求解M 阵M=A d；M=v;M;v；%m=0;f=0;for x-a:hh：b;ifx=a；p=l；el s ep-c e il(x-s(1)/(b-

22、a)/n);endff2 =0；加=0；ff4=0;m=m+1;ffl=l/h(p)*(s(p+l)-x)”*M(p)/6；ff2=1/h(p)*(x-s(p)-3*M(p+/)/6;ff3=(Rf(s(p+1)Rf(s(p)/h(p)-h(p)(M (p+1)-M如)/6产U-s(p)；ff4-Rf(s(p)M(p)*h(p)*h(p)/6;f(m)=ffl+ff2+ff3+ff4；e nd%作出插值图形x=a:h h:b；P lo t(xf f,k)h old ongr i d on结果分析和讨论：本实验采用函数/()=一 二 进 行 数 值 插 值，插值区间为-1,给定节点为1 +25%

23、xj=-l+j h,h=0.1J=O,-,o下 面 分 别 给 出La grange插值，三次样条插值，线性插值的函数曲线和数据表。图 中 只 标 出Lagrange插值的十次多项式的曲线,其它曲线没有标出,从数据表中可以看出具体的误差。表中，入。为Lagran g e 插值的1 0 次多项式,S/(x),S,。分别代表n=1 0,40的三次样条插值函数，Xl0(x),X40 6)分别代表=10,40的线性分段插值函数。xf(x)Lio(x)Sio(x)S40(X)X i o(x)X4o(*；-1.00 0.5 4 0.5 4 0.54 0.540.54 0.54-0.00 0.391.20

24、0.4 00.39 0.1 0 0.39-0 .00 0.411 .26 0.5 80.4 1 0.65 0.41-0.00 0.4 4 0.7 8 20.79 0.44 0.21 0.4 4-0.(X)0.76 0.7 6 0.76 0.7 6 0.76 0.76-0.00 0.78-0.234 0.44 0.7 8 0.82 0.7 8-0.0 0 0.2 I-0.22 6 0.6 6 0.21 0.8 8 0.2 1-0.00 0.4 9-0.1 8 0.49 0.490.940.4 9-0.0 00.0 0 0.0 00.000.0 00.000.00-0.0 00.1880.2 57

25、 0.113 0.1 8 80.000.188-0.()00.27 60.0 3 0.7 30 0.27 60.0000.2 76-O.(X)0.8 2 50.67 0.883 0.8250.00 00.825-0.0 00.0 00.00 0.00 0.000.0 00.00-0.000.2 4 60.37 60.6 40.246 0.000.2460.310.80 0.60().310.0 00.31-0.0 00.0 20.890.270.02 0.()00.02-0.000.0 00.000.0 00.00 0.0 00.(X)-0.0000.000.4 0 0.3 10.00 0.0

26、 00.000.000.0 00.900.280.(X)O.(X)0.00-0.()00.940.730.100.9410.0 00.94 101.001 .(X)1.00 1.001 .0 00.000.9410.7 300.9 40.000.9 410.000.00 0.900.280.0 0().0 00.0000.0 00.400.0 00.00 0.000.0 00.000.000.0 00.0 00.00 0.0 00.000.0 20.89 0.27 0.0 2 0.00 0.020.0 00.3 10.8 00.600.3 1().00 0.3 10.000.24 60.376

27、0.6 40.24 60.00 0.2 4 60.000.0 00.000.0 00.0 0 0.00 0.000.000.82 50.670.8 83 0.8250.0 00 0.8250.00 0,27 6 0.0 3 0.730 0.27 6 0.0 0 0 0.2 7 60.0 0 0,1 88 0.2 5 7 0.1 1 3 0.1 88 0.00 0.1 880.0 0 0.00 0.00 0.00 0.00 0.00 O.(X)0.0 0 0.49-0.18 0.49 0.49 0.94 0.4 90.0 0 0.2 1 -0.2 26 0.6 6 0.210.8 8 0.2 1

28、0.00 0.7 8-0.23 4 0.44 0.7 8 0.8 2 0.780.0 00.7 6 0.76 0.76 0.7 6 0.7 6 0.7 60.00 0.4 4 0.782 0.7 9 0.4 4 0.2 1 0.4 40.00 0.4 1 1.2 6 0.58 0.41 0.65 0.4 10.00 0.3 9I.20 0.40 0.39 0.1 0 0.391.0 0 0.54 0.54 0.540.54 0.54 0.5 4从以上结果可以看到，用三次样条插值和线性分段插值，不会出现多项式插值是出现的 Runge现象，插值效果明显提高。进一步说，为了提高插值精度，用三次样条插

29、值和线性分段插值是可以增长插值节点的办法来满足规定，而用多项式插值函数时,节点数的增长必然会使多项式的次数增长，这样会引起数值不稳定，所以说这两种插值要比多项式插值好的多。并且在给定节点数的条件下，三次样条插值的精度要优于线性分段插值，曲线的光滑性也要好一些。实验报告四题目：多项式最小二乘法摘要：对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观测和测量给出离散的一些点，再来求出具体的函数解析式。又由于测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。前言:(目的和意义)1 .学

30、习使用最小二成法的原理2 .了解法方程的特性数学原理：对于给定的测量数据(为，力)3=1,2,，n),设函数分布为J=o特别的，取巴(X)为多项式(pf(x)=xj(j=0,1,m)则根据最小二乘法原理,可以构造泛函 机”(。0,4 1,.4)=2(工i=l J=0令理 =()(攵=0,7,777)dak则可以得到法方程（死,仰）（%，。0）（%，。0）（九。0）（外,必）%（于*1）（。0 必）（必，。,“）_an,_、于冲m）求该解方程组，则可以得到解,勾,，耳,，因此可得到数据的最小二乘解】/（X）j=0程序设计：本实验采用M a t la 6 的M文献编写。其中多项式函数%=x，写成

31、c Z Jo 的方式,如下f u nc t i ony=fai（x,j）y=l;for i=1 :jy=x.*y;e nd写成如上形式即可，下面给出主程序。多项式最小二乘法源程序clear%给定测量数据点6 1）s=3 4 5 6 7 8 9;f=2.0 1 2.98 3.50 5.0 2 5.4 7 6.0 2 7.0 5；%计算给定的数据点的数目n=l e ngt h（f）；%给定需要拟合的数据的最高次多项式的次数m=1 0；%程序主体fo r k=0：m;g=zeros(lfm+l);fo r j=0：m；t=0;f or计算内积(fai(si),fai(si)t=t+fai(s(i),

32、j)fai(s(i),k)；endg(J+l)=t;endA(k+1,:)=g;%法方程的系数矩阵t=0;fo r /=/;计算内积(f(si),fai(si)t=t+f(i)*f a i(s(i),k)；endb(k+1J)=t;e nda-A /?%求出多项式系数x=s(1):0.01：s(n)Jy=0;fo r i=0:m;y=y+a(i+1)*fai(x,i)；e ndpl o t(x,y)%作出拟合成的多项式的曲线g rid onhold onpl O t(s,f,rx 9%在上图中标记给定的点结果分析和讨论：例用最小二乘法解决下面的实验数据.X-3寸5986fi2.012.988.

33、9 05,02L Q .46.02T 05并作出了(x)的近似分布图。分别采用一次，二次和五次多项式来拟合数据得到相应的拟合多项式为：yl=-0.386 4 3+0.82750X；y 2=-1.0302 4+1.06 8 9 3 x-0.0 2 0 23 x 2；y 5=-50.75309+5 1.5 3527x-19.65947 x 2+3.66585x3-0.32886x4+0.011 3 7x5;分别作出它们的曲线图,图中点划线为y l 曲线,实线为y 2 曲线，虚线为y 5 曲线。X，为给定的数据点。从图中可以看出并不是多项式次数越高越好，次数高了，曲线越能给定点处和实际吻合，但别的地

34、方就很差了。因此，本例选用一次和两次的多项式拟合应当就可以T o8实验报告五题目：Rom/?erg 积分法摘要：对于实际的工程积分问题，很难应用N e%。沏公式去求解。因此应用数值方法进行求解积分问题已有着很广泛的应用，本 文 基 于 积 分 法 来 解 决 一 类积分问题。前言：(目的和意义)1.理解和掌握Romberg积分法的原理；2.学会使用Rom berg积分法；3.明确Rom 6 e 吆积分法的收敛速度及应用时容易出现的问题。数学原理：考虑积分/(7)=,欲求其近似值，通常有复化的梯形公式、S加psi。公式和Cotes公式。但是给定一个精度,这些公式达成规定的速度很缓慢。如何提高收敛

35、速度，自然是人们极为关心的课题。为此，记 为 将 区 间la,b l进行等分的复化的梯形公式计算结果，记T2,&为将区间口力 进行2”等分的复化的Sim psion公式计算结果，记T3.k为将区间出,刈 进行/等分的复化的C ores公式计算结果。根据R ich ard so n 外推加速方法，可以得到收敛速度较快的R。仍 erg积分法。其具体的计算公式为：1.准备初值，计算Tu=f(a)+f(b)2.按梯形公式的递推关系，计算KT M=-1 Trj-,b ci r b a 八 .k+-r E /(a +尹”+05)3.按 Rombe r g 积分公式计算加速值Tm,k-m4 加 i T _

36、T4 1m=2,k4.精度控制。对给定的精度R,若I*T n|1)结果分析和讨论：进行具体的积分时,精度取R=le-8.1.求积分广，3公。精确解/=2 4 99 9 6 7 6.运营程序得Romberg积分法的函数表为1.Oe+007*4.0 0 3.0 0 2.0 02.0 0 2.0 0 02.00 0 0由函数表知Ro mb e r g 积分给出的结果为2.4999676*10”,与精确没有误差,精度很高。2.求 积 分 精 确 解/=/3=/.1 1。Ji X运营程序得兄。m。e r g 积分法的函数表为1.331.6 671.1671.681.031.4 6 1.591.1 11.

37、001.3 51.481.2 71.3001.261.711.49 1.251.37001.0 01.331.061 .640001 .1 31.5 01.7900001.931.46000001.1 900000从积分表中可以看出程序运营的结果为1.1 9,取 8 位有效数字，满足规定。3,求积分，叫 改。J。x直接按前面方法进行积分，会发现系统报错，出现了。为除数的现象。出现这种情况的因素就是当x=0 时,被积函数分母出现了 0,假如用一个适当的小数 （最佳不要小于程序给定的最小误差值，但是不能小于机器的最大精度）来代替，可以避免这个问题。本实验取 =R,可得函数表为:0.5 90.90

38、0.17 0.890.4 30.9460.60 0.4 6 0.950.680.92 0.3 80.2 20.260000000.180000故该函数的积分为0.18,取 8 位有效数字。4.求积分（sin x2dx本题的解析解很难给出,但运用Romberg积分可以很容易给出近似解,函数表为:0.95 0.24 0.3 2 2 0.3140.4 90.560.3 050.8 8 0.18 0.000.5 90.540.670.390.9 60.6 50.690.080.8 40.110.6 2故该函数的积分为0.6取 8 位有效数字。00002,00000000000结论:Rombe r g

39、积分通常规定被积函数在积分区间上没有奇点。如有奇点，且奇点为第一间断点，那么采用例3 的方法，还是可以求出来的，否则，必须采用其它的积分方法。当然，油e rg 积分的收敛速度还是比较快的。实验报告六题目：常微分方程初值问题的数值解法摘要：本实验重要采用经典四阶的R-K 方法和四阶/出加s 预测一校正方法来求解常微分方程的数值解。前言：(目的和意义)通过编写程序，进行上机计算，使得对常微分方程初值问题的数值解法有更深刻的理解，掌握单步法和线性多步法是如何进行实际计算的及两类方法的合用范围和优缺陷，特别是对这两类方法中最有代表性的方法:不/方法和/出侬方法及预测-校正方法有更好的理解。通过这两种方

40、法的配合使用，掌握不同的方法如何配合在一起，解决实际问题。数学原理：对于一阶常微分方程初值问题dy/djc=f(x,y)3(1)I 武与人先的数值解法是近似计算中很中的一部分。常微分方程的数值解法通常就是给出定义域上的个等距节点，求出所相应的函数值%。通常其数值解法可分为两大类：1.单步法:这类方法在计算的值时，只需要知道为,+八X”和 y.即可，就可算出。这类方法典型有欧拉法和R-K法。2.多步法:这类方法在计算上m 的值时，除了需要知道尤”+八龙 和外值外，还需要知道前“步的值。典型的方法如A加侬 法。经典的A-K法是一个四阶的方法。它最大的优点就是它是单步法，精度高，计算过程便于改变步长

41、。其缺陷也很明显，计算量大，每前进一步就要计算四次函数值外它的具体的计算公式如式(2)所示。四阶加包始预测-校正方法是一个线性多步法，它是由/曲加s显式公式回+1 =%+(/+2 K 2+2(+?)人/6&=f(x,yn)+i =f(x“+i，%+i)校正%=%+(9嘉+1 9/“-5加+八 乃/2 4修正 y“+i =cn+l-(c+1-p,”|)1 9 /2 7 0求导 力+1由于开始时无预测值和校正值可以运用，故令p o=c o=O,以后按上面进行计算。此方法的优点是可以减少计算量;缺陷是它不是自开始的，需要先知道前面的四个点的值%，3，必，必，因此不能单独使用。此外，它也不便于改变步长

42、。程序设计：本实验采用Ma t 1 a b的M 文献编写。其中待求的微分方程写成加4。的方式,如下f unc t io n yy=g(x,y);yy=-x 今 x-y*y;写成如上形式即可，下面给出主程序。经 典 四 阶 的 方 法 源 程 序clea r%步长选取h=0.1;%初始条件，即x=0时,y=l。y(l)=1;%求解区间a-0；b=2;%迭代公式f o r x=a:h：b-h;k 1=g(xfy(x-a)/h+1);%/og(x)=y(1),下同。k2=g(x+h/2,y(x a)/h+l)+h/2*k 1);k 3=g(x+h/2,y(xa)/h+1)+h/2*k2)；k4=g(

43、x+h,y(x-a)/h+I)+h*k3);y(x-a)/h+2)=y(xa)/h+1)+h*(k 1+2*k 2+2 k3+k4)/6;end四 阶 a 侬预测一校正方法源程序%步长选取h=0.1;%初始条件y(1)=1；%求解区间a=0;b=2;%应用RK迭代公式计算初始值y 0,y 1,y2,y 3for x=a:h:a+2h;k 1-g(xfy(x-a)/h+1);k2=g(x+h/2,y(x-a)/h+1)+h/2kl)；k3=g(x+h/2,y(x-a)/h+l)+h/2 k 2)；k 4=g(x+h,y(x-a)/h+1)+h*k3);y(x-a)/h+2)=y(x-a)/h+l

44、)+h*(kl+2*k2+2*k 3+k4)/6;end%应用预测校正法求解c(4)=0；%校正初值P 4 片o；%预测初值f(J)=g(a+0*h,y(l)；g(x)=y(x),且将该值存在数组/中。f(2)=g(a+l*h,y(2)；f(3)=g(a+2 h,y(3);f(4)=g(a+3*h,y(4)；for n=4：(b-a)/h;%P预测p(n+l)=y Cn)+h/24*(55*f(n)-59 f(n-l)+37*f(n-2)9*f (n-3);%M修正m(n+l)=p(n+1)-f-2 51/270*(c(n)-p(n);%E求导f (n+1)=g(a+(n/l-l)*h,m(n

45、-/1);%C校正c(n+1)=y(n)+lt/24*(9*f(n+1)+1 9 *f(n)5*f(n-l)hf(n-2)；%M修正y(n-/l)=c(n 1 9/270*(c(n+1)-p(n+l);%E求导f(n+I)=g(a-h(n-F1-1)*h,y(n+1);e n d结果分析和讨论：计算实例一：对 初 值 问 题)=K-y取 步 长h=o.i；计算在-1,0上的数值解。、y(T)=o在计算机上,运用所编的程序进行了运算，结果如下,其中第一列为在区间上的等分点,第二列为运用R-K法的计算结果，第三列为/da/zzs预测一校正法计算结果。由于本题的解析解很难求出，无法看出精度如何，为此

46、进行第二实例计算。第一实例计算结果-1.00 0 0-0.00 0.4 0 0.40-0.0 0 0.12 8 0.12 8-0.0 0 0.268 0.26 8-0.00 0.59 0.37-0.0 0 0.51 0.3 2-0.00 0.1 4 0.56-0.00 0.77 0.22-0.00 0.2 88 0.2 8 8-0.0 0 0.14 0.1 10 0.6 2 0.37计 算 实 例 二 对 初 值 问 题y 取 步 长h=o.i；计算在 o,is上的数值解。本题I y(0)=i的解析解为y=J1 +2 X。同样在计算机上,运用所编的程序进行了运算，结果如下，其中第一列为在区间上

47、的等分点,第二列为运用7?-K法的计算结果，第三列为4而m s预测-校正法计算结果，第四列为精确解。第二实例计算结果0 1.00 1.00 1.000.0 0 1.09 1.09 1.330.0 0 1.59 9 1 .5 99 1.9920.0 01.391.391.350.0 0 1.34 7 1.34 7 1.34 70.00 1.409 1.448 1.4 100.00 1.9 9 1.66 1.130.00 1.54 9 1.549 1.5490.00 1.99 1.4 2 1.7 10.00 1.2 6 1.70 1.1 51.00 1.57 1.8 6 1.881 .00 1.7 9 1.49 1 .8 31.00 1.91 1.83 1.5 81.00 1.80 1.72 1.0 31.00 1.1 2 1.08 1.7 91.00 2.27 1.24 2.OO根据计算结果，发现两种方法的结果与精确解很接近，精度均达成5位有效数字，但是H-K法运算是占用的内存要比/而勿S预测-校正法多，即其对计算机的规定要高一些，这与理论分析吻合。当然，四阶力就勿s预测-校正法的启动要由R-K法给出，即/，必，2，的值需预先知道。