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

1数值分析上机实验报告华南理工大学08级计算机科学与技术2班1.三次样条插值（初值条件1）：P52.9、给定函数的函数表和边界条件，求三次样条插值函数，并求的近似值。函数表7576777879802.7682.8332.9032.9793.0623.153源代码：yangtiao.cpp#includeiostream.h#includemath.hvoidmain(){intchoice=0;intn=2;doublexx,*x,*y,*a,*b,*a1,*b1,*h,*m;cout请输入插值节点个数n：endl;cinn;x=newdouble[n];y=newdouble[n];a=newdouble[n];b=newdouble[n];a1=newdouble[n];b1=newdouble[n];h=newdouble[n-1];m=newdouble[n+1];cout请输入n个插值的节点(xi，yi)：endl;for(inti=0;in;i++){cinx[i]y[i];}for(intj=0;jn-1;j++){h[j]=x[j+1]-x[j];}cout请输入待估点xx：endl;cinxx;cout请选择边界条件：endl;cinchoice;switch(choice){case1:{doubletemp1,temp2;a[0]=0;a[n-1]=1;cout请输入边界条件的两个一阶微商值s'(x1)与s'(xn)：endl;cintemp1temp2;b[0]=2*temp1;b[n-1]=2*temp2;break;}2case2:{a[0]=1;a[n-1]=0;b[0]=3/h[0]*(y[1]-y[0]);b[n-1]=3/h[n-2]*(y[n-1]-y[n-2]);break;}}for(intk=1;kn-1;k++){a[k]=h[k-1]/(h[k-1]+h[k]);b[k]=3*((1-a[k])/h[k-1]*(y[k]-y[k-1])+a[k]/h[k]*(y[k+1]-y[k]));}a1[0]=-a[0]/2;b1[0]=b[0]/2;for(intl=1;ln;l++){a1[l]=-a[l]/(2+(1-a[l])*a1[l-1]);b1[l]=(b[l]-(1-a[l])*b1[l-1])/(2+(1-a[l])*a1[l-1]);}m[n]=0;for(j=n-1;j=0;j--){m[j]=a1[j]*m[j+1]+b1[j];}//判别xx所在区间并输出结果cout\n插值结果为：;for(k=0;kn-1;k++){if(x[k]=xx&&x[k+1]xx){doubleoutput=0;output=(1+2*(xx-x[k])/(x[k+1]-x[k]))*pow(((xx-x[k+1])/(x[k]-x[k+1])),2)*y[k]+(1+2*(xx-x[k+1])/(x[k]-x[k+1]))*pow(((xx-x[k])/(x[k+1]-x[k])),2)*y[k+1]+(xx-x[k])*pow(((xx-x[k+1])/(x[k]-x[k+1])),2)*m[k]+(xx-x[k+1])*pow(((xx-x[k])/(x[k+1]-x[k])),2)*m[k+1];coutoutputendl;break;}}deletex;deletey;deletea;deleteb;deletea1;deleteb1;deleteh;deletem;}运行结果截图：32.三次样条插值（初值条件2）：P52.10、给定函数的函数表和边界条件，求三次样条插值函数，并求的近似值。函数表0.250.30.390.450.530.50.54770.62450.67080.728源代码：yangtiao.cpp(同上)运行结果截图：3.自动选取步长梯形法：P9747、使用自动选取步长梯形法计算积分的近似值。（给定ε=0.01）源代码：SelfSelLength.cpp#includeiostream.h#includemath.hdoublefun(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;k=n;k++){S=S+fun(R_a+(2*k-1)*h/n);}R1=R0/2+S*h/n;E=fabs(R1-R0);n=2*n;}while(E3*e);returnR1;}voidmain(){doublea,b,e;cout请依次输入待求积分函数的下界a、上界b及精度要求e:endl;cinabe;cout自动选取步长梯形法可求得积分:SelfSelLength(a,b,e)endl;}运行结果截图：4.Romberg求积法：P9758、使用Romberg求积法计算积分的近似值。（给定ε=0.01，且取）源代码：Romberg.cpp#includeiostream.h#includemath.hstaticdoubleTri[128][128];doublefun(doublea){returnsqrt(a);}doubleRomberg(doubleR_a,doubleR_b,doublee){Tri[0][0]=(R_b-R_a)/2*(fun(R_a)+fun(R_b));intk=0;doubleE;do//每当误差值不符合要求时，计算下一行的Tri[][]值{k++;doubletemp=0;//计算T[0][k]的数值for(inti=1;i=pow(2,k-1);i++){temp=temp+fun(R_a+(2*i-1)*(R_b-R_a)/pow(2,k));}Tri[0][k]=0.5*(Tri[0][k-1]+(R_b-R_a)/pow(2,k-1)*temp);for(intm=1;m=k;m++){Tri[m][k-m]=(pow(4,m)*Tri[m-1][k-m+1]-Tri[m-1][k-m])/(pow(4,m)-1);}E=fabs(Tri[k][0]-Tri[k-1][0]);}while(Ee);returnTri[k][0];}voidmain(){doublea,b,e;cout请依次输入待求积分函数的下界a、上界b及精度要求e:endl;cinabe;cout按Romberg求积法可求得积分:Romberg(a,b,e)endl;}运行结果截图：5.列主元高斯消去法：6测试矩阵为：=源代码：Guess_Elimination.cpp//列主元高斯消去法#includeiostream.h#includemath.h#defineN4//矩阵的维数，可按需更改staticdoubleA[N][N]={2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2};//系数矩阵staticdoubleB[N]={1,0,1,0};//右端项staticdoubleX[N];inti,j,k;//计数器voidmain(){for(k=0;kN-1;k++){//选取最大主元intindex=k;for(i=k;iN;i++){if(fabs(A[index][k])fabs(A[i][k])){index=i;}}//交换行doubletemp;for(i=k;iN;i++){temp=A[index][i];A[index][i]=A[k][i];A[k][i]=temp;}temp=B[index];B[index]=B[k];B[k]=temp;for(i=k+1;iN;i++){doubleT=A[i][k]/A[k][k];B[i]=B[i]-T*B[k];for(j=k+1;jN;j++){A[i][j]=A[i][j]-T*A[k][j];}}}X[N-1]=B[N-1]/A[N-1][N-1];for(i=N-2;i=0;i--){7doubleTemp=0;for(intj=i+1;jN;j++)Temp=Temp+A[i][j]*X[j];X[i]=(B[i]-Temp)/A[i][i];}cout线性方程组的解（X1，X2，X3......Xn）为：endl;for(i=0;iN;i++){coutX[i];}}运行结果截图：6.列主元LU分解法：测试矩阵为：=源代码：LU_Decomposition.cpp#includeiostream.h#includemath.h#defineN4//矩阵维数，可自定义staticdoubleA[N][N];//系数矩阵staticdoubleB[N];//右端项staticdoubleY[N];//中间项staticdoubleX[N];//输出staticdoubleS[N];//选取列主元的比较器inti,j,k;//计数器voidmain(){cout请输入线性方程组（ai1，ai2，ai3......ain,yi）：endl;for(i=0;iN;i++){for(intj=0;jN;j++)cinA[i][j];cinB[i];}for(k=0;kN;k++){//选列主元intindex=k;8for(i=k;iN;i++){doubletemp=0;for(intm=0;mk;m++){temp=temp+A[i][m]*A[m][k];}S[i]=A[i][k]-temp;if(S[index]S[i]){index=i;}}//交换行doubletemp;for(i=k;iN;i++){temp=A[index][i];A[index][i]=A[k][i];A[k][i]=temp;}temp=B[index];B[index]=B[k];B[k]=temp;//构造L、U矩阵for(j=k;jN;j++){doubletemp=0;for(intm=0;mk;m++){temp=temp+A[k][m]*A[m][j];}A[k][j]=A[k][j]-temp;//先构造U一行的向量}for(i=k+1;iN;i++){doubletemp=0;for(intm=0;mk;m++){temp=temp+A[i][m]*A[m][k];}A[i][k]=(A[i][k]-temp)/A[k][k];//再构造L一列的向量}}//求解LY=BY[0]=B[0];for(i=1;iN;i++){doubletemp=0;for(intj=0;ji;j++)9{temp=temp+A[i][j]*Y[j];}Y[i]=B[i]-temp;}//求解UX=YX[N-1]=Y[N-1]/A[N-1][N-1];for(i=N-2;i=0;i--){doubletemp=0;for(intj=i+1;jN;j++){temp=temp+A[i][j]*X[j];}X[i]=(Y[i]-temp)/A[i][i];}//打印Xcout线性方程组的解（X1，X2，X3......Xn）为：endl;for(i=0;iN;i++){coutX[i];}}运行结果截图：7.简单迭代法（Jacobi迭代）：测试矩阵为：=取精度为esp=0.001，最