数值分析上机作业.doc

数值分析上机实习题- 13 -*第一题 程序*一 程序要求:已知A=12.3841202.115237-1.0610741.112336-0.1135840.7187191.7423823.067813-2.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.1084374.101011-3.7418562.101023-0.718280-0.037585-0.1135842.1897362.0314544.10101119.8979190.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.7138463.123789-2.2134743.0678131.1123480.336993-0.7182802.121314-1.1034563.12378930.7193344.446782-2.0317433.123124-1.010103-0.0375851.7843170.238417-2.2134744.44678240.00001b=(2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230 4.719345,-5.6784392);1.用Householder变换，把A化为三对角阵B(并打印B).2.用超松弛法求解Bx=b(取松弛因子=1.4,X(0)=0,迭代9次)3.用列主元素消去法求解Bx=b。二解题算法1.HOUSEHOLDER算法步骤： 令A0=A, aij(1)=aij，已知Ar-1 即Ar-1=(aij(r) Sr=(air(r )2)1/2 r=Sr2+|a(r)r+1,r|Srur=0,，0,a(r)r+1,r+Sign(a(r)r+1,r)Sr,a(r)r+2,r,anr(r)T yr=Ar-1ur/r kr= urTyr/2r qr=yr-krur Ar=Ar-1-(qrurT+urqrT)， r=(1,2,，n-2)2. SOR解题算法 其基本思想是在高斯方法已求出x(m),x(m-1)的基础上，经过重新组合的新序列，而此新序列收敛速度加快。其算式是：xi(m)=(1-)xi(m-1)+(bijxi(m)+ xj(m-1)+gi)其中是超松弛因子，当1时，可以加快收敛速度。3. 列主元素消去法程序算法 对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。然后进行行变换，把矩阵化为上三角阵，再进行回代，求出方程的解。三 主要程序#include #define N 9main()float eN+1,fN+1,x1N+1,A0NN,BNN,yN,uN,wN,x0N;float xN,qN,q1N,u1N,ma,a1N,s,s2,k,b1N;static float bN=2.1874369,33.992318,-25.173417,0.8471695,1.784317,-86.612343,1.1101230,4.719345,-5.56784392;int i,j,r,sign;float ANN=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.784163,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.742384,-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; clrscr(); for(i=0;iN;i+) for(j=0;jN;j+) A0ij=Aij; for(r=0;rN-2;r+) s=0; for(i=r+1;i0) sign=1; if(Ar+1r0) sign=-1; for(i=0;iN;i+) if(i=r) ui=0; else if(i=r+1) ui=Ar+1r+sign*s; else ui=Air; for(i=0;iN;i+) yi=0; for(j=0;jN;j+) yi=yi+Aij*uj; yi=yi/ma; k=0; for(i=0;iN;i+) k=k+ui*yi; k=k/(2*ma); for(i=0;iN;i+) qi=yi-k*ui; for(i=0;iN;i+) for(j=0;jN;j+) Aij=Aij-qi*uj-qj*ui; for(i=0;iN;i+) for(j=0;jN;j+) Bij=Aij; for(i=0;iN;i+) for(j=0;jN;j+) printf(%8.4f,Bij); printf(n); for(i=0;iN;i+) x0i=0; xi=0; for(i=0;iN;i+) x0=-(B01/B00)*x1+b0/B00; x0=1.4*x0-0.4*x00;for(j=1;jN-1;j+) xj=-(Bjj-1/Bjj)*xj-1-(Bjj+1/Bjj)*x0j+1+bj/Bjj; xj=1.4*xj-0.4*x0j; xN-1=-(BN-1N-2/BN-1N-1)*x0N-2+bN-1/BN-1N-1; xN-1=1.4*xN-1-0.4*x0N-1; for(j=0;jN;j+) x0j=xj; printf(n x:n); for(i=0;iN;i+) printf(x%d=%10.6fn,i,xi); for(i=0;iN-1;i+) b1i+1=bi; a1i+2=Bi+1i; ei+1=Bii+1;b1i+1=Bii; b1N=BN-1N-1; q10=0; u10=0; a11=0; fN=bN-1; b1N=BN-1N-1; for(i=1;i0;i-) x1i=q1i*x1i+1+u1i; printf(n x:n); for(i=1;iN+1;i+) printf(x%d=%10.6fn,i,x1i); getch(); 第一题程序结果：12.3841 -4.8931 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -4.8931 25.3984 6.4941 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 6.4941 20.6115 8.2439 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 8.2439 23.4228-13.8801 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000-13.8801 29.6983 4.5345 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 4.5345 16.0061 4.8814 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 4.8814 26.0133 -4.5037 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -4.5037 21.2541 4.5045 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 4.5045 14.5341 x: x0= 1.073156 x1= 2.273598 x2= -2.853854x3= 2.295630x4= 2.109620x5= -6.432651x6= 1.349386x7= 0.649445x8= -0.555408 x:x1= 0.000042x2= 0.000105x3= -0.000380x4= 0.000867x5= 0.001238x6= -0.005451x7= 0.016724x8= 0.090691x9= -0.411195*第三题 程序*/第三题：用三次样条插值法求插值及导数/一.程序要求已知十点函数值及起点和终点的导数值,要求用三次样条插值求f(4.563)及f(4.563)的近似值.x12345f(x)00.693147181.09861231.38629441.6094378x678910f(x)1.79175951.94591012.0794452.19722462.3025851f(x)f(x)=1f(x)=0.1二.程序算法对于均匀分划的插值函数取B样条函数为基函数，则任意样条函数可用线形组合表示，写成矩阵的形式： 用三对角阵的追赶法求解 追赶法步骤为以下递推式 y=(x-xj)/h 用上述得出的解代入公式，并根据3（y）的取值 得该题的解 f(x)=s(x) f(x)=s(x)且 s(x)= 且3(x)的表达式为以下所示 三次样条插值法求插值及导数#include /*本程序利用已知数据,采用样条插值公式,求解在某点的函数值及导数*/#includemain()float func1(float t1),func2(float t2); float b12=2,0,6*0.69314718,6*1.0986123,6*1.3862944, 6*1.6094378,6*1.7917595,6*1.9459101, 6*2.0794415,6*2.1972246,6*2.3025851,0.2; float c12=0,0,0,0,0,0,0,0,0,0,0,0; float x12=0,1,2,3,4,5,6,7,8,9,10,11; int i,j,k; int a1212=-1,0,1,0,0,0,0,0,0,0,0,0, 1,4,1,0,0,0,0,0,0,0,0,0, 0,1,4,1,0,0,0,0,0,0,0,0, 0,0,1,4,1,0,0,0,0,0,0,0, 0,0,0,1,4,1,0,0,0,0,0,0, 0,0,0,0,1,4,1,0,0,0,0,0, 0,0,0,0,0,1,4,1,0,0,0,0, 0,0,0,0,0,0,1,4,1,0,0,0, 0,0,0,0,0,0,0,1,4,1,0,0, 0,0,0,0,0,0,0,0,1,4,1,0, 0,0,0,0,0,0,0,0,0,1,4,1, 0,0,0,0,0,0,0,0,0,-1,0,1; float b1,b2; float s1=0,s2=0; for (k=0;k=500;k+) for (i=0;i=11;i+) b1=0; b2=0; for (j=0;j=i-1;j+) b1=b1+aij*cj; for (j=i+1;j=11;j+) b2=b2+aij*cj; ci=(bi-b1-b2)/aii; for (j=0;j=2) m=0; else if(fabs(t1)=2.0/3.0) n=0; else if(-0.5=fabs(t2)&fabs(t2)=0.5) n=-t2*t2+3.0/4.0; else n=0.5*t2*t2-(3.0/2.0)*fabs(t2)+9.0/8.0; return(n); 第三题程序结果： 用三次样条插值公式求某点的近似函数值及导数的解为： 在4.563处的插值为:s(4.563)=1.517932 在3.000处的导数为:s(3.00)=0.334969*第四题程序*/*第四题：用NEWTON法求方程的根*/用Newton 法求方程 x7-28*x4+14=0在（0.1, 1.9）中的近似根（初始值为区间端点，迭代6次或误差小于0.00001）用牛顿迭代法求7次方程的根，当N=1时 所以设其初始值为1.9，方程为： x7-28x4+14=0 牛顿迭代公式为 f(x1)是曲线f（x）在x1点处的切线斜率 因此 x2是上述切线与x轴的交点。求出x2后再找出f（x2）,既f(x)在x2处的切线此切线与轴交于x3。如此一次一次的迭代，逼近x的真实根。当前后两个求出 的差=时，就认为求出了近似的根#include /*本程序采用牛顿迭代法求方程的根*/main()int n; float A,B,x2,x1; n=1; x2=1.2; do /*迭代过程*/ x1=x2; A=x1*x1*x1*x1*x1*x1*x1-28*x1*x1*x1*x1+14; B=7*x1*x1*x1*x1*x1*x1-112*x1*x1*x1; x2=x1-A/B; n=n+1; while(fabs(x2-x1)=1e-5); printf(%dn,n); printf(%fn,x2);getch(); 第四题程序结果： NEWTON迭代法求方程的近似根为：0.845497*第五题程序*/*数值分析第五题：数值积分的Romberg算法*/用ROMBERG算法 当就停止运算。#include /*本程序采用Romberg迭代法求积分的值*/#includemain()float sc(float z); float t112,t29,t38,t47,t56; float t65,t74,t83,t92; float g10=0,0,0,0,0,0,0,0,0,0; int i,j,k; t10=36*sin(1)+27*pow(3,1.4)*22*sin(9); for(j=1;j=8;j+) for(i=1;i=pow(2,j-1);i+) gj-1=gj-1+sc(1+(2*i-1)*pow(2,1-j); t1j=0.5*(t1j-1+pow(2,2-j)*gj-1); for(k=1;k=8;k+) t2k-1=(4*t1k-t1k-1)/(4-1); for(k=1;k=7;k+) t3k-1=(4*4*t2k-t2k-1)/(4*4-1); for(k=1;k=6;k+) t4k-1=(pow(4,3)*t3k-t3k-1)/(pow(4,3)-1); for(k=1;k=5;k+) t5k-1=(pow(4,4)*t4k-t4k-1)/(pow(4,4)-1); for(k=1;k=4;k+) t6k-1=(pow(4,5)*t5k-t5k-1)/(pow(4,5)-1); for(k=1;k=3;k+) t7k-1=(pow(4,6)*t6k-t6k-1)/(pow(4,6)-1); for(k=1;k=2;k+) t8k-1=(pow(4,7)*t7k-t7k-1)/(pow(4,7)-1); t90=(pow(4,8)*t81-t80)/(pow(4,8)-1); printf(%fn,t90); getch(); float sc(float z) float n; n=pow(3,z)*pow(z,1.4)*(5*z+7)*si