《数值分析》上机实验报告

1、数值分析上机实验报告?数值分析?上机实验报告1. 用Newton法求方程X7-X4+14=0在(0.1, 1.9)中的近似根(初始近似值取为区间端点，迭代6次或误差小于0.00001)。1.1理论依据：设函数在有限区间a, b上二阶导数存在，且满足条件1. f(x)f(b) 02. f(x)在区间a, b上不变号3f(x)0;4f (c) | f .(x) |,其中 c是a,b中使 mir(| f .(a), f .(b) |)到达的一个 b a那么对任意初始近似值x0a,b,由Newton迭代过程f (xk)Xk1(Xk) Xk k ,k 0,1,2,3f'(Xk)所生的迭代序列Xk

2、平方收敛于方程f(x) 0在区间a,b上的惟一解令74f (x) x 28x14, f (0.1)0, f (1.9)0f (x) 7x6112x3 7x3(x316)0f (x)42x5336x242x2(x38)0f (1.9) f (1.9)0故以1.9为起点f(Xk)Xk 1 Xkf (Xk)x0 1.9如此一次一次的迭代，逼近X的真实根。当前后两个的差=£时， 就认为求出了近似的根。本程序用 Newt on法求代数方程(最高次数不 大于10)在(a,b)区间的根。1.2 C语言程序原代码:#i nclude<stdio.h>#in clude<math.h

3、> mai n()double x2,f,f1;double x1=1.9;/ 取初值为 1.9dox2=x1;f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3);/限制循环次数x仁 x2-f/f1;while(fabs(x1-x2)>=0.00001|x1<0.1); printf("计算结果:x=%fn",x1);1.3运行结果：1.4 MATLAB上机程序fun cti on y=Newt on( f,df,x0,eps,M) d=0;for k=1:Mif feval(df,x0)=

4、0 d=2;breakelse x1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1)v=eps d=1;breakend endif d=1y=x1;elseif d=0y='迭代M次失败'；elsey= '奇异 'endfunction y=df(x)y=7*xA6-28*4*xA3;Endfunction y=f(x)y=xA7-28*xA4+14;End>> x0=1.9;>> eps=0.00001;&g

5、t;> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论：1. 使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围 内得到比较理想的计算结果。 此程序的缺乏之处是， 所要求解的方程必须满足上 述定理的四个条件，但是第二和第四个条件在计算机上比较难以实现。2. Newton迭代法是一个二阶收敛迭代式，他的几何意义Xi+1是Xi的切线与 x轴的交点，故也称为切线法。它是平方收敛的，但它是局部收敛的，即要求初始 值与方程的根充分接近，所以在计算过程中需要先确定初

6、始值。3. 此题在理论依据局部，讨论了区间(0.1,1.9)两端点是否能作为Newton迭代 的初值，结果发现 0.1 不满足条件，而 1 . 9满足，能作为初值。另外，该程序简 单，只有一个循环，且为顺序结构，故采用 do-while 循环。当然也可以选择 for 和 while 循环。2. 函数值如下表:X12345f(x)o0.693147181.o9861231.38629441.6094378X67891of(x)1.79175951.94591012.0794452.19722462.3025851f(x)f (1)=1f (10)=0.1试用三次样条插值求f(4.563)及f (

7、4.563)的近似值2.1理论依据S(x) Mj(Xj x)31 _6hj 1(x Xj 1)3Mj L- (yj1 Mj6hj 1.2h, 1 X, x1士)(+ )6 hj 1(力.2hj 1 x Xj 1 叫士)(十-)6 hj 1这里hj 1XjXj，所以只要求出就能得出插值函数S(x)。MoM1doa这里的方法为:dNdodj_6( _hoho6 (yj1 yj hj 1 hjhj61h yN h (yN h N 1hN 1hj 1hj 1 hjyoyo)yjhj 1yN 1 )丄)(j 1,2,L ,Ni)最终归结为求解一个三对角阵的解。用追赶法解三对角阵的方法如下:idi, li

8、 ii 1bi 1li iC ,AbiG82b2OC2O3n 1O1l11l2 1O1Oln 112 2O On 1n 1nLU0 1 canbnn 1Ld卄1LUxd,即,假设记M，那么由Ld得Uxn11d111X11l21MMOOMMO OMMn 1n 1MMln1ndnnXnn综上可得求解方程Ax=d的算法:i 1 di 1 1i 1 i i 1,2,3,L ,n 1Xn, Xi丄沁，i n 1,L ,2,1ni2.2 C语言程序代码:#i nclude<stdio.h>#in clude<math.h> void mai n()int i,j,m, n,k,p;

9、double q10,p10,s4,g4,x0,x1,g0=1,g9=0.1;double s1010;double a10,b10,c10,d10,e10,x10,h9,u9,r9;double f10=0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851; printf("请依次输入 xi:n");for(i=0;i<=9;i+)scanf("%lf",&ei);/求 h 矩阵for(n=0;n<=8;

10、n+)h n=e n+1-e n;d0=6*(f1-f0)/h0-g0)/h0;d9=6*(g9-(f9-f8)/h8)/h8;for(j=0;j<=7;j+) dj+1=6*(fj+2-fj+1)/hj+1-(fj+1-fj)/hj)/(hj+hj+1);for(m=1;m<=8;m+) um=hm-1/(hm-1+hm);for(k=1;k<=8;k+)rk=hk/(hk-1+hk);for(i=0;i<=9;i+) / 求 u 矩阵 for(p=0;p<=9;p+)sip=0;if(i=p)sip=2;s01=1;s98=1;for(i=1;i<=8;

11、i+)sii-1=ui;sii+1=ri;printf(" 三对角矩阵为 :n");for(i=0;i<=9;i+)for(p=0;p<=9;p+)/求 r 矩阵 printf("%5.2lf",sip);if(p=9)printf("n");printf(" 根据追赶法解三对角矩阵得： n");a0=s00;b0=d0;for(i=1;i<9;i+)ci=sii-1/ai-1;/求 d 矩阵ai=sii-si-1i*ci;bi=di-ci*bi-1;if(i=8)p10=bi;q10=ai; x

12、9=p10/q10;printf("M10=%lfn",x9);for(i=9;i>=1;i-) xi-1=(bi-1-si-1i*xi)/ai-1;printf("M%d=%lfn",i,xi-1);printf("可得s(x)在区间4,5上的表达式；n");printf("将 x=4.563 代入得：n");x0=5-4.563;x1=4.563-4;s4=x3*pow(x0,3)/6+x4*pow(x1,3)/6+(f3-x3/6)*(5-4.563)+(f4-x4/6)*(4.563 -4);g4=-

13、x3*pow(x0,2)/2+x4*pow(x1,2)/2-(f3-x3/6)+(f4-x4/6);printf("计算结果：f(4.563)的函数值是：lfnf(4.563)的导数值是：lfn",s4,g4);2.3 运行结果：2.4 问题讨论1. 三次样条插值效果比Lagrange插值好，没有Runge现象，光滑性较好。2. 此题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。3. 编程过程中由于定义的数组比较多，需要仔细弄清楚各数组所代表的参 数，要注意各下标代表的含义，特别是在用追赶法计算的过程中 。33.用 Romberg算法求 3xx1'4(5

14、x 7)sinx2dx(允许误差0.00001).3.1理论依据：Romberg算法的计算步骤如下:(1) 先求出按梯形公式所得的积分值f(a)f(b)(2) 把区间2等分，求出两个小梯形面积之和，记为 T,，即久晋f(a) f(b) 2f (穿)这样由外推法可得T2(0) , t2(0)(1)扌聲1 (；)24T, T,(0)n(3) 把区间再等分(即22等分)，得复化梯形公式T,，由T,与T外推可(2) (1) 2 (1) (0)得T2 红 匚，T3(0) 4 T22 T2 ，如此，假设已算出2k等分的复化梯形公式4 141TV，那么由Richards on外推法，构造新序列4mT (k)

15、 Tk 1)Tnmk11)4 Tmm Tm ，m=1,2,1, k=1,2,l-m+141最后求得Tl(0)(4) T|(0)丁黑或 |T|(0)Ti(01 |就停止计算，否那么回到(3)，计算T1(l 1)，般可用如下算法:T1(0)T1(l)(k 1)T m 1b a f(a)21 (i 1) b 尹 P i14mT,k) T,k1)f(b)211afa(2i1,2,L ,l , k 1,2,L ,l m 1其具体流程如下，并全部存入第一列通常计算时，用固定l=N来计算，一般1=4或5即能到达要求。3.2 C语言程序代码：#in clude<math.h>#i nclude&l

16、t;stdio.h>double f(double x)/计算 f(x)的值double z;z=pow(3,x)*pow(x,1.4)*(5*x+7)*si n(x*x);return(z);mai n() double t2020,s,e=0.00001,a=1,b=3;下为的差是int i,j,l,k;t01=(b-a)*(f(b)+f(a)/2;/romberg 算法t11=(b-a)*(f(b)+2*f(b+a)/2)+f(a)/4; t02=(a*t11-t01)/(4-1);j=3;for(l=2;fabs(t0j-1-t0j-2)>=e;l+)for(k=1,s=0

17、;k<=pow(2,l-1);k+)s+=f(a+(2*k-1)*(b-a)/pow(2,l);/判断前后两次所得的T(0)否符合要求，如果符合精度要求那么停止循环tl1=(tl-11+(b-a)*s/pow(2,l-1)/2;for(i=l-1,j=2;i>=0;i-,j+) tij=(pow(4,j-1)*ti+1j-1-tij-1)/(pow(4,j-1)-1);if(t01<e)prin tf("t=%0.6fn",t01);elseprintf(" 用 Romberg算法计 算函数 所得近 似结果 nf(x)=%0.6fn",

18、t0j-1);3.3运行结果:'D WC+ -EXERCISE3XDtbug3.eRThe yeisLiil't is "440,5301 Piess bfi¥to continue3.4 MATLAB上机程序function T,n=mromb(f,a,b,eps)if narginv4,eps=1e-6; endh=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b);n=1;J=O;err=1;while (err>eps)J=J+1;h=h/2;S=0;for i=1: nx=a+h*(2*i-1);S=S+feval

19、(f,x);endR(J+1,1)=R(J,1)/2+h*S;for k=1:JR(J+1,k+1)=(4Ak*R(J+1,k)-R(J,k)/(4Ak-1);enderr=abs(R(J+1,J+1)-R(J+1,J);n=2*n;endR；T=R(J+1,J+1);format lo ngf=(x)(3.Ax)*(x.A1.4)*(5*x+7)*si n(x*x);T, n=mromb(f,1,3,1.e-5)3.5问题讨论：1. Romberge算法的优点是：把积分化为代数运算，而实际上只需求Ti(i),以后用递推可得.算法简单且收敛速度快，一般4或5次即能到达要求。2. Romberg

20、e 算法的缺点是 :对函数的光滑性要求较高， 计算新分点 的值时，这些数值的个数成倍增加。3. 该程序较为复杂，涉及函数定义，有循环，而且循环中又有判 断，编写时需要注意该判断条件是处于循环中， 当到达要求时跳出循 环，终止运算。4. 函数的定义可放在主函数前也可在主程序后面。本程序采用的 后置方式。4.用定步长四阶Runge-Kutta求解dy1 /dt 1dy2 /dt y3dy3 /dt 1000 1000y2 100y3y1(0) 0y2(0) 0y3(0) 0h=0.0005,打印 yi(0.025) , yi(0.045) , yi(0.085) , yi(0.1),(i=1， 2

21、， 3)4.1 理论依据：Ru nge_Kutta采用高阶单步法，这里不是先按Taylor公式展开，而是先写成tn处附近的值的线性组合(有待定常数)再按 Taylor公式展开，然后确定待定常数，这就是 Runge-Kutta 法的思想方法。此题采用四阶古典的Run ge-Kutta公式：Yn 1 Yn K1 3K2 3K3 K4/8K1hF ( xn ,Yn )K2hF(xnh/3,YnhK1/3)K3hF(xn2h/3,YnhK1 /3hK2)K4hF(xnh,Yn hK1 hK2hK3)4.2 C 语言程序代码：#include<stdio.h>void fun(double

22、x4,double y4,double h) y1=1*h;y2=x3*h;y3=(1000-1000*x2-100*x2-100*x3)*h;/微分方程向量函数 void main() double Y54,K54,m,z4,e=0.0005;double y5=0,0.025,0.045,0.085,0.1;int i,j,k;for(i=1;i<=3;i+)Y1i=0;for(i=1;i<=4;i+)for(j=1;jv=3;j+)Kij=O;for(k=1;k<=5;k+)for(m=yk-1;m<=yk;m=m+e)for(i=1;i<=3;i+)zi=

23、Yki;fun (z,K1,e);for(i=1;i<=3;i+)zi=Yki+e*K2i/2;依此求 K1,K2K3 的值fun (z,K2,e); for(i=1;i<=3;i+)zi=Yki+e*K2i/2;fun( z,K3,e);for(i=1;i<=3;i+)zi=Yki+e*K3i;fun (z,K4,e);/求for(i=1;i<=3;i+)Yki=Yki+(K1i+2*K2i+2*K3i+K4i)/6;YiN+1的值if(k!=5)for(i=1;i<=3;i+)Yk+1i=Yki;printf("计算结果:n");for(i

24、=1;i<5;i+)for(j=1;j<=3;j+)pri ntf("y%d%4.3f=%-10.8f,",j,yi,Yij);if(j=3)prin tf("n");prin tf("n");4.3运行结果：025000 yl<0.a45>=0.045000 pl<0-085>=0-a8590S yl<S.10G>=S丄司日y2<0.025>-6.149501 y2<0.045>=0.310966 Lf2<B-B85> =0-559342y3<

25、0.025>-S_33575& y3<0.045>=7-536549 p3<a-685=4-946488Ppeis Artii k电y to continue4.4 问题讨论：1.定步长四阶 Runge-kutta 方法是一种高阶单步法法稳定，精度较 高，误差小且程序相对简单，存储量少。不必求出起始点的函数值， 可根据精度的要求修改步长，不会由于起始点的误差造成病态。2.本程序可以通过修改主程序所调用的函数中的表达式来实现 对其它函数的任意初值条件求微分计算。3. 程序中运用了大量的 for 循环语句，因为该公式中涉及大量的求 和，且有不同的函数和对不同的数值求

26、值，编程稍显繁琐。所以编写 过程中一定要注意各循环的次数，以免出错。5.12.384122.115237-1.0610741.112336- 0.1135840.7187191.7423823.0678132.0317432.11523719.141823- 3.125432- 1.0123452.1897361.563849- 0.7841651.1123483.123124-1.061074- 3.12543215.5679143.1238482.0314541.836742- 1.0567810.336993- 1.0101031.112336- 1.0123453.12384827.1

27、084374.101011- 3.7418562.101023- 0.71828- 0.037585A- 0.1135842.1897362.0314544.10101119.8979180.431637- 3.1112232.1213141.7843170.7187191.5638491.836742- 3.7418560.4316379.789365- 0.103458- 1.1034560.2384171.742382- 0.784165-1.0567812.101023- 3.111223- 0.10345814.7138465 3.123789- 2.2134743.0678131.

28、1123480.336993- 0.718282.121314-1.1034563.12378930.7193344.446782- 2.0317433.123124-1.010103- 0.0375851.7843170.238417- 2.2134744.44678240.00001b(2.187436933.992318- 25.1734170.846716951.784317- 86.6123431.11012304.719345 - 5.6784392 )用列主元消去法求解 Ax=b。5.1 理论依据：列主元素消元法是在应用Gauss消元法的根底上，凭借长期经验 积累提出的， 是线性

29、方程组一般解法， 目的是为防止在消元计算中使 误差的扩大， 甚至严重损失了有效数字使数据失真， 而在每次初等变 换前对矩阵作恰当的调整，以提高Gauss消元法的数字稳定性，进而 提高计算所得数据的精确度。 即在每主列中取绝对值最大的元素作主 元，再做对应的行交换然后消元求解的方法。具体做法如下：将方阵A和向量b写成C=( A，b)。将C的第1列中第1行的 元素与其下面的此列的元素逐一进行比较，找到最大的元素Cji，将第 j 行的元素与第 1 行的元素进行交换，然后通过行变换，将第 1 列中 第2到第n个元素都消成0。将变换后的矩阵c的第二列中第二行 的元素与其下面的此列的元素逐一进行比较，找到

30、最大的元素ck(12) ，将第k行的元素与第2行的元素进行交换，然后通过行变换，将第 2 列中第3到第n个元素都消成0。以此方法将矩阵的左下局部全都消成 0 后再求解。最终形式如下：a(n) a11 0LOLa(n)a1nMg1M(A,b )MOOMM0L0a(n) nngn5.2 C 语言程序代码(1)比较该列的元素的绝对值的大小，将绝对值最大的元素通过行 变换使其位于主对角线上；(2)进行高斯消去法变换，把系数矩阵化成上三角形，然后回代求#include "math.h"#include "stdio.h"void Householder(doubl

31、e A99);void expunction(double A99,double b9,double x9);void main()double A99=12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.8367

32、42,-1.056781,0.336993,-1.010103,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0. 037585,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,1.742382,-0.78

33、4165,-1.056781,2.101023,-3.111223,-0.103458,14.713847,3.123789,- 2.213474,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.4 46782,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,4 0.00001;double b9= 2.1874369,33.992318,-25.173417,0.84671695,1.7

34、84317,-86.612343,1.1101230,4.71 9345,-5.6784392;double x9=0.0;int i,j;Householder(A);printf("n The Results of X are:n");expunction(A,b,x);for(i=1;i<10;i+)printf("X%1d=%fn",i,xi-1);void Householder(double A99)double q9,u9,y9,s,a,kr;int i,j,k;for(i=0;i<7;i+)s=0;for(j=i+1;j<

35、;9;j+) s+=Aji*Aji;s=sqrt(s);a=s*s+fabs(Ai+1i)*s;for(j=0;j<9;j+)if(j<=i) uj=0;else if(j=i+1) uj=Aji+Aji/fabs(Aji)*s;else if(j>i+1) uj=Aji;for(k=0;k<9;k+)yk=0;for(j=0;j<9;j+)yk+=Akj*uj;yk/=a;kr=0;for(k=0;k<9;k+)kr+=yk*uk;kr/=2*a;for(k=0;k<9;k+)qk=yk-kr*uk;for(k=0;k<9;k+)for(j=0

36、;j<9;j+)Akj-=uk*qj+uj*qk;void expunction(double A99,double b9,double x9)int i,j,k;double B910;double z3;double t1=0,t2=0,t3=0;for(i=0;i<8;i+)if(Ai+1i>Aii)for(j=i,k=0;j<i+3;j+,k+)zk=Aij;Aij=Ai+1j;Ai+1j=zk;t1=bi;bi=bi+1;bi+1=t1;t2=Ai+1i;for(j=i;j<i+3;j+)Ai+1j=Ai+1j-Aij*t2/Aii;bi+1=bi+1-

37、bi*t2/Aii; x8=b8/A88;for(i=7;i>=0;i-)for(j=i+1;j<9;j+)t3=t3+Aij*xj;xi=(bi-t3)/Aii;t3=0;5.3运行结果*' D=VC + + EXERC5 + ewe"The JEes u 1of K Ai*e =1-1-G757992=2 27S7443=-2 合F95：L占"1=2 2930775-2-112&346=-6-4238377-1-3579238=0-6342449=-0-587266i*ess a.ny ko_y to continue5.4 MATLAB上机程序un cti on x=mgauss2(A,b,flag)if narginv 3,flag=0;e ndn=len gth(b);for k=1:( n-1)ap,p=max(abs(A(k: n, k);p=p+k-1;if p>kA(k p,:)=A(p k,:);b(k p,:)=b(p k,:);endm=A(k+1: n,k)/A(k,k);A(k+1: n,k+1: n)=A(k