数值分析上机报告.doc

_第一题：1、已知A与b（1）用Househloser变换，把A化为三对角阵（并打印B）。（2）用超松弛法求解Bx=b（取松弛因子=1.4，x（0）=0，迭代9次）。（3）用列主元素消去法求解Bx=b。一、分析如下：（3）用列主元素消去法求解Bx=b。将方阵A和向量b写成C=（A b）。将C的第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素，将第j行的元素与第1行的元素进行交换，然后通过行变换，将第1列中第2到第n个元素都消成0。将变换后的矩阵的第二列中第二行的元素与其下面的此列的元素逐一进行比较，找到最大的元素，将第k行的元素与第2行的元素进行交换，然后通过行变换，将第2列中第3到第n个元素都消成0。以此方法将矩阵的左下部分全都消成0。最终形式如下：（A b）二、程序：#include math.h#include stdio.h /* 标准的基本库函数 头文件 */#define ge 8void main()int sign(double x); /*该函数为sign函数*/double t99= /* 数组赋初值 */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,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.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-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 p,v,s,w;int i,j,h,m,g;double u9,x19,y9,q9,b1910,x9;double d9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392;/*-开始迭代-*/for(j=0;j7;+j) /*Household 变换 */s=0.0;for(i=j+1;i0)?(s*s+s*tj+1j):(s*s-s*tj+1j);for(g=0;g9;+g)if (g=j)ug=0;else if (g=j+1)ug=tj+1j+s*sign(tj+1j);else ug=tgj;for(m=0;m9;+m)ym=0;for(h=0;h9;+h)ym=ym+tmh*uh;ym=ym/v;p=0;for(i=0;i9;+i)p=p+ui*yi;p=0.5*p/v;for(i=0;i9;+i)qi=yi-p*ui;for(h=0;h9;+h)for(m=0;m9;+m)tmh=tmh-(qm*uh+um*qh);printf(Household:h);for(i=0;i9;+i)for(j=0;j9;+j)if (j%9=0)printf(n);printf(%-9.5f,tij);printf(n);w=1.4; /*超松弛法*/for(i=0;i9;i+) /*开始循环 */x1i=0;for(i=0;i9;i+)for(j=0;j9;j+)if(i=j)b1ij=0;else b1ij=-tij/tii;for(i=0;i9;i+)b1i9=di/tii;for(h=0;h9;h+)for(i=0;i9;i+)s=0;for(j=0;j9;j+)s=s+b1ij*x1j;s=s+b1i9;x1i=x1i*(1-w)+w*s;for(i=0;i9;i+)if (i=5)printf(n);printf(x%d=%-10.6f,i,x1i); /*输出结果*/printf(n);u0=t00; /*以下是消去法*/y0=d0;for(i=1;i=0;i-)xi=(yi-tii+1*xi+1)/ui;for(i=0;i=(1e-6)?1:-1); return(z); /* 主函数返回z 告诉系统 程序结束 */ 三、运行结果：四问题讨论：题中矩阵B是对称正定阵，且是三对角的，所以选择合适的松驰因子，收敛速度是很快的。已经求出x的第m解，第m-1基础上，经过重新组合得新的序列，而在此新序列收敛速度加快。第三题：3、已知函数值如下表：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)的近似值。一、分析如下：这里 ，所以只要求出，就能得出插值函数S（x）。求的方法为：这里最终归结为求解一个三对角阵的解。用追赶法解三对角阵的方法如下： ， 综上可得求解方程Ax=d的算法：二、程序如下：#include #include main()double th1010;int i;double d10;double y10=0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851,xm10=1,2,3,4,5,6,7,8,9,10;double h=1;a;double x10;double b10,l10;double c10,c110,c210;float f10,f110;for(i=0;i10;i+)thii=2;th01=1;for(i=1;i9;i+)thii+1=0.5;thii-1=0.5;th98=1; *输入矩阵值*d0=6*(y1-y0)/h-1)/h;for(i=1;i9;i+)di=6*(yi+1-yi)/h-(yi-yi-1)/h)/(2*h);d9=6*(0.1-(y9-y8)/h)/h; *计算d值*b0=th00;c0=d0;for(i=1;i=0;i-)xi=(ci-thii+1*xi+1)/bi; *用追赶法求解三样条插值法中的M*printf(please input:); *输入x值* scanf(%lf,&a);for(i=1;i10;i+) c1i=(xi-1*pow(h,2)/6+yi-xi*pow(h,2)/6-yi-1)/h; *求出c1值*c2i=yi-xi*pow(h,2)/6; *求出c1值* fi=xi-1*pow(xmi-a),3)/6+xi*pow(a-xmi-1),3)/6+c1i*(a-xmi)+c2i; *f(x)表达式* f1i=xi*pow(a-xmi-1),2)/2+c1i-xi-1*pow(xmi-a),2)/2; *f(x)表达式* printf(n); *以下为判断求输入的自变量x的表达式* if(0=a&a=1) printf(f=%8.7f,f1=%8.7f,f1,f11); else if(1=a&a=2) printf(f=%8.7f,f1=%8.7f,f1,f11); else if(2=a&a=3) printf(f=%8.7f,f1=%8.7f,f2,f12); else if(3=a&a=4) printf(f=%8.7f,f1=%8.7f,f3,f13); else if(4=a&a=5) printf(f=%8.7f,f1=%8.7f,f4,f14); else if(5=a&a=6) printf(f=%8.7f,f1=%8.7f,f5,f15); else if(6=a&a=7) printf(f=%8.7f,f1=%8.7f,f6,f16); else if(7=a&a=8) printf(f=%8.7f,f1=%8.7f,f7,f17); else if(8=a&a=9) printf(f=%8.7f,f1=%8.7f,f8,f18); else printf(f=%8.7f,f1=%8.7f,f9,f19);三、计算结果如下： f(4.563)=1.517932 f(4.563)=0.249350 四、问题讨论：此题最终能实现对任意一组数据的三次样条插值计算得出区间内任意一个点的函数值及其导数值。第四题：4、用Newton法求方程 在（0.1 , 1.9）中的近似根（初始近似值取为区间端点，迭伐6次或误差小于0.000 01）。一、 分析如下：对满足下列条件的函数：i. f(a)f(b)0；ii. f”(x)在区间a,b上不变号；iii. f(x)0；iv. |f(c)|/(b-a)|f(c)|，其中c 是a,b中使min(|f(a)|,|f(b)|)达到的一个。则对任意初始近似值，迭代方程为：二、程序如下：#include #includevoid main ( ) double x1=1.9,x2,f, f;dof=pow(x1,7)-28*pow(x1,4)+14;f=7*pow(x1,6)-28*4*pow(x1,3);x2=x1-f/f;=x1-x2;x1=x2;while (fabs()=pow(10,-6);printf(“%f”,x2); 三、计算结果如下：equation is : 1.00 x7 - 28.00 x4 + 14.00 = 0x = 0.8454966四、问题讨论：当x=1.9时运行结果在0.1，1.9之间，是要求的解；当x=0.1时运行结果不在0.1，1.9之间，不符合要求。这并不是说当x=0.1时所得解不是原方程的解，而是所得解不在所要求的范围内。程序对不同的方程求解，或有不同的误差要求，只需对do while 循环中量做适当的修改即可，程序没有考虑f(x)=0的情况，有可能产生死循环。第五题：5、用Romberg算法求（允许误差=0.000 01）。一、本题分析如下：Romberg算法的计算步骤如下：（1）先求出按梯形公式所得的积分值（2）把区间2等分，求出两个小梯形面积之和，记为，即这样由外推法可得,。（3）把区间再等分（即等分），得复化梯形公式，由与外推可得，如此，若已算出等分的复化梯形公式，则由Richardson外推法，构造新序列， m=1,2,l, k=1,2,l-m+1,最后求得。（4）或就停止计算，否则回到（3），计算，一般可用如下算法：其具体流程如下，并全部存入第一列 Romberg算法的一些说明：通常计算时，用固定L=N来计算，一般L=4或5即能达到要求。应用Romberg算法时，可以先给出所要求的误差，当就停机。Romberg算法的优点是：（1）把积分化为代数运算，而实际上只需求，以后用递推可得。（2）算法简单且收敛速度快。Romberg算法的缺点是：（1）对函数的光滑性要求较高。（2）计算新分点的值时，这些数值的个数成倍增加。二 、程序如下：#include stdio.h#includemath.hmain() int i,k,e; double ya,yb,a,b,tt,h,m,eps; double t1010; a=1; b=3; h=b-a; m=0; tt=0; eps=0.001; ya=rem(a); yb=rem(b); t01=(ya+yb)*h/2; for(k=1;k=9;k+) for(i=1;i=pow(2,k-1);i+) m=a+(i-0.5)*h/pow(2,k-1); tt=tt+rem(m) ; m=0; tk1=tk-11+h*tt/pow(2,k-1);tk1=tk1/2; tt=0; for(k=1;k=8;k+) for(i=0;i=9-k;i+) tik+1=(pow(4,k)*ti+1k-tik)/(pow(4,k)-1); if(abs(t0k+1-t0k)eps) break; printf(Romberg算法求得的积分结果为: %e n,t0k+1); rem(x)double x;double y; y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return(y);三、计算结果如下：The result: y=440.535928四、问题讨论：本程序可通过一个一维数组将进行存储，但考虑到对程序的可读性，故将存储在一个二维数组里。第六题：6、用定步长四阶Runge-Kutta法求解h=0.000 5，打印。一、分析如下：二、 程序如下：#include#include#define N 3float yN=0,0,0,kNN+1;float ydN+1N,aN+1=0.025,0.045,0.085,0.10;/*以上是预处理和定义全局变量*/main()int i,l,j=0,n=1; float h=0.0005;dofor(l=0;lN;l+)for(i=0;i=N;i+)/*计算每一个k值*/kli=0;k00=h;k10=h*y2;k20=h*(1000-1000*y1-100*y2);k01=h;k11=h*(y2+(1/3.0)*k00);k21=h*(1000-1000*(y1+(1/3.0)*k10)-100*(y2+(1/3.0)*k10);k02=h;k12=h*(y2+(1/3.0)*k10+k11);k22=h*(1000-1000*(y1+(1/3.0)*k20+k21)-100*(y2+(1/3.0)*k20+k21);k03=h;k13=h*(y2+k10-k11+k12);k23=h*(1000-1000*(y1+k20-k21+k22)-100*(y2+k20-k21+k22);for(i=0;iN;i+) yi=yi+(ki0+3*ki1+3*ki2+ki3)/8;