数值分析计算实习作业一

《数值分析B》大作业一ZY1515105 樊雪松一．算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] 。

在数组MatrixC[5][501]中检索A的带内元素a ij的方法是：A的带内元素a ij=C中的元素c i-j+2,j。

2.求解λ1，λ501，λs1、首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。

λmin即为λs；如果λ max>0,则λ501=λmax；如果λmax<0,则λ1=λmax。

2、使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ’max，如果λ max>0,则λ1=λ’max+p；如果λmax<0,则λ501=λ’max+p。

3、求解A的与数μk=λ1+k（λ501-λ1）/40 的最接近的特征值λik （k=1,2，…，39）。

使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λ ik。

4、求解A的（谱范数）条件数cond（A）2和行列式detA。

cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。

求解矩阵A的行列式，可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。

二．源程序#include<stdio.h>#include<math.h>#include<conio.h>//定义A中元素double C[5][501];double a[501];double b;double c;//声明所有函数void YaSuoJZ(double C[5][501],double a[501],double b,double c) ;//压缩矩阵函数double mifa(double C[5][501]); //幂法函数void daizhuangLU(double A[5][501]); //带状矩阵的LU分解double fanmifa(double C[5][501]);//反幂法函数//最值函数int max2(int x,int y);int max3(int x,int y,int z);int min(int x,int y);//最值函数int max2(int x,int y) //求2个数的最大值{int z;z=x>y?x:y;return(z);}int max3(int x,int y,int z) //求3个数的最大值{int w;w = z > max2(x,y)? z:max2(x,y);return(w);}int min(int x,int y) //求2个数的最小值{int z;z=x>y?y:x;return(z);}//将矩阵A压缩存储在矩阵C中void YaSuoJZ(double C[5][501],double a[501],double b,double c) {int i;for(i=0;i<=500;i++){if(i>=2) C[0][i]=c;else C[0][i]=0;if(i>=1) C[1][i]=b;else C[1][i]=0;if(i<=499) C[3][i]=b;else C[3][i]=0;if(i<=498) C[4][i]=c;else C[4][i]=0;C[2][i]=a[i];}}//幂法函数：用幂法求矩阵模最大的特征值double mifa(double C[5][501]){double u[501];double y[501]={0},η=0;double β,βk=0;double ε=1;// ε为精度double sumu=0,sumAY=0;int i,j,k=1; //k为循环次数for (i=0;i<=500;i++) //取任一非零向量u0u[i] = 1.0;while(ε>=1e-12){for(i=0;i<=500;i++) //求u(k-1)的2范数ηsumu=sumu+u[i]*u[i];η=sqrt(sumu);sumu=0;for(i=0;i<=500;i++) //求y(k-1)y[i]=u[i]/η;for(i=0;i<=500;i++) //求u(k)的各分量u[i]{for(j=max2(0,i-2);j<=min(i+2,500);j++)sumAY=sumAY+C[i-j+2][j]*y[j];u[i]=sumAY;sumAY=0;}//求幂法中的βkβ=βk; //将β(k-1)放在β中βk=0;for(i=0;i<=500;i++) //求βkβk=βk+y[i]*u[i];if(k>=2)ε=fabs(βk-β)/fabs(βk);k++;}return(βk);}//带状矩阵的LU分解void daizhuangLU(double A[5][501]){int i,j,k,m,t;double sumukj=0,sumlik=0;for(k=0;k<=500;k++){for(j=k;j<=min(k+2,500);j++) //求ukj并存在A[k-j+2][j]中{for(t=max3(0,k-2,j-2);t<=k-1;t++)sumukj=sumukj+A[k-t+2][t]*A[t-j+2][j];A[k-j+2][j]=A[k-j+2][j]-sumukj;sumukj=0;}if(k<500)for(i=k+1;i<=min(k+2,500);i++) //求lik并存在A[i-k+2][k]中{for(m=max3(0,i-2,k-2);m<=k-1;m++)sumlik=sumlik+A[i-m+2][m]*A[m-k+2][k];A[i-k+2][k]=(A[i-k+2][k]-sumlik)/A[2][k];sumlik=0;}}}//反幂法函数：用反幂法求矩阵的模最小的特征值double fanmifa(double M[5][501]){double u[501];double y[501]={0},x[501],η=0;double fβ,fβk=0;double ε=1;double fsumu=0,sumLX=0,sumUu=0;int i,t,m,k=1;for(i=0;i<=500;i++) //任取一非零向量u0u[i]=1;daizhuangLU(M); //对A进行LU分解A=LU，Au(k)=y(k-1)等价于Uu(k)=x和Lx=y(k-1) while(ε>=1e-12){for(i=0;i<=500;i++) //求u(k-1)的2范数ηfsumu=fsumu+u[i]*u[i];η=sqrt(fsumu);fsumu=0;for(i=0;i<=500;i++) //求y(k-1)y[i]=u[i]/η;for(i=0;i<=500;i++) //求中间向量xx[i]=y[i];for(i=1;i<=500;i++){for(t=max2(0,i-2);t<=i-1;t++)sumLX=sumLX+M[i-t+2][t]*x[t];x[i]=x[i]-sumLX;sumLX=0;}u[500]=x[500]/C[2][500]; //求u(k)的各分量u[i]for(i=499;i>=0;i--){for(m=i+1;m<=min(i+2,500);m++)sumUu=sumUu+M[i-m+2][m]*u[m];u[i]=(x[i]-sumUu)/M[2][i];sumUu=0;}//求反幂法中的βkfβ=fβk; //将fβ(k-1)放在fβ中fβk=0;for(i=0;i<=500;i++) //求fβkfβk=fβk+y[i]*u[i];if(k>=2)ε=fabs(1/fβk-1/fβ)/fabs(1/fβk);k++;}return(1/fβk);}//主函数void main(){int i,j,k;double λ1,λ501,λm,λm1,λm2,λs,λ,p;double cond,detA=1;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);b=0.16;c=-0.064;YaSuoJZ(C,a,b, c); //将矩阵A中元素压缩存储在C中λm1=mifa(C); //对A用幂法求出模最大的特征值λm1λs=fanmifa(C); //对A用反幂法求出模最小的特征值λsYaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中for(j=0;j<=500;j++) //对A进行平移，平移量为λm1，平移后矩阵元素压缩存储在C中C[2][j]=C[2][j]-λ?m1;λm=mifa(C);λm2=λm1+λm; //λm1与λm2是矩阵的最大最小特征值if(λm1>λm2) //判断A最大最小特征值{λ501=λm1;λ1=λm2;}else{λ501=λm2;λ1=λm1;}printf("数值分析计算实习第一题\n\n ZY1515105 樊雪松\n\n （1）A的最大最小以及模最小的特征值\n");printf("A的最小特征值λ1=%.13e\n",λ1);printf("A的最大特征值λ501=%.13e\n",λ501);printf("A的模最小特征值λs=%.13e\n",λs);printf("\n（2）与数μk最接近的特征值\n");printf("\t要求接近的值\t\t\t实际求得的特征值\n");YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中for(k=1;k<=39;k++){p=λ1+k*(λ501-λ1)/40;for(j=0;j<=501;j++)C[2][j]=C[2][j]-p;λ=fanmifa(C)+p;printf("μ%d=%.13e λ%d=%.13e\n",k,p,k,λ);YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中}printf("\n（3）计算A的条件数cond(A)和行列式detA\n");cond=λm1/λs;daizhuangLU(C);for(j=0;j<=500;j++)detA=detA*C[2][j];printf("A的条件数cond(A)=%.13e\n",cond);printf("A的行列式detA=%.13e\n",detA);getch();}三、运行结果数值分析计算实习第一题ZY1515105 樊雪松（1）A的最大最小以及模最小的特征值A的最小特征值λ1=-1.0700113615018e+001A的最大特征值λ501=9.7246340987773e+000A的模最小特征值λs=-5.5579107942295e-003（2）与数μk最接近的特征值要求接近的值实际求得的特征值μ1=-1.0189494922173e+001 λ1=-1.0182934033146e+001 μ2=-9.6788762293280e+000 λ2=-9.5857074250676e+000 μ3=-9.1682575364831e+000 λ3=-9.1726724239280e+000 μ4=-8.6576388436383e+000 λ4=-8.6522840078976e+000 μ5=-8.1470201507934e+000 λ5=-8.0934838086753e+000 μ6=-7.6364014579485e+000 λ6=-7.6594054076924e+000 μ7=-7.1257827651036e+000 λ7=-7.1196846486912e+000 μ8=-6.6151640722588e+000 λ8=-6.6117643393973e+000 μ9=-6.1045453794139e+000 λ9=-6.0661032265951e+000 μ10=-5.5939266865690e+000 λ10=-5.5851010526284e+000 μ11=-5.0833079937241e+000 λ11=-5.1140835298122e+000 μ12=-4.5726893008792e+000 λ12=-4.5788721768651e+000 μ13=-4.0620706080344e+000 λ13=-4.0964709262599e+000 μ14=-3.5514519151895e+000 λ14=-3.5542112157508e+000 μ15=-3.0408332223446e+000 λ15=-3.0410900181333e+000 μ16=-2.5302145294997e+000 λ16=-2.5339703111304e+000 μ17=-2.0195958366549e+000 λ17=-2.0032307695635e+000μ18=-1.5089771438100e+000 λ18=-1.5035576112274e+000μ19=-9.9835845096511e-001 λ19=-9.9355860600754e-001μ20=-4.8773975812023e-001 λ20=-4.8704267388496e-001μ21=2.2878934724645e-002 λ21=2.2317362495748e -002μ22=5.3349762756952e-001 λ22=5.3241747420686e -001μ23=1.0441163204144e+000 λ23=1.0528989626935e+000μ24=1.5547350132593e+000 λ24=1.5894458818809e+000μ25=2.0653537061042e+000 λ25=2.0603304602743e+000μ26=2.5759723989490e+000 λ26=2.5580755970728e+000μ27=3.0865910917939e+000 λ27=3.0802405093071e+000μ28=3.5972097846388e+000 λ28=3.6136208676923e+000μ29=4.1078284774837e+000 λ29=4.0913785104506e+000μ30=4.6184471703285e+000 λ30=4.6030353782791e+000μ31=5.1290658631734e+000 λ31=5.1329242838984e+000μ32=5.6396845560183e+000 λ32=5.5949063480833e+000μ33=6.1503032488632e+000 λ33=6.0809338570269e+000μ34=6.6609219417080e+000 λ34=6.6803540921116e+000μ35=7.1715406345529e+000 λ35=7.2938774481266e+000μ36=7.6821593273978e+000 λ36=7.7171117142356e+000μ37=8.1927780202427e+000 λ37=8.2252200140502e+000μ38=8.7033967130876e+000 λ38=8.6486660651935e+000μ39=9.2140154059324e+000 λ39=9.2542003445750e+000（3）计算A 的条件数cond(A)和行列式detAA 的条件数cond(A)=1.9252042739022e+003A 的行列式detA=2.7727861417521e+118四、结果分析设A 的n 个线性无关的特征向量为1x ，2x ,…,n x ，其相对应的特征值满足的关系为n λλλλ≥≥≥> 321。

数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。

再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。

2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。

3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。

4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。

二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。

2019-2020 第1学期数值分析上机实习题总目标：会算，要有优化意识。

（以下程序要求以附件1例题代码格式给出）1. 对给定的线性方程组Ax b =进行迭代求解。

（1）给出Jacobi 迭代的通用程序。

（2）给出Gauss-Seidel 迭代的通用程序。

调用条件：系数矩阵A ，右端项b ，初值0x ，精度要求ε。

输出结果：方程组的近似解。

给定线性方程组211122241125x --⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪---⎝⎭⎝⎭，和122711122215x -⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，取初值0x 为0， 分别利用Jacobi 迭代和G-S 迭代进行求解，观察并解释其中的数学现象。

2. 利用紧凑格式（即直接分解法或逐框运算法）对给定的矩阵A 进行Doolittle 分解，并用其求线性方程组的解。

调用条件：矩阵A 。

输出结果：单位下三角矩阵L 和上三角矩阵U 。

给定矩阵1112A ⎛⎫= ⎪⎝⎭，利用以下算法：1）将A 作Doolittle 分解11A LU =，2）令211A U L =，并对2A 作Doolittle 分解222A L U =，3）重复2）的过程令11n n n A U L --=，并对n A 作Doolittle 分解n n n A L U =，2,3,4,n =， 观察n L ，n U ，n A 的变化趋势，思考其中的数学现象。

3. 给定函数21(),12511f x x x -≤+≤=，取164,8,n =，用等距节点21,i i n x =-+ 0,1,,1i n =+对原函数进行多项式插值和五次多项式拟合，试画出插值和拟合曲线，并给出数学解释。

4. 给出迭代法求非线性方程()0f x =的根的程序。

调用条件：迭代函数()x ϕ，初值0x输出结果：根的近似值k x 和迭代次数k给定方程32()10f x x x =--=，用迭代格式1k x +=0 1.5x =附近的根，要使计算结果具有四位有效数字，利用估计式*1||1||k k k L x x x x L -≤---，或估计式*10||1||kk L x x x x L-≤--来判断需要的迭代次数，分别需要迭代多少次？两者是否有冲突？5. 利用数值求积算法计算()ba f x dx ⎰。

数值分析计算实习题-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN《数值分析》计算实习题姓名：学号：班级：第二章1、程序代码Clear;clc;x1=[ ];y1=[ ];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i-1)*df(i);endP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=;q1=q(1,:)*[^3;^2;;1];q1=vpa(collect(q1),5)q2=q(1,:)*[^3;^2;;1];q2=vpa(collect(q2),5)q3=q(1,:)*[^3;^2;;1];q3=vpa(collect(q3),5)q4=q(1,:)*[^3;^2;;1];q4=vpa(collect(q4),5)%求解并化简多项式2、运行结果P4 =*x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - + q1 =- *x^3 + *x^2 - *x +q2 =- *x^3 + *x^2 - *x + q3 =- *x^3 + *x^2 - *x + q4 =- *x^3 + *x^2 - *x +3、问题结果4次牛顿差值多项式4()P x = *x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - +三次样条差值多项式()Q x0.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1323232321.33930.803570.40714 1.04,[0.2,0.4]1.3393 1.60710.88929 1.1643,[0.4,0.6]1.3393 2.4107 1.6929 1.4171,[0.6,0.8]1.3393 3.21432.8179 1.8629,[0.8,1.0]x x x x x x x x x x x x x x x x ⎧-+-+∈⎪-+-+∈⎪⎨-+-+∈⎪⎪-+-+∈⎩第三章1、程序代码Clear;clc; x=[0 1]; y=[1 ];p1=polyfit(x,y,3)%三次多项式拟合 p2=polyfit(x,y,4)%四次多项式拟合 y1=polyval(p1,x);y2=polyval(p2,x);%多项式求值plot(x,y,'c--',x,y1,'r:',x,y2,'y-.')p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。

《数值分析》计算作业院系：航空科学与工程学院学号： SY1005512姓名：王天龙日期： 2010年10月31日计算实习说明书目的：训练运用计算机进行科学与工程计算的能力。

要求：1．独立进行算法设计、程序设计和上机运算，并得出正确的结果。

2．编制程序时全部采用双精度，要求按题目的要求设计输出，并执行打印。

3．只能根据题目给出的信息并且只允许一次计算得出全部结果。

题目：第一题 设有501×501的矩阵123499500501a b c b a b cc b a b c A c b a b c c b a b c ba ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦其中0.1(1.640.024)s i n (0.2)0.64 (125i i a i i e i =--= ，，，；0.16b =；0.064c =-。

矩阵A 的特征值12501()ii λ= ，，，满足 125011501||min ||S i i λλλλλ≤≤<<<= ，试求：1.1λ，501λ和S λ的值。

2.A 的与数5011140k kλλμλ-=+最接近的特征值(1239)ik k λ= ，，，。

3.A 的（谱范数）条件数2()cond A 和行列式det A 。

说明：1．在所有的算法中，凡是要给出精度水平ε的，都取1210ε-=。

2．选择算法时，应使A 的所有零元素都不存储。

3．打印以下内容： （1）算法的设计方案。

（2）全部源程序（要求注明主程序和每个子程序的功能）。

（3）特征值1λ，501λ，S λ和(1239)ik k λ= ，，，以及2()cond A ，det A 的值。

4.采用e 型输出所有计算结果，并至少显示12位有效数字。

一、程序算法的设计算法设计方案如下：二、全部源程序编程软件：Fortran：三、计算结果1.特征值1λ，501λ，S λ1-.107001135582E+02λ=，501 .972463398616E+01λ=，-.555823879237E-02S λ=2. (1239)ik k λ= ，，，如下表所示（ZK 代表ik λ）3.A 的条件数2()cond A 和行列式det A 的值2() .192509100000E+04cond A =，det .277059968428+119A =五、讨论这里选取的初始向量为X(i)=1,X={x1,x2,x3,,,,,,,x501},当初始向量与特征向量较近时，收敛较快，若初始向量与特征向量正交，则求解可能失真。

插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-24221063/63504000*x^2+95549/72072*x-1/3048192000*x^8-2168879/43545600 0*x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395 008000*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数 S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x) 输出结果为：S =23491/304472833/8*x+76713/*x^2+6867/42624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线010203040506070-2020406080100可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。

《数值分析》课程设计—作业实验一1.1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。

由于旅途的颠簸，大家都很疲惫，很快就入睡了。

第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。

第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。

解：一、问题分析：对于本题，比较简单，我们只需要判断原来椰子的个数及每个人私藏了一份之后剩下的是否能被5除余1，直到最后分完。

二、问题求解：通过matlab 建立M 文件，有如下程序：或者对于第一个程序，n 取2000；对于第二个程序，n 取20001，就能得到我们想要的结果，即原先一共有15621个椰子，最终平均每人得4092个椰子。

n=input('input n:');forx=1:n p=5*x+1;for k=1:5 p=5*p/4+1;end if p==fix(p) break ; end enddisp([x,p]) input n:20001023 15621function fentao(n)a=cat(1,7);for j=n:-1:1 a(1)=j;i=1; while i<7a(i+1)=4*(a(i)-1)/5; i=i+1;endif a(7)==fix(a(7)) a, end end end>> fentao(20001) a =15621 12496 9996 7996 6396 5116 4092（本文档内的有些运行结果，限于篇幅，使文档结构更和谐、紧凑，已做相关的改动，程序代码没变）1.2 当0,1,2,,100n = 时，选择稳定的算法计算积分1d 10n xn xe I x e--=+⎰.解：一、问题分析：由1d 10n xn xe I x e--=+⎰知： 110101==+⎰dx I I以及： )1(110101011)1(1nnxxnxxn n n endx edx ee eI I ----+-+-==++=+⎰⎰得递推关系：⎪⎩⎪⎨⎧--=-=-+n n n I e n I I I 10)1(1101101，但是通过仔细观察就能知道上述递推公式每一步都将误差放大十倍，即使初始误差很小，但是误差的传播会逐步扩大，也就是说用它构造的算法是不稳定的，因此我们改进上述递推公式（算法）如下：⎪⎪⎩⎪⎪⎨⎧--=-=+-))1(1(101)1(101110n n n I e n I I I通过比较不难得出该误差是逐步缩小的，即算法是稳定的。

问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。

分析：本问题是已知五个点，由这五个点求一三次样条插值函数。

边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。

对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。

⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。

《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det 要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。

变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。

对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。

如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。

一、生成矩阵1.5个100阶对角方阵。

根据下述问题，设计5个100阶对角矩阵，对角元分布如下：表1.1 5个对角阵的对角元分布生成的5个对角矩阵中，前四个矩阵为正定矩阵，第五个矩阵为非正定矩阵。

其中，前三个矩阵条件数都为100。

2.生成10个100阶矩阵生成10个100阶矩阵，做QR分解，生成10个Q矩阵。

由于生成的Q,即对应的特征向量，对最后的求解结果影响不大，故此处随机生成了10个100阶矩阵，且为了每次的结果一样，将矩阵存到文件中。

一、实验结果（1）用共轭梯度法，Lanczos法，MINRES法进行计算，分析特征值分布对算法的影响。

计算得到的结果如下表所示：其中的为e w最后一次迭代之后求得的解与真实值的差，即e w=x w−x exact.从表格中可以看出以下几点：●在矩阵特征值相同时，特征值的分布对收敛性有影响。

在特征值偏向较小特征值的时候，CG法和Minres算法迭代次数增加，而Lanczos算法更是出现了不收敛的情况。

相对的，当特征值偏向较大特征值的时候，三种算法的迭代次数都均匀分布的有所下降。

故可有结论，特征值偏向最大值时有较好收敛性。

●矩阵条件数对收敛性有影响。

由三种算法的第3，4组数据能看出，3对应矩阵的条件数大于4矩阵的条件数，且3的迭代次数均小于4的迭代次数，故可有结论，条件数大的正定矩阵有较好的收敛性。

●矩阵是否正定对收敛性有影响。

由于CG法不适用于非正定矩阵，故只对Minres法和Lanczos法分析。

Lanczos法在实验中的非正定矩阵的运算中不收敛，且Minres法的迭代次数也较之前有了大幅度的提高。

故可有结论，正定矩阵较非正定矩阵有较好的收敛性。

（2）取定一个对角矩阵，用10个单位对称阵做实验，观察特征向量对算法的影响。

此处为了观察特征向量的影响，选定了最普通的1号对角阵。

结果如下表所示：表2.2 共轭梯度法对于不同特征向量的结果表2.4 Minres法对于不同特征向量的结果从以上三个表格数据可以看出，矩阵对于不同特征向量的变化收敛性变化不大，在10个不同特征向量下，迭代次数变化最大在3次之内，可以认为对迭代次数没有影响。

《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。

1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。

2)使用反幂法求λs ，其中需要解线性方程组。

因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。

（2）与140k λλμλ-5011=+k 最接近的特征值λik 。

通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。

（3）2cond(A)和det A 。

1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。

2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。

由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。

2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2] ;}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。

北京航空航天大学《数值分析》计算实习报告第一大题学 院：自动化科学与电气工程学院 专 业： 控制科学与工程 学 生 姓 名： 学 号： 教 师： 电 话： 完 成 日 期： 2015年11月6日北京航空航天大学Beijing University of Aeronautics and Astronautics实习题目：第一题 设有501501⨯的实对称矩阵A ，其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。

矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有1.求1λ，501λ和s λ的值。

2.求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。

3.求A 的(谱范数)条件数2)A (cond 和行列式detA 。

说明：1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。

2.选择算法时，应使矩阵A 的所有零元素都不储存。

3.打印以下内容： （1）全部源程序；（2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。

4.采用e 型输出实型数，并且至少显示12位有效数字。

一、算法设计方案1、求1λ，501λ和s λ的值。

由于||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤，可知绝对值最大特征值必为1λ和501λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则501λ=λ。

将矩阵A 进行一下平移：I -A A'λ= （1）对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。

s λ为按模最小特征值，||min ||5011i i s λλ≤≤=，可对A 使用反幂法求得。

数值分析—计算实习作业一学院：17系专业：精密仪器及机械姓名：张大军学号：DY14171142014-11-11数值分析计算实现第一题报告一、算法方案算法方案如图1所示。

（此算法设计实现完全由本人独立完成）图1算法方案流程图二、全部源程序全部源程序如下所示#include <iostream.h>#include <iomanip.h>#include <math.h>int main(){double a[501];double vv[5][501];double d=0;double r[3];double uu;int i,k;double mifayunsuan(double *a,double weiyi);double fanmifayunsuan(double *a,double weiyi);void yasuo(double *A,double (*C)[501]);void LUfenjie(double (*C)[501]);//赋值语句for(i=1;i<=501;i++){a[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);}//程序一：使用幂方法求绝对值最大的特征值r[0]=mifayunsuan(a,d);//程序二：使用幂方法求求平移λ[0]后绝对值最大的λ，得到原矩阵中与最大特征值相距最远的特征值d=r[0];r[1]=mifayunsuan(a,d);//比较λ与λ-λ[0]的大小，由已知得if(r[0]>r[1]){d=r[0];r[0]=r[1];r[1]=d;}//程序三：使用反幂法求λr[2]=fanmifayunsuan(a,0);cout<<setiosflags(ios::right);cout<<"λ["<<1<<"]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[0]<<endl;cout<<"λ["<<501<<"]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[1]<<endl;cout<<"λ[s]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[2]<<endl;//程序四：求A的与数u最接近的特征值for(k=1;k<40;k++){uu=r[0]+k*(r[1]-r[0])/40;cout<<"最接近u["<<k<<"]"<<"的特征值为"<<setiosflags(ios::scientific)<<setprecision(12)<<fanmifayunsuan(a,uu)<<endl;}//程序五：谱范数的条件数是绝对值最大的特征值除以绝对值最小的特征值的绝对值cout<<"cond(A)2="<<fabs(r[0]/r[2])<<endl;//程序六：A的行列式的值就是A分解成LU之U的对角线的乘积yasuo(a,vv);LUfenjie(vv);uu=1;for(i=0;i<501;i++){uu=uu*vv[2][i];}cout<<"Det(A)="<<uu<<endl;return 1;}double mifayunsuan(double *a,double weiyi){int i,k;double b=0.16;double c=-0.064;double ee,w,v1,v2,mm,sum;double u[501];double y[505]={0};for(i=0;i<501;i++)u[i]=1;//给u赋初值if (weiyi!=0){for (i=0;i<501;i++)a[i]-=weiyi;}ee=1;k=0;//使得初始计算时进入循环语句while(ee>1e-12){mm=0;for(i=0;i<501;i++){mm=mm+u[i]*u[i];}w=sqrt(mm);for(i=0;i<501;i++){y[i+2]=u[i]/w;//注意此处编程与书上不同，之后会解释它的巧妙之处1 }for(i=0;i<501;i++){u[i]=c*y[i]+b*y[i+1]+a[i]*y[i+2]+b*y[i+3]+c*y[i+4];//1显然巧妙之处凸显出来}sum=0;for(i=0;i<501;i++){sum+=y[i+2]*u[i];}v1=v2;v2=sum;//去除特殊情况，减少漏洞if(k==0){k++;}else{ee=fabs(v2-v1)/fabs(v2);}}if (weiyi!=0){for (i=0;i<501;i++)a[i]+=weiyi;}//还原A矩阵return (v2+weiyi);}double fanmifayunsuan(double *a,double weiyi){int i,k;double b=0.16;double c=-0.064;double ee,w,v1,v2,mm,sum;double u[501];double y[501];double C[5][501];void yasuo(double *A,double (*C)[501]);void LUfenjie(double (*C)[501]);void qiuU(double (*C)[501],double *y,double *u);//把A阵压缩到C阵中for(i=0;i<501;i++)u[i]=1;//给u赋初值if (weiyi!=0){for (i=0;i<501;i++)a[i]-=weiyi;}yasuo(a,C);LUfenjie(C);ee=1;k=0; //使得初始计算时进入循环语句while(ee>1e-12){mm=0;for(i=0;i<501;i++){mm=mm+u[i]*u[i];}w=sqrt(mm);for(i=0;i<501;i++){y[i]=u[i]/w;}qiuU(C,y,u);sum=0;for(i=0;i<501;i++){sum+=y[i]*u[i];}v1=v2;v2=sum;//去除特殊情况，减少漏洞if(k==0){k++;}else{ee=fabs(1/v2-1/v1)/fabs(1/v2);}}if (weiyi!=0){for (i=0;i<501;i++)a[i]+=weiyi;}//还原A矩阵return (1/v2+weiyi);}void yasuo(double *A,double (*C)[501]){double b=0.16;double c=-0.064;int i;for(i=0;i<501;i++){C[0][i]=c;C[1][i]=b;C[2][i]=A[i];C[3][i]=b;C[4][i]=c;}}void LUfenjie(double (*C)[501]){int k,t,j;int r=2,s=2;double sum;int minn(int ,int );int maxx(int ,int );for(k=0;k<501;k++){for(j=k;j<=minn(k+s,501-1);j++){if(k==0)sum=0;else{sum=0;for(t=maxx(k-r,j-s);t<k;t++){sum=sum+C[k-t+s][t]*C[t-j+s][j];}}C[k-j+s][j]=C[k-j+s][j]-sum;}for(j=k+1;j<=minn(k+r,501-1);j++){if(k<501-1){if(k==0)sum=0;else{sum=0;for(t=maxx(j-r,k-s);t<k;t++){sum=sum+C[j-t+s][t]*C[t-k+s][k];}}C[j-k+s][k]=(C[j-k+s][k]-sum)/C[s][k];}}}}void qiuU(double (*C)[501],double *y,double *u){int i,t;double b[501];double sum;int r=2,s=2;int minn(int ,int );int maxx(int ,int );for(i=0;i<501;i++){b[i]=y[i];}for(i=1;i<501;i++){sum=0;for(t=maxx(0,i-r);t<i;t++){sum=sum+C[i-t+s][t]*b[t];}b[i]=b[i]-sum;}u[500]=b[500]/C[s][500];for(i=501-2;i>=0;i--){sum=0;for(t=i+1;t<=minn(i+s,500);t++){sum=sum+C[i-t+s][t]*u[t];}u[i]=(b[i]-sum)/C[s][i];}}int minn(int x,int y){int min;if(x>y)min=y;elsemin=x;return min;}int maxx(int b,int c){int max;if(b>c){if(b>0)max=b;elsemax=0;}else{if(c>0)max=c;elsemax=0;}return max;}三、特征值以及的值λ[1]=-1.070011361502e+001 λ[501]=9.724634098777e+000λ[s]=-5.557910794230e-003最接近u[1]的特征值为-1.018293403315e+001最接近u[2]的特征值为-9.585707425068e+000最接近u[3]的特征值为-9.172672423928e+000最接近u[4]的特征值为-8.652284007898e+000最接近u[5]的特征值为-8.0934********e+000最接近u[6]的特征值为-7.659405407692e+000最接近u[7]的特征值为-7.119684648691e+000最接近u[8]的特征值为-6.611764339397e+000最接近u[9]的特征值为-6.0661********e+000最接近u[10]的特征值为-5.585101052628e+000最接近u[11]的特征值为-5.114083529812e+000最接近u[12]的特征值为-4.578872176865e+000最接近u[13]的特征值为-4.096470926260e+000最接近u[14]的特征值为-3.554211215751e+000最接近u[15]的特征值为-3.0410********e+000最接近u[16]的特征值为-2.533970311130e+000最接近u[17]的特征值为-2.003230769563e+000最接近u[18]的特征值为-1.503557611227e+000最接近u[19]的特征值为-9.935586060075e-001最接近u[20]的特征值为-4.870426738850e-001最接近u[21]的特征值为2.231736249575e-002最接近u[22]的特征值为5.324174742069e-001最接近u[23]的特征值为1.052898962693e+000最接近u[24]的特征值为1.589445881881e+000最接近u[25]的特征值为2.060330460274e+000最接近u[26]的特征值为2.558075597073e+000最接近u[27]的特征值为3.080240509307e+000最接近u[28]的特征值为3.613620867692e+000最接近u[29]的特征值为4.0913********e+000最接近u[30]的特征值为4.603035378279e+000最接近u[31]的特征值为5.132924283898e+000最接近u[32]的特征值为5.594906348083e+000最接近u[33]的特征值为6.080933857027e+000最接近u[34]的特征值为6.680354092112e+000最接近u[35]的特征值为7.293877448127e+000最接近u[36]的特征值为7.717111714236e+000最接近u[37]的特征值为8.225220014050e+000最接近u[38]的特征值为8.648666065193e+000最接近u[39]的特征值为9.254200344575e+000cond(A)2=1.925204273902e+003 Det(A)=2.772786141752e+118四、现象讨论在大作业的程序设计过程当中，初始向量的赋值我顺其自然的设为第一个分量为1，其它分量为0的向量，计算结果与参考答案存在很大差别，计算结果对比如下图2所示（左侧为正确结果，右侧为错误结果），导致了我花了很多的时间去检查程序算法。

数值分析大作业第一题学号：姓名：学院：授课教师：2016.10.24一、算法设计方案1、求λ1，λ501和λsa)利用幂法计算出矩阵A 按模最大的特征值，设其为λm 。

b)令矩阵B =A −λm I （I 为单位矩阵），同样利用幂法计算出矩阵B 按模最大的特征值λm ′。

c)令λm ′′=λm ′+λm 。

由计算过程可知λm 和λm ′′分别为矩阵A 所有特征值按大小排序后，其中两端的值。

即，λ1=min⁡{λm ,λm ′′}，λ501=max⁡{λm ,λm ′′}。

d) 利用反幂法计算λs 。

其中通过使用doolittle 三角分解法求解线性方程组AX =I 的形式求得矩阵A 的逆A −1。

进行doolittle 分解，可以节省内存空间而把矩阵A 压缩，将A 中带上的元素存储为数组C [5][501]。

2、求A 的与数μk =λ1+k λ501−λ140最接近的特征值λi k (k =1,2, (39)e) 令矩阵D k =A −μk I ，利用反幂法计算出矩阵D k 按模最小的特征值λi k ′,则λi k =λi k ′+μk 。

其中，同样使用doolittle 三角分解法求解线性方程组D k X =I 的形式求得矩阵D k 的逆D k −1,其中k =1,2, (39)3）求矩阵A 的条件数cond(A )2和行列式det(A )f) 计算矩阵的条件数cond(A)2=|λm λs |。

g)计算矩阵的行列式det(A )为A 三角分解后得到的上三角矩阵U 的对角线元素之积。

二、程序#include<iostream>#include<cmath>#include <iomanip>using namespace std;double LU(double C[6][502]);//LU 分解double solve(double CC[6][502],double array[502],double U[502]);//解线性方程组double mefa(double arr[][502]);//幂法double fanmefa(double array[][502]);//反幂法int max1(int a,int b,int c);int max2(int a ,int b);int min1(int a,int b);int main(){double a[502],d[6][502],d1[6][502],lam1,lam501,temp,lams,lamki,uk[40],condA,detA=1; //d 为压缩矩阵,d1为进行原点平移时求lam1和lam501的压缩矩阵int i,j,k;for(i=1;i<=501;i++){a[i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);//求解主对角线上的元素}for (i=1;i<=5;i++)for (j=1;j<=501;j++){if(i==1){if ((j==1)||(j==2))d[i][j]=0;else d[i][j]=-0.064;}else if (i==2){if (j==1)d[i][j]=0;else d[i][j]=0.16;}else if (i==3)d[i][j]=a[j];else if(i==4){if (j==501)d[i][j]=0;else d[i][j]=0.16;}else if(i==5){if ((j==500)||(j==501))d[i][j]=0;else d[i][j]=-0.064;}}//The purpose of this part is to assign values to matrix d[][],and d[][] is a condensed matrix for A.lam1=mefa(d);for (i=1;i<=5;i++)for (j=1;j<=501;j++){if(i==3)d1[3][j]=d[3][j]-lam1;else d1[i][j]=d[i][j];}//设置d1的值,进行原点平移lam501=mefa(d1);lam501=lam501+lam1;if(lam1>lam501){temp=lam1;lam1=lam501;lam501=temp;}//将lam1和lam501中小的值赋给lam1for(i=1;i<=39;i++)uk[i]=lam1+i*(lam501-lam1)/40;//求解uk的值cout<<"近似特征值的结果为:"<<endl;for(i=1;i<=39;i++){for(k=1;k<=5;k++){for(j=1;j<=501;j++){if(k==3)d1[k][j]=d[k][j]-uk[i];else d1[k][j]=d[k][j];}}//求解uk时的原点平移LU(d1);lamki=fanmefa(d1);cout<<"λk"<<i<<'=';cout<<setiosflags(ios::scientific)<<setprecision(12)<<lamki+uk[i];if(i%2==0)cout<<endl;else cout<<" ";}cout<<endl;LU(d); //进行LU分解lams=fanmefa(d); //进行反幂法condA=fabs(lam1/lams);//求矩阵的条件数for(i=1;i<=501;i++)detA=detA*d[3][i];//求行列式的值cout<<"矩阵的最小特征值为：λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<lam1<<endl;cout<<"矩阵的最大特征值为：λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<lam501<<endl;cout<<"按模最小的特征值为：λs="<<setiosflags(ios::scientific)<<setprecision(12)<<lams<<endl;cout<<"矩阵的条件数为：cond(A)="<<setiosflags(ios::scientific)<<setprecision(12)<<condA<<endl;cout<<"行列式的值为：det(A)="<<setiosflags(ios::scientific)<<setprecision(12)<<detA<<endl;return(0);}double mefa(double arr[][502])//幂法求解按模最大的特征值{double u[502],y[502],eita,beta=0,firstbeta=0,sum,sum1;int i,j,k;for(i=1;i<=501;i++)u[i]=1;for (k=1;k<10000;k++){sum=0;for(j=1;j<=501;j++)sum=sum+u[j]*u[j];eita=sqrt(sum);for(j=1;j<=501;j++)y[j]=u[j]/eita;for(j=1;j<=501;j++){if(j==1)u[j]=arr[3][1]*y[1]+arr[2][2]*y[2]+arr[1][3]*y[3];if(j==2)u[j]=arr[4][1]*y[1]+arr[3][2]*y[2]+arr[2][3]*y[3]+arr[1][4]*y[4];if((j>=3)&&(j<=499))u[j]=arr[5][j-1]*y[j-2]+arr[4][j-1]*y[j-1]+arr[3][j]*y[j]+arr[2][j+1]*y[j+1]+arr[1][j+2]*y[j+2];if(j==500)u[j]=arr[5][498]*y[498]+arr[4][499]*y[499]+arr[3][500]*y[500]+arr[2][501] *y[501];if(j==501)u[j]=arr[5][499]*y[499]+arr[4][500]*y[500]+arr[3][501]*y[501];}//幂法中利用矩阵相乘求解uk的值firstbeta=beta;sum1=0;for(j=1;j<=501;j++)sum1=sum1+y[j]*u[j];beta=sum1;if(fabs((beta-firstbeta)/beta)<=1e-12) break;}return(beta);}double fanmefa(double array[][502])//利用反幂法求解矩阵按模最小的特征值{double u[502],y[502],eita,firstbeta=0,beta=0,sum,sum1;int i,j,k;for(i=1;i<=501;i++)u[i]=1;for (k=1;k<10000;k++){sum=0;for(j=1;j<=501;j++)sum=sum+u[j]*u[j];eita=sqrt(sum);for(j=1;j<=501;j++)y[j]=u[j]/eita;solve(array,y, u); //调用解线性方程组的程序解出ukfirstbeta=beta;sum1=0;for(j=1;j<=501;j++)sum1=sum1+y[j]*u[j];beta=sum1;if(fabs((firstbeta-beta)/firstbeta)<=1e-12) break;//}return(1/beta);}double LU(double C[6][502])//此函数对矩阵进行LU分解{double sum1,sum2;int i,j,k,t;for (k=1;k<=501;k++){for (j=k;j<=min1(k+2,501);j++){sum1=0;for (t=max1(1,k-2,j-2);t<=k-1;t++)sum1=sum1+C[k-t+3][t]*C[t-j+3][j];C[k-j+3][j]=C[k-j+3][j]-sum1;}for(i=k+1;i<=min1(k+2,501);i++){sum2=0;for(t=max1(1,i-2,k-2);t<=k-1;t++)sum2=sum2+C[i-t+3][t]*C[t-k+3][k];C[i-k+3][k]=(C[i-k+3][k]-sum2)/C[3][k];}}return(0);}double solve(double CC[6][502],double array[502],double U[502]) //此函数的功能是解线性方程组{double sum1,sum2,b[502];int i,t;//b[1]=array[1];for(i=1;i<=501;i++)b[i]=array[i];for(i=2;i<=501;i++){sum1=0;for(t=max2(1,i-2);t<=i-1;t++)sum1=sum1+CC[i-t+3][t]*b[t];b[i]=b[i]-sum1;}U[501]=b[501]/CC[3][501];for(i=500;i>=1;i--){sum2=0;for(t=i+1;t<=min1(i+2,501);t++)sum2=sum2+CC[i-t+3][t]*U[t];U[i]=(b[i]-sum2)/CC[3][i];}return(0);}int max1(int a,int b,int c){int t=a;if(b>t)t=b;if(c>t)t=c;return(t);}int max2(int a ,int b){int t=a;if(b>t)t=b;return(t);}int min1(int a,int b) {int t=a;if(b<t)t=b;return(t);}三、计算结果计算结果如下所示：1、矩阵的最小特征值为：λ1=-1.070011361502e+001矩阵的最大特征值为：λ501=9.724634098777e+000按模最小的特征值为：λs=-5.557910794230e-0032、近似特征值的结果为:λk1=-1.018293403315e+001 λk2=-9.585707425068e+000λk3=-9.172672423928e+000 λk4=-8.652284007898e+000λk5=-8.0934********e+000 λk6=-7.659405407692e+000λk7=-7.119684648691e+000 λk8=-6.611764339397e+000λk9=-6.0661********e+000 λk10=-5.585101052628e+000λk11=-5.114083529812e+000 λk12=-4.578872176865e+000λk13=-4.096470926260e+000 λk14=-3.554211215751e+000λk15=-3.0410********e+000 λk16=-2.533970311130e+000λk17=-2.003230769563e+000 λk18=-1.503557611227e+000λk19=-9.935586060075e-001 λk20=-4.870426738850e-001λk21=2.231736249575e-002 λk22=5.324174742069e-001λk23=1.052898962693e+000 λk24=1.589445881881e+000λk25=2.060330460274e+000 λk26=2.558075597073e+000λk27=3.080240509307e+000 λk28=3.613620867692e+000λk29=4.0913********e+000 λk30=4.603035378279e+000λk31=5.132924283898e+000 λk32=5.594906348083e+000λk33=6.080933857027e+000 λk34=6.680354092112e+000λk35=7.293877448127e+000 λk36=7.717111714236e+000λk37=8.225220014050e+000 λk38=8.648666065193e+000λk39=9.254200344575e+0003、矩阵的条件数为：cond(A)=1.925204273902e+003行列式的值为：det(A)=2.772786141752e+118四、讨论：初始向量的选取对计算结果的影响当选取初始向量为(1,1,1,1,...,1),1,2,3,...,501iu i==时，得到以上结果。

1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-/*x^2+95549/72072*x-1/00*x^8-2168879/0*x^4+19/0*x^7+657859/*x^3+33983/ 0*x^5-13003/00*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x)输出结果为：S =/6464-2399/88*x+/1984*x^2+2656867/624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。