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

1、数值分析上机实验报告 ?数值分析?上机实验报告1.用Newton法求方程X7-X4+14=0在0.1，1.9中的近似根初始近似值取为区间端点，迭代6次或误差小于0.00001。1. 1 理论依据：设函数在有限区间a，b上二阶导数存在，且满足条件令故以1.9为起点如此一次一次的迭代，逼近x的真实根。当前后两个的差<=时，就认为求出了近似的根。本程序用Newton法求代数方程(最高次数不大于10)在a,b区间的根。1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main()double x2,f,f1; double x1=

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

3、/feval(df,x0); end e=abs(x1-x0); x0=x1; if e<=eps&&abs(feval(f,x1)<=eps d=1;break endendif d=1 y=x1;elseif d=0 y='迭代M次失败'else y= '奇异'endfunction y=df(x)y=7*x6-28*4*x3;Endfunction y=f(x)y=x7-28*x4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newto

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

5、Newton迭代的初值，结果发现0.1不满足条件，而1.9满足，能作为初值。另外，该程序简单，只有一个循环，且为顺序结构，故采用do-while循环。当然也可以选择for和while循环。2.函数值如下表：x12345f(x)00.693147181.09861231.38629441.6094378x678910f(x)1.79175951.94591012.0794452.19722462.3025851f(x)f(1)=1f(10)=0.1试用三次样条插值求f(4.563)及f(4.563)的近似值。2.1 理论依据这里 ，所以只要求出，就能得出插值函数Sx。求的方法为：这里最终归结为求

6、解一个三对角阵的解。用追赶法解三对角阵的方法如下： ， 综上可得求解方程Ax=d的算法：2.2 C语言程序代码：#include<stdio.h>#include<math.h>void main()int i,j,m,n,k,p;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.079

7、445,2.1972246,2.3025851; printf("请依次输入xi:n"); for(i=0;i<=9;i+) scanf("%lf",&ei); /求h矩阵 for(n=0;n<=8;n+) hn=en+1-en; 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

8、); 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;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"); p

9、rintf("根据追赶法解三对角矩阵得：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; x9=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上

10、的表达式;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=-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. 三次样条插值效果比Lag

11、range插值好，没有Runge现象，光滑性较好。2. 此题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。3. 编程过程中由于定义的数组比拟多，需要仔细弄清楚各数组所代表的参数，要注意各下标代表的含义，特别是在用追赶法计算的过程中。3.用Romberg算法求.3.1 理论依据：Romberg算法的计算步骤如下：1先求出按梯形公式所得的积分值2把区间2等分，求出两个小梯形面积之和，记为，即这样由外推法可得,。3把区间再等分即等分，得复化梯形公式，由与外推可得，如此，假设已算出等分的复化梯形公式，那么由Richardson外推法，构造新序列， m=1,2,l, k=1,2,l-m+1,

12、最后求得。4或<就停止计算，否那么回到3，计算，一般可用如下算法：其具体流程如下，并全部存入第一列 通常计算时，用固定l=N来计算，一般l=4或5即能到达要求。3.2 C语言程序代码：#include<math.h>#include<stdio.h>double f(double x) /计算f(x)的值 double z; z=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return(z);main() double t2020,s,e=0.00001,a=1,b=3; int i,j,l,k; t01=(b-a)*(f(b)+f

13、(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;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

14、+1j-1-tij-1)/(pow(4,j-1)-1); if(t01<e) printf("t=%0.6fn",t01); else printf("用Romberg算法计算函数所得近似结果为：nf(x)=%0.6fn",t0j-1);3.3 运行结果：3.4 MATLAB上机程序function T,n=mromb(f,a,b,eps)if nargin<4,eps=1e-6;endh=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b);n=1;J=0;err=1;while (err>eps) J=J+

15、1;h=h/2;S=0; for i=1:n x=a+h*(2*i-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; for k=1:J R(J+1,k+1)=(4k*R(J+1,k)-R(J,k)/(4k-1); end err=abs(R(J+1,J+1)-R(J+1,J); n=2*n;endR;T=R(J+1,J+1);format longf=(x)(3.x)*(x.1.4)*(5*x+7)*sin(x*x);T,n=mromb(f,1,3,1.e-5)3.5 问题讨论：1.Romberge算法的优点是:把积分化为代数运算,而实际上只需求T

16、1(i),以后用递推可得.算法简单且收敛速度快,一般4或5次即能到达要求。2.Romberge算法的缺点是:对函数的光滑性要求较高，计算新分点的值时，这些数值的个数成倍增加。3.该程序较为复杂，涉及函数定义，有循环，而且循环中又有判断，编写时需要注意该判断条件是处于循环中，当到达要求时跳出循环，终止运算。4.函数的定义可放在主函数前也可在主程序后面。本程序采用的后置方式。4. 用定步长四阶Runge-Kutta求解h=0.0005，打印yi(0.025) ， yi(0.045) ， yi(0.085) ， yi(0.1) ，(i=1，2，3)4.1 理论依据：Runge_Kutta采用高阶单步

17、法，这里不是先按Taylor公式展开，而是先写成处附近的值的线性组合有待定常数再按Taylor公式展开，然后确定待定常数，这就是Runge-Kutta法的思想方法。此题采用四阶古典的Runge-Kutta公式：4.2 C语言程序代码：#include<stdio.h>void fun(double 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

18、.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;j<=3;j+) Kij=0; for(k=1;k<=5;k+) for(m=yk-1;m<=yk;m=m+e) for(i=1;i<=3;i+) zi=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);

19、 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=1;i<5;i+) for(j=1;j<=3;j+) printf("y%d%4.3f=%-10.8f,",j,yi,Yij); if(j=3) printf("n");pri

20、ntf("n");4.3 运行结果：4.4 问题讨论：1.定步长四阶Runge-kutta方法是一种高阶单步法法稳定，精度较高，误差小且程序相对简单，存储量少。不必求出起始点的函数值，可根据精度的要求修改步长，不会由于起始点的误差造成病态。2.本程序可以通过修改主程序所调用的函数中的表达式来实现对其它函数的任意初值条件求微分计算。3.程序中运用了大量的for循环语句，因为该公式中涉及大量的求和，且有不同的函数和对不同的数值求值，编程稍显繁琐。所以编写过程中一定要注意各循环的次数，以免出错。5. 用列主元消去法求解Ax=b。5.1 理论依据：列主元素消元法是在应用Gauss消

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

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

23、in()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.836742,-1.056781,0.336993,-1.010103, 1.112336,-1.012345,3.123848,2

24、7.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.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713847,3.12378

25、9,-2.213474,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001; double b9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392; double x9=0.0; i

26、nt 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<9;j+) s+=Aji*Aji; s=sqrt(s); a=s*s+fabs(Ai+1i)*s; fo

27、r(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;j<9;j+) Akj-=uk*qj+uj*qk; void

28、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-bi*t2/Aii; x8=b8/A88;

29、 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 运行结果5.4 MATLAB上机程序unction x=mgauss2(A,b,flag)if nargin<3,flag=0;endn=length(b);for k=1:(n-1) ap,p=max(abs(A(k:n,k); p=p+k-1; if p>k A(k p,:)=A(p k,:); b(k p,:)=b(p k,:); end m=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1); if flag=0,Ab=A,b,endendx=zeros(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endformat longA=12.38412,2.115237,-1.061074,1.112336