C++数值分析列主元高斯消去法列主元LU分解法Jacobi迭

1、数值分析上机实验报告华南理工大学08级计算机科学与技术2班1,三次样条插值（初值条件1）：P52.9、给定函数的函数表和边界条件，求三次样条插值函数，并求的近似值。函数表7576777879802.7682.8332,9032.9793.0623.153源代码：yangtiao.cpp#include#includevoidmain()(intchoice=0;intn=2;doublexx,*x,*y,*a,*b,*a1,*b1,*h,*m;cout请输入插值节点个数n：n;x=newdoublen;y=newdoublen;a=newdoublen;b=newdoublen;a1=newd

2、oublen;b1=newdoublen;h=newdoublen-1;m=newdoublen+1;cout请输入n个插值的节点(xi,yi)：endl;for(inti=0;ixiyi;)for(intj=0;jn-1;j+)(hj=xj+1-xj;)cout请输入彳f估点xx：xx;cout请选择边界条件：choice;switch(choice)(case 1:(doubletemp1,temp2;a0=0;an-1=1;cout请输入边界条件的两个一阶微商值s(x1)与s(xn):temp1temp2;b0=2*temp1;bn-1=2*temp2;break;case 2:(a0=

3、1;an-1=0;b0=3/h0*(y1-y0);bn-1=3/hn-2*(yn-1-yn-2);break;)for(intk=1;kn-1;k+)(ak=hk-1/(hk-1+hk);bk=3*(1-ak)/hk-1*(yk-yk-1)+ak/hk*(yk+1-yk);)a10=-a0/2;b10=b0/2;for(intl=1;l=0;j-)(mj=a1j*mj+1+b1j;)/判别xx所在区间并输出结果coutn插值结果为：;for(k=0;kn-1;k+)(if(xkxx)(doubleoutput=0;output=(1+2*(xx-xk)/(xk+1-xk)*pow(xx-xk+

4、1)/(xk-xk+1),2)*yk+(1+2*(xx-xk+1)/(xk-xk+1)*pow(xx-xk)/(xk+1-xk),2)*yk+1+(xx-xk)*pow(xx-xk+1)/(xk-xk+1),2)*mk+(xx-xk+1)*pow(xx-xk)/(xk+1-xk),2)*mk+1;coutoutput：0.250.50.30.54770-390.624S0.450.67380.530.728请输!入待估点：0.35请选择边界条件请输入边界条件的两个一阶微商值与&1j10.6868插值结果为：99160?Pressanykeytocontinue3 .自动选取步长梯形法：P977

5、、使用自动选取步长梯形法计算积分的近似值。(给定=0.01)源代码：SelfSelLength.cpp#include#includedoublefun(doublea)return2/(1+a*a);doubleSelfSelLength(doubleR_a,doubleR_b,doublee)doubleh=(R_b-R_a)/2;doubleR1=(fun(R_a)+fun(R_b)*h;intn=1;doubleR0;doubleS;doubleE;do/每当误差值不符合要求时，计算下一个result值R0=R1;S=0;for(intk=1;k3*e);returnR1;voidma

6、in()doublea,b,e;cout请依次输入待求积分函数的下界a、上界b及精度要求e:abe;cout自动选取步长梯形法可求得积分：SelfSelLength(a,b,e)endl;运行结果截图：4 .Romberg求积法：P978、使用Romberg求积法计算积分的近似值。（给定e=0.01,且取）源代码：Romberg.cpp#include#includestaticdoubleTri128128;doublefun(doublea)returnsqrt(a);doubleRomberg(doubleR_a,doubleR_b,doublee)Tri00=(R_b-R_a)/2*(

7、fun(R_a)+fun(R_b);intk=0;doubleE;do/每当误差值不符合要求时，计算下一行的Tri口值k+;doubletemp=0;/计算T0k的数值for(inti=1;i=pow(2,k-1);i+)temp=temp+fun(R_a+(2*i-1)*(R_b-R_a)/pow(2,k);Tri0k=0.5*(Tri0k-1+(R_b-R_a)/pow(2,k-1)*temp);for(intm=1;me);returnTrik0;voidmain()doublea,b,e;cout请依次输入待求积分函数的下界a、上界b及精度要求e:abe;cout按Romberg求积法

8、可求得积分:Romberg(a,b,e)endl;运行结果截图：请依次输入特求积分函数的下界击上界h及精度要求八L90.01按刖mber求积法可求得积分;173328?i*ecsanykeytocontinue5 .列主元高斯消去法:测试矩阵为：源代码：Guess_Elimination.cpp/列主元高斯消去法#include#include#defineN4/矩阵的维数，可按需更改系数矩阵staticdoubleANN=2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/staticdoubleBN=1,0,1,0;/右端项staticdoubleXN;inti,

9、j,k;/计数器voidmain()for(k=0;kN-1;k+)/选取最大主元intindex=k;for(i=k;iN;i+)if(fabs(Aindexk)fabs(Aik)index=i;/交换行doubletemp;for(i=k;iN;i+)temp=Aindexi;Aindexi=Aki;Aki=temp;temp=Bindex;Bindex=Bk;Bk=temp;for(i=k+1;iN;i+)doubleT=Aik/Akk;Bi=Bi-T*Bk;for(j=k+1;j=0;i-)doubleTemp=0;for(intj=i+1;jN;j+)Temp=Temp+Aij*Xj

10、;Xi=(Bi-Temp)/Aii;cout线性方程组的解(X1,X2,X3.Xn)为：endl;for(i=0;iN;i+)coutXi;运行结果截图：XlfX2,X3XnL.21.41.60.8Pressanykeytocontinue6 .列主兀LU分解法:测试矩阵为:源代码：LU_Decomposition.cpp#include#include#defineN4/矩阵维数，可自定义/staticdoubleANN;staticdoubleBN;staticdoubleYN;staticdoubleXN;staticdoubleSN;系数矩阵右端项中间项输出选取列主元的比较器inti,

11、j,k;/计数器voidmain()(ai1,ai2,ai3ain,yi):endl;cout请输入线性方程组for(i=0;iN;i+)for(intj=0;jAij;cinBi;for(k=0;kN;k+)/选列主元intindex=k;for(i=k;iN;i+)(doubletemp=0;for(intm=0;mk;m+)(temp=temp+Aim*Amk;)Si=Aik-temp;if(SindexSi)(index=i;)/交换行doubletemp;for(i=k;iN;i+)(temp=Aindexi;Aindexi=Aki;Aki=temp;)temp=Bindex;Bin

12、dex=Bk;Bk=temp;/构造L、U矩阵for(j=k;jN;j+)(doubletemp=0;for(intm=0;mk;m+)(temp=temp+Akm*Amj;)Akj=Akj-temp;/先构造U一行的向量)for(i=k+1;iN;i+)(doubletemp=0;for(intm=0;mk;m+)(temp=temp+Aim*Amk;)Aik=(Aik-temp)/Akk;/再构造L一列的向量)/求解LY=BY0=B0;for(i=1;iN;i+)(doubletemp=0;for(intj=0;j=0;i-)(doubletemp=0;for(intj=i+1;jN;j+

13、)(temp=temp+Aij*Xj;)Xi=(Yi-temp)/Aii;)/打印Xcout线性方程组的解(X1,X2,X3.Xn)为：endl;for(i=0;iN;i+)(coutXi;)运行结果截图:线性方程组的解(XI,X2,X3Xn)为:口.21.41.60.8Pfcssanykeytocontinue7.简单迭代法(Jacobi迭代)：测试矩阵为：=取精度为esp=0.001,最大迭代次数M=100源代码：Jacobi.cpp#include#include#defineN4/矩阵的维数，可按需更改系数矩阵staticdoubleANN=2,-1,0,0,-1,2,-1,0,0,-

14、1,2,-1,0,0,-1,2;/staticdouble丫N;staticdoubleXN;staticdoubleGN;staticdoubleBN=1,0,1,0;/右端项/输出比较项staticdoubleYN;/输出项/X=BX+G的G矩阵inti,j,k;/计数器doubleeps=0.001;intM=100;booldistance()/求两输出项的差的范数是否满足精度要求doubletemp=0;for(i=0;ieps)returnfalse;elsereturntrue;/满足精度要求则结束程序voidmain()/形成迭代矩阵B,存放到A中for(i=0;iN;i+)i

15、f(fabs(Aii)eps)cout打印失败endl;return;doubleT=Aii;for(j=0;jN;j+)Aij=-Aij/T;Aii=0;Gi=Bi/T;intcounter=0;while(counterM)(/迭代for(i=0;iN;i+)(doubletemp=0;for(j=0;jN;j+)(temp=temp+Aij*Yj;Xi=Gi+temp;if(distance()=true)break;else(/交换X,Y向量；for(i=0;iN;i+)(Yi=Xi;counter+;/打印Xcout迭代次数为：counter次。该线性方程组的解(X1,X2,X3Xn

16、)为：endl;for(i=0;iN;i+)(coutXi;运行结果截图:8.Seidel迭代:测试矩阵为：=取精度为esp=0.001,最大迭代次数M=100源代码：Seidel.cpp#include#include#defineN4/矩阵的维数，可按需更改staticdoubleANN=2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/系数矩阵staticdoubleBN=1,0,1,0;/staticdoubleYN;staticdoubleXN;staticdoubleGN;inti,j,k;/doubleeps=0.001;intM=100;/X=BX+

17、G计数器右端项输出比较项输出项的G矩阵booldistance()/求两输出项的差的范数是否满足精度要求doubletemp=0;for(i=0;ieps)returnfalse;elsereturntrue;/满足精度要求则结束程序voidmain()/形成迭代矩阵B,存放到A中for(i=0;iN;i+)if(fabs(Aii)eps)cout打印失败endl;return;doubleT=Aii;for(j=0;jN;j+)Aij=-Aij/T;Aii=0;Gi=Bi/T;intcounter=0;while(counterM)/迭代for(i=0;iN;i+)doubletemp=0;

18、for(j=0;jN;j+)(temp=temp+Aij*Xj;)Xi=Gi+temp;)if(distance()=true)break;else(/交换X,Y向量;for(i=0;iN;i+)(Yi=Xi;)counter+;)cout迭代次数为：counter次。该线性方程组的解(X1,X2,X3Xn)为：endl;for(i=0;iN;i+)(coutXi)3idelDebugSeic迭代江数为：次.该线性方程组的解1；X31.199561.399421.599530.799766PessanyKeytocontinue9.松弛法（SOR迭代）测试矩阵为：取松弛系数w=1.46=取精度为esp=0.001,最大迭代次数M=100源代码：SOR.cpp#include#include系数矩阵#defineN4/矩阵的维数，可按需更改staticdoubleANN=2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2;/staticdoubleBN=1,0,1,0