北航数值分析全部三次大作业

北航《数值分析》计算实习题目一一、算法设计：要计算矩阵的最大最小特征值，可先通过幂法求得模最大的特征值λa，再根据λa进行原点平移求出下一特征值λb，比较两值的大小，循环下去，最终可以确定最小特征值λ1和最大特征值λ501。要计算矩阵按模最小的特征值λs，可以用反幂法求得，此处需要求解中间过渡向量，采用LU分解法解带状方程组，故要设计Doolite分解。要求解与给定数值接近的特征值，可以对该数做漂移量，新数组特征值倒数的绝对值满足反幂法的要求，故通过反幂法即可求得。要计算cond(A)，可由公式cond(A)=│λmax/λmin│求得，其中λmax和λmin分别矩阵A的模为最大和模为最小的特征值，λmax可以通过幂法求得，而λmin可已通过反幂法求得。要计算detA，可对A作LU分解，即A=LU，则│A│=│L││U│，其中L为单位下三角阵，U为上三角阵，则detA即为U矩阵主对角元素之积。由于矩阵A中的零元素较多，为了节省存储空间，可以将矩阵A中的元素存储在5*501的二维数组中。编译环境：vc++6.0编译函数：幂法，反幂法，Doolite分解二、全部源程序：#include<iostream.h>#include<math.h>#include<stdio.h>#include<stdlib.h>intMax(intv1,intv2);intMin(intv1,intv2);voidzhuancun(doubleA[5][501]);doublemifa(doubleA[5][501]);voiddoolite(doubleA[5][501],doublex[501],doubleb[501]);doublefanmifa(doubleA[5][501]);doubleDet(doubleA[5][501]);voidmain(){ doubleb=0.16,c=-0.064,p,q; inti,j; doubleA[5][501]; zhuancun(A,b,c);//进行A的赋值 cout.precision(12);//定义输出精度 doublelamda1,lamda501,lamdas; doublek=mifa(A); if(k>0)//判断求得最大以及最小的特征值.如果K>0，则它为最大特征值值，//并以它为偏移量再用一次幂法求得新矩阵最大特征值，即为最大//与最小的特征值的差 { lamda501=k; for(i=0;i<=500;i++) A[2][i]=A[2][i]-k; lamda1=mifa(A)+lamda501; for(i=0;i<=500;i++) A[2][i]=A[2][i]+k; } else//如果K<=0，则它为最小特征值值，并以它为偏移量再用一次幂法 //求得新矩阵最大特征值，即为最大与最小的特征值的差{ lamda1=k; for(i=0;i<=500;i++) A[2][i]=A[2][i]-k; lamda501=mifa(A)+lamda1; for(i=0;i<=500;i++) A[2][i]=A[2][i]+k; } lamdas=fanmifa(A); FILE*fp=fopen("result.txt","w"); fprintf(fp,"λ1=%.12e\n",lamda1); fprintf(fp,"λ501=%.12e\n",lamda501); fprintf(fp,"λs=%.12e\n\n",lamdas); fprintf(fp,"\t要求接近的值\t\t\t实际求得的特征值\n"); for(i=1;i<=39;i++)//反幂法求得与给定值接近的特征值 { p=lamda1+(i+1)*(lamda501-lamda1)/40; for(j=0;j<=500;j++) A[2][j]=A[2][j]-p; q=fanmifa(A)+p; for(j=0;j<=500;j++) A[2][j]=A[2][j]+p; fprintf(fp,"μ%d:%.12eλi%d:%.12e\n",i,p,i,q); } doublecond=fabs(mifa(A)/fanmifa(A)); doubledet=Det(A); fprintf(fp,"\ncond(A)=%.12e\n",cond); fprintf(fp,"\ndetA=%.12e\n",det);}/*定义2个判断大小的函数*/intMax(intv1,intv2){ return((v1>v2)?v1:v2);}intMin(intv1,intv2){ return((v1<v2)?v1:v2);}/*将矩阵值转存在一个数组里*/voidzhuancun(doubleA[5][501],doubleb,doublec){ inti=0,j=0; A[i][j]=0,A[i][j+1]=0; for(j=2;j<=500;j++) A[i][j]=c; i++; j=0; A[i][j]=0; for(j=1;j<=500;j++) A[i][j]=b; i++; for(j=0;j<=500;j++) A[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1)); i++; for(j=0;j<=499;j++) A[i][j]=b; A[i][j]=0; i++; for(j=0;j<=498;j++) A[i][j]=c; A[i][j]=0,A[i][j+1]=0;}/*幂法求解模最大的特征值*/doublemifa(doubleA[5][501]){ ints=2,r=2,m=0,i,j; doubleb2,b1=0,sum,u[501],y[501]; for(i=0;i<=500;i++){ u[i]=1.0; } do { sum=0; if(m!=0)b1=b2; m++; for(i=0;i<=500;i++) sum+=u[i]*u[i]; for(i=0;i<=500;i++) y[i]=u[i]/sqrt(sum); for(i=0;i<=500;i++) { u[i]=0; for(j=Max(i-r,0);j<=Min(i+s,500);j++) u[i]=u[i]+A[i-j+s][j]*y[j]; } b2=0; for(i=0;i<=500;i++) b2=b2+y[i]*u[i]; } while(fabs(b2-b1)/fabs(b2)>=exp(-12));//设置计算精度 returnb2;}/*带状DOOLITE分解，并且求解出方程组的解*/voiddoolite(doubleA[5][501],doublex[501],doubleb[501]){ inti,j,k,t,s=2,r=2; doubleB[5][501],c[501]; for(i=0;i<=4;i++) { for(j=0;j<=500;j++) B[i][j]=A[i][j]; } for(i=0;i<=500;i++) c[i]=b[i]; for(k=0;k<=500;k++) { for(j=k;j<=Min(k+s,500);j++) { for(t=Max(0,Max(k-r,j-s));t<=k-1;t++) B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j]; } for(i=k+1;i<=Min(k+r,500);i++) { for(t=Max(0,Max(i-r,k-s));t<=k-1;t++) B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k]; B[i-k+s][k]=B[i-k+s][k]/B[s][k]; } } for(i=1;i<=500;i++) for(t=Max(0,i-r);t<=i-1;t++) c[i]=c[i]-B[i-t+s][t]*c[t]; x[500]=c[500]/B[s][500]; for(i=499;i>=0;i--) { x[i]=c[i]; for(t=i+1;t<=Min(i+s,500);t++) x[i]=x[i]-B[i-t+s][t]*x[t]; x[i]=x[i]/B[s][i]; }}/*反幂法求解模最大的特征值*/doublefanmifa(doubleA[5][501]){ ints=2,r=2,m=0,i; doubleb2,b1=0,sum=0,u[501],y[501]; for(i=0;i<=500;i++) { u[i]=1.0; } do { if(m!=0)b1=b2; m++; sum=0; for(i=0;i<=500;i++) sum+=u[i]*u[i]; for(i=0;i<=500;i++) y[i]=u[i]/sqrt(sum); doolite(A,u,y); b2=0; for(i=0;i<=500;i++) b2+=y[i]*u[i]; } while(fabs(b2-b1)>=fabs(b1)*exp(-12)); return1/b2;}/*矩阵A的LU分解，其中U的主线乘积即为矩阵的DET*/doubleDet(doubleA[5][501]){ inti,j,k,t,s=2,r=2; for(k=0;k<=500;k++) { for(j=k;j<=Min(k+s,500);j++) { for(t=Max(0,Max(k-r,j-s));t<=k-1;t++) A[k-j+s][j]=A[k-j+s][j]-A[k-t+s][t]*A[t-j+s][j]; } for(i=k+1;i<=Min(k+r,500);i++) { for(t=Max(0,Max(i-r,k-s));t<=k-1;t++) A[i-k+s][k]=A[i-k+s][k]-A[i-t+s][t]*A[t-k+s][k]; A[i-k+s][k]=A[i-k+s][k]/A[s][k]; } } doubledet=1; for(i=0;i<=500;i++) det*=A[s][i]; returndet;}三、运行结果:λ1=-1.069936345952e+001λ501=9.722283648681e+000λs=-5.557989086521e-003要求接近的值 实际求得的特征值μ1:-9.678281104107e+000λi1:-9.585702058251e+000μ2:-9.167739926402e+000λi2:-9.172672423948e+000μ3:-8.657198748697e+000λi3:-8.652284007885e+000μ4:-8.146657570993e+000λi4:-8.093482780052e+000μ5:-7.636116393288e+000λi5:-7.659405420574e+000μ6:-7.125575215583e+000λi6:-7.119684646576e+000μ7:-6.615034037878e+000λi7:-6.611764337314e+000μ8:-6.104492860173e+000λi8:-6.066103126985e+000μ9:-5.593951682468e+000λi9:-5.585101045269e+000μ10:-5.083410504763e+000λi10:-5.114083539196e+000μ11:-4.572869327058e+000λi11:-4.578872177367e+000μ12:-4.062328149353e+000λi12:-4.096473385708e+000μ13:-3.551786971648e+000λi13:-3.554211216942e+000μ14:-3.041245793943e+000λi14:-3.041090018133e+000μ15:-2.530704616238e+000λi15:-2.533970334136e+000μ16:-2.020163438533e+000λi16:-2.003230401311e+000μ17:-1.509622260828e+000λi17:-1.503557606947e+000μ18:-9.990810831232e-001λi18:-9.935585987809e-001μ19:-4.885399054182e-001λi19:-4.870426734583e-001μ20:2.200127228676e-002λi20:2.231736249587e-002μ21:5.325424499917e-001λi21:5.324174742068e-001μ22:1.043083627697e+000λi22:1.052898964020e+000μ23:1.553624805402e+000λi23:1.589445977158e+000μ24:2.064165983107e+000λi24:2.060330427561e+000μ25:2.574707160812e+000λi25:2.558075576223e+000μ26:3.085248338516e+000λi26:3.080240508465e+000μ27:3.595789516221e+000λi27:3.613620874136e+000μ28:4.106330693926e+000λi28:4.091378506834e+000μ29:4.616871871631e+000λi29:4.603035354280e+000μ30:5.127413049336e+000λi30:5.132924284378e+000μ31:5.637954227041e+000λi31:5.594906275501e+000μ32:6.148495404746e+000λi32:6.080933498348e+000μ33:6.659036582451e+000λi33:6.680354121496e+000μ34:7.169577760156e+000λi34:7.293878467852e+000μ35:7.680118937861e+000λi35:7.717111851857e+000μ36:8.190660115566e+000λi36:8.225220016407e+000μ37:8.701201293271e+000λi37:8.648665837870e+000μ38:9.211742470976e+000λi38:9.254200347303e+000μ39:9.722283648681e+000λi39:9.724634099672e+000cond(A)=1.925042185755e+003detA=2.772786141752e+118四、实验分析：实验发现，当初始向量u0中0元素的个数较多时，计算结果与准确值相差较大；而当增加初始向量u0中1元素的个数，发现其结果更加靠近准确值。所以，初始向量的选取会对计算结果产生一定的影响。另外，初始向量选择不当时将导致迭代中x1的系数等于零。但是，由于舍入误差的影响，经若干步迭代后，按照基向量展开时，x1的系数可能不等于零，把这一向量看作初始向量，用幂法继续求向量序列，仍然会得出预期的结果，不过收敛速度较慢。如果收敛很慢，可改换初始向量，所以初始向量的选取还会影响到迭代的次数。北航《数值分析》计算实习题目二一．算法设计方案:1.矩阵的QR分解。把矩阵A分解为一个正交矩阵Q与一个上三角矩阵R的乘积，称为矩阵A的正三角分解，简称QR分解。QR分解的算法如下：记，并记，令（n阶单位矩阵）对于r=1,2,…,n-1执行若全为零，则令，转(5);否则转(2)计算(若,则取)令计算继续当此算法执行完后就得到正交矩阵和上三角矩阵且有。2.矩阵的

拟上三角化。对实矩阵A的拟上三角化具体算法如下：记，并记的第r列到第n列的元素为。对于执行若全为零，则令,转(5);否则转(2)。计算 令计算继续算法执行完后，就得到与原矩阵A相似的拟上三角矩阵。3.带双步位移的QR方法求实矩阵全部特征值的具体算法如下：使用矩阵的拟上三角化的算法把矩阵化为拟上三角矩阵；给定精度水平和迭化最大次数L。记，令。如果，则得到A的一个特征值，置（降阶），转（4）；否则转（5）。如果,则得到A的一个特征值，转(11)；如果,则直接转(11);如果转（3）。求二阶子阵的两个特征值和，即计算二次方程的两个根和如果,则得到A的两个特征值和,转(11)；否则转(7)。如果,则得到A的两个特征值和，置(降阶),转(4);否则转(8)。如果k=L,则计算终止，未得到A的全部特征值；否则转9。记,计算置，转（3）。A的全部特征值已计算完毕，停止计算。4.算法的主过程:根据公式生成原始矩阵对矩阵A进行拟上三角化用带双步位移的QR方法求矩阵的特征值对每个特征值λ，解出线性方程组的所有非零解即得A的属于λ的全部特征向量。二．全部源程序：#include<stdio.h>#include<math.h>#include<stdlib.h>#defineL2500#definen11#definee0.000000000001inti,j,s,p,k,ik;doubleA[n][n],q[n][n],r[n][n],rq[n][n],I[n][n];doubleP[n],W[n],u[n],Q[n];doubledr,cr,hr,ar,tr;intnR,nC;doubles1,t,x,tzR[n],tzC[2][n],sum,M[n][n],v[n];voidjuzhenA()//生成矩阵A{for(i=1;i<11;i++){for(j=1;j<11;j++) {if(j!=i)A[i][j]=sin(0.5*i+0.2*j);elseA[i][j]=1.5*cos(i+1.2*j);}}}voidhess(){for(s=1;s<n-2;s++){for(ar=0.0,i=s+2;i<n;i++) ar+=A[i][s]*A[i][s];//printf("ar=%1.12e\n",ar);if(ar==0) continue; else{ ar+=A[s+1][s]*A[s+1][s]; dr=sqrt(ar);//printf("dr=%1.12e\n",dr); if(A[s+1][s]>0)cr=-dr; elsecr=dr;hr=cr*cr-cr*A[s+1][s];for(i=1;i<=s;i++)//生成u向量u[i]=0.0; u[s+1]=A[s+1][s]-cr; for(i=s+2;i<n;i++) u[i]=A[i][s]; for(j=1;j<n;j++)//P {for(P[j]=0.0,i=1;i<n;i++) P[j]+=A[i][j]*u[i]/hr;} for(tr=0.0,i=1;i<n;i++)//tr {tr+=P[i]*u[i]/hr;} for(i=1;i<n;i++)//Q {for(Q[i]=0.0,j=1;j<n;j++) Q[i]+=A[i][j]*u[j]/hr;} for(i=1;i<n;i++)//W {W[i]=Q[i]-tr*u[i];} for(i=1;i<n;i++)//生成Ar+1 {for(j=1;j<n;j++) A[i][j]-=(W[i]*u[j]+u[i]*P[j]);}}}for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(A[i][j])<e)A[i][j]=0.0;}printf("拟上三角化后A(n-1):\n"); for(i=1;i<n;i++) {for(j=1;j<n;j++)printf("%1.12e,",A[i][j]); printf("\n");}}voidAQR()//矩阵A的QR分解{doubleu[n],w[n],F[n]; for(s=1;s<n;s++) for(p=1;p<n;p++) r[s][p]=A[s][p]; for(s=1;s<n-1;s++) { for(dr=0.0,i=s;i<n;i++) dr+=r[i][s]*r[i][s]; dr=sqrt(dr); if(dr==0) continue; else {if(A[s][s]>0)cr=-dr; elsecr=dr; hr=cr*cr-cr*r[s][s]; for(i=1;i<s;i++)//u u[i]=0; u[s]=r[s][s]-cr; for(i=s+1;i<n;i++) u[i]=r[i][s]; for(i=1;i<n;i++)//F {for(F[i]=0.0,j=1;j<n;j++) F[i]+=r[j][i]*u[j]/hr;} for(i=1;i<n;i++)//r {for(j=1;j<n;j++) r[i][j]=r[i][j]-u[i]*F[j];} for(i=1;i<n;i++)//w {for(w[i]=0.0,j=1;j<n;j++) w[i]+=q[i][j]*u[j];} for(i=1;i<n;i++)//Q {for(j=1;j<n;j++) q[i][j]-=w[i]*u[j]/hr;} } } for(i=1;i<n;i++) {for(j=1;j<n;j++) {for(rq[i][j]=0.0,s=1;s<n;s++) rq[i][j]+=r[i][s]*q[s][j];}} printf("生成的Q阵：\n"); for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(q[i][j])<e)q[i][j]=0.0;} for(i=1;i<n;i++) {for(j=1;j<n;j++)printf("%1.12e,",q[i][j]); printf("\n");} printf("生成的R阵：\n"); for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(r[i][j])<e)r[i][j]=0.0;} for(i=1;i<n;i++) {for(j=1;j<n;j++)printf("%1.12e,",r[i][j]); printf("\n");} printf("生成的RQ阵：\n"); for(i=1;i<n;i++){for(j=1;j<n;j++)if(fabs(rq[i][j])<e)rq[i][j]=0.0;} for(i=1;i<n;i++) {for(j=1;j<n;j++)printf("%1.12e,",rq[i][j]); printf("\n");}}voidQRmeth(){ intK=1,m=10; nR=0,nC=0;loop3:if(fabs(A[m][m-1])<=e) {nR++; tzR[nR]=A[m][m]; m--; gotoloop4; } elsegotoloop5;loop4:if(m==1) { nR++; tzR[nR]=A[1][1]; gotoloop9; } elseif(m==2) { s1=A[m-1][m-1]+A[m][m]; t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m]; x=s1*s1-4*t; if(x>=0) { nR++; tzR[nR]=(s1+sqrt(x))/2; nR++; tzR[nR]=(s1-sqrt(x))/2; } else { nC++; tzC[0][nC]=s1/2; tzC[1][nC]=sqrt(-x)/2; nC++; tzC[0][nC]=s1/2; tzC[1][nC]=-sqrt(-x)/2; } gotoloop9; } elsegotoloop3;loop5:{ if(fabs(A[m-1][m-2])<=e) { s1=A[m-1][m-1]+A[m][m]; t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m]; x=s1*s1-4*t; if(x>=0) { nR++; tzR[nR]=(s1+sqrt(x))/2; nR++; tzR[nR]=(s1-sqrt(x))/2; } else { nC++; tzC[0][nC]=s1/2; tzC[1][nC]=sqrt(-x)/2; nC++; tzC[0][nC]=s1/2; tzC[1][nC]=-sqrt(-x)/2; } m--; m--; gotoloop4; } elsegotoloop6; }loop6:{if(K==L) printf("计算终止，未能得到全部特征值\n"); elsegotoloop7; }loop7:{ s1=A[m-1][m-1]+A[m][m]; t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m]; for(i=1;i<=m;i++)//生成M for(j=1;j<=m;j++) { for(sum=0.0,s=1;s<=m;s++) sum+=A[i][s]*A[s][j]; M[i][j]=sum-s1*A[i][j]+t*I[i][j]; } for(s=1;s<=m;s++) { for(ar=0.0,i=s+1;i<=m;i++) ar+=M[i][s]*M[i][s]; if(ar==0) continue; else { ar+=M[s][s]*M[s][s]; dr=sqrt(ar); if(M[s][s]>0) cr=-dr; elsecr=dr; hr=cr*cr-cr*M[s][s]; for(i=1;i<s;i++) u[i]=0.0; u[s]=M[s][s]-cr; for(i=s+1;i<=m;i++) u[i]=M[i][s]; for(j=1;j<=m;j++)//v for(v[j]=0.0,i=1;i<=m;i++) v[j]+=M[i][j]*u[i]/hr; for(i=1;i<=m;i++)//Mr+1 for(j=1;j<=m;j++) M[i][j]-=u[i]*v[j]; for(j=1;j<=m;j++)//P for(P[j]=0.0,i=1;i<=m;i++) P[j]+=A[i][j]*u[i]/hr; for(i=1;i<=m;i++)//Q for(Q[i]=0.0,j=1;j<=m;j++) Q[i]+=A[i][j]*u[j]/hr; for(tr=0.0,i=1;i<=m;i++)//tr tr+=P[i]*u[i]/hr; for(i=1;i<=m;i++)//W W[i]=Q[i]-tr*u[i]; for(i=1;i<n;i++)//Ar+1 for(j=1;j<n;j++) A[i][j]-=(W[i]*u[j]+u[i]*P[j]); } } gotoloop8; }loop8:{k++;gotoloop3;}loop9:{ ; } printf("矩阵的全部特征值为：\n"); for(i=1;i<=nR;i++) printf("%1.12e\n",tzR[i]); for(i=1;i<=nC;i++) { printf("%1.12e+0*i",tzC[0][i]); if(tzC[1][i]>=0) printf("+*i%1.12e\n",tzC[1][i]); elseprintf("*i%1.12e\n",tzC[1][i]); } }voidgauss()//列主元的高斯消元法求解特征向量{ doublech,m[n],x[n][n]; for(p=1;p<=nR;p++) { juzhenA(); for(i=1;i<n;i++) for(j=1;j<n;j++) A[i][j]-=tzR[p]*I[i][j]; for(k=1;k<n;k++) { for(ik=k,i=k;i<n;i++) if(A[ik][k]<A[i][k]) ik=i; for(j=k;j<n;j++) { ch=A[ik][j]; A[ik][j]=A[k][j]; A[k][j]=ch; } for(i=k+1;i<n;i++) { m[i]=A[i][k]/A[k][k]; for(j=k+1;j<n;j++) A[i][j]-=m[i]*A[k][j]; } } x[p][n-1]=1.0; for(k=n-2;k>0;k--) { for(ar=0.0,j=k+1;j<=n;j++) ar+=A[k][j]*x[p][j]; } } for(p=1;p<=nR;p++) { printf("对应特征值%1.12e的特征向量为:\n",tzR[p]); for(i=1;i<n;i++) printf("%1.12e\n",x[p][i]); }}voidmain(){for(i=1;i<n;i++) {for(j=1;j<n;j++) {if(i==j) {I[i][j]=1;q[i][j]=1;} else{I[i][j]=0;q[i][j]=0;}}}juzhenA();hess();AQR();QRmeth();gauss();}三．实验运行结果：（1）拟上三角化后的矩阵A(n-1):-8.827516758830e-001,-9.933136491826e-002,-1.103349285994e+000,-7.600443585637e-001,1.549101079914e-001,-1.946591862872e+000,-8.782436382927e-002,-9.255889387184e-001,6.032599440534e-001,1.518860956469e-001,-2.347878362416e+000,2.372370104937e+000,1.819290822208e+000,3.237804101546e-001,2.205798440320e-001,2.102692662546e+000,1.816138086098e-001,1.278839089990e+000,-6.380578124405e-001,-4.154075603804e-001,0.000000000000e+000,1.728274599968e+000,-1.171467642785e+000,-1.243839262699e+000,-6.399758341743e-001,-2.002833079037e+000,2.924947206124e-001,-6.412830068395e-001,9.783997621284e-002,2.557763574160e-001,0.000000000000e+000,0.000000000000e+000,-1.291669534130e+000,-1.111603513396e+000,1.171346824096e+000,-1.307356030021e+000,1.803699177750e-001,-4.246385358369e-001,7.988955239304e-002,1.608819928069e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,1.560126298527e+000,8.125049397515e-001,4.421756832923e-001,-3.588616128137e-002,4.691742313671e-001,-2.736595050092e-001,-7.359334657750e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-7.707773755194e-001,-1.583051425742e+000,-3.042843176799e-001,2.528712446035e-001,-6.709925401449e-001,2.544619929082e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-7.463453456938e-001,-2.708365157019e-002,-9.486521893682e-001,1.195871081495e-001,1.929265617952e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-7.701801374364e-001,-4.697623990618e-001,4.988259468008e-001,1.137691603776e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,7.013167092107e-001,1.582180688475e-001,3.862594614233e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,4.843807602783e-001,3.992777995177e-001,（2）生成的Q阵：-3.519262579534e-001,4.427591982245e-001,-6.955982513607e-001,6.486200753651e-002,3.709718861896e-001,1.855847143605e-001,1.628942319628e-002,-1.181053169648e-001,-5.255375383724e-002,5.486582943568e-002,-9.360277287361e-001,-1.664679186545e-001,2.615299548560e-001,-2.438671728934e-002,-1.394774360893e-001,-6.977585391241e-002,-6.124472142964e-003,4.440505443139e-002,1.975907909728e-002,-2.062836970533e-002,0.000000000000e+000,-8.810520554692e-001,-3.989762796959e-001,3.720308728479e-002,2.127794064090e-001,1.064463557221e-001,9.343171079758e-003,-6.774200464527e-002,-3.014340698675e-002,3.146955080444e-002,0.000000000000e+000,0.000000000000e+000,-5.371806806439e-001,-1.234945854205e-001,-7.063151608719e-001,-3.533456368496e-001,-3.101438948264e-002,2.248676491598e-001,1.000601783527e-001,-1.044622748702e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,9.892235468621e-001,-1.239411731211e-001,-6.200358589825e-002,-5.442272839461e-003,3.945881637235e-002,1.755813350011e-002,-1.833059462907e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,5.323610690264e-001,-6.733900344896e-001,-5.910581205868e-002,4.285425323867e-001,1.906901343193e-001,-1.990794495295e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-6.059761505747e-001,9.165783032819e-002,-6.645586508974e-001,-2.957110877580e-001,3.087207462557e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-9.933396625117e-001,-9.690440311939e-002,-4.311990584470e-002,4.501694411183e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,5.410088006061e-001,-5.817540838226e-001,6.073480580545e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-7.221591336735e-001,-6.917269588876e-001,生成的R阵：2.508342744917e+000,-2.185646885493e+000,-1.314609070786e+000,-3.558787493836e-002,-2.609857850388e-001,-1.283121847090e+000,-1.390878610606e-001,-8.712897972161e-001,3.849367902971e-001,3.353802899665e-001,0.000000000000e+000,-1.961603277854e+000,2.407523727633e-001,7.054714572823e-001,5.957204318279e-001,5.526978774676e-001,-3.268209924413e-001,-5.769498668365e-002,2.871129330189e-001,-8.895128754189e-002,0.000000000000e+000,0.000000000000e+000,2.404534601993e+000,1.706758096328e+000,-4.239566704091e-001,3.405332305815e+000,-1.050017655853e-001,1.462257102734e+000,-6.684487469283e-001,-4.027646209664e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,1.577122080722e+000,6.399535133956e-001,3.468127872427e-001,-5.701786649768e-002,4.014788054433e-001,-2.222476176311e-001,-6.317059236442e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-1.447846997770e+000,-1.415724007744e+000,-2.806139044665e-001,-2.817910521892e-001,-4.611434881851e-002,1.996629079956e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,1.231641451542e+000,1.619701003419e-001,1.962638275504e-001,5.350035621760e-001,-1.509273424767e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,7.753441914209e-001,3.464514508820e-001,-4.312226803504e-001,-1.234643696237e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,1.296312940612e+000,-4.288053318338e-001,2.737334158165e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-6.707396440648e-001,-4.842320121884e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-7.168323926323e-002,生成的RQ阵：1.163074414164e+000,2.632670934508e+000,-1.772796003272e+000,-8.668899138521e-002,3.300503471047e-001,1.455162371214e+000,9.730650448593e-001,-4.873031174655e-001,-7.756411630489e-001,3.249201979113e-001,1.836115060851e+000,1.144286420080e-001,-9.880381403133e-001,5.589725694767e-001,4.694190067101e-002,-2.978478237007e-001,-1.617130577649e-002,6.936977702522e-001,1.367670571405e-001,-1.419099231519e-002,0.000000000000e+000,-2.118520153533e+000,-1.876189745783e+000,-5.407071940597e-001,1.171538359721e+000,-2.550323020223e+000,-1.691577936540e+000,1.229951613262e+000,1.387947777212e+000,-8.667502917242e-001,0.000000000000e+000,0.000000000000e+000,-8.471995127808e-001,4.382910468318e-001,-1.008632199185e+000,-7.959374261495e-001,-4.769258865577e-001,4.072683083890e-001,4.096390493527e-001,-3.363378940862e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-1.432244342447e+000,-5.742284908055e-001,1.213151477723e+000,3.457508625575e-001,-4.749853573124e-001,-3.176158274191e-001,4.294507015032e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,6.556779598004e-001,-9.275250974463e-001,-2.529079844054e-001,6.905949216976e-001,-2.374430675823e-002,2.429781119781e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-4.698400884876e-001,-2.730776009527e-001,-7.821296259798e-001,9.580964936399e-002,7.846239841323e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-1.287679058937e+000,-3.576058900348e-001,-4.116725408808e-003,-3.914268216423e-001,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,-3.628760503545e-001,7.398980975354e-001,-7.241608309576e-002,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,0.000000000000e+000,5.176670596524e-002,4.958522909877e-002,（3）矩阵A的全部特征值为：6.360627875746e-001+0*i5.650488993501e-002+0*i9.355889078188e-001+0*i1.577548557113e+000+0*i-1.484039822259e+000+0*i3.383039617436e+000+0*i-9.805309562902e-001+1.139489127430e-001*i-9.805309562902e-001-1.139489127430e-001*i-2.323496210212e+000+8.930405177200e-001*i-2.323496210212e+000-8.930405177200e-001*i（4）相应于实特征值的特征向量：对应特征值6.360627875746e-001的特征向量为:-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+0611.000000000000e+000对应特征值5.650488993501e-002的特征向量为:-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+0611.000000000000e+000对应特征值9.355889078188e-001的特征向量为:-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+0611.000000000000e+000对应特征值1.577548557113e+000的特征向量为:-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+0611.000000000000e+000对应特征值-1.484039822259e+000的特征向量为:-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+0611.000000000000e+000对应特征值3.383039617436e+000的特征向量为:-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+061-9.255963134932e+0611.000000000000e+000北航《数值分析》计算实习题目三算法设计方案:先利用二维数表，通过分片二次插值，得到。然后再由数列和，用牛顿法迭代求出非线性方程组，得到共组值（这些均在二维数表的范围内）。由插值函数，求得关于、的函数。曲面拟合的基底采用幂函数。系数矩阵的形式为。这其中有矩阵求逆的运算。求逆可通过解线性方程组的方法实现。，求n个线性方程组，得到即。具体来说是：（1）使用牛顿迭代法(简单迭代法不收敛)，对原题中给出的，，()的11*21组分别求出原题中方程组的一组解，于是得到一组和对应的。（2）对于已求出的，使用分片二次代数插值法对原题中关于的数表进行插值得到。于是产生了z=f(x,y)的11*21个数值解。（3）从k=1开始逐渐增大k的值，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，得到每次的。当时结束计算，输出拟合结果。（4）计算的值并输出结果，以观察逼近的效果。其中。编译环境：Windows7系统，VC++6.02.全部源程序:#defineN6//矩阵的阶；#defineM4//方程组未知数个数；#defineEPSL1.0e-12//迭代的精度水平；#defineMAXDIGIT1.0e-219#defineSIGSUM1.0e-7#defineL500//迭代最大次数；#defineK10//最小二乘法拟合时的次数上限；#defineX_NUM11#defineY_NUM21#defineOUTPUTMODE1//输出格式：0--输出至屏幕，1--输出至文件；#include<stdio.h>#include<math.h>doublemat_u[N]={0,0.4,0.8,1.2,1.6,2.0},mat_t[N]={0,0.2,0.4,0.6,0.8,1.0},mat_ztu[N][N]={{-0.5,-0.34,0.14,0.94,2.06,3.5},{-0.42,-0.5,-0.26,0.3,1.18,2.38},{-0.18,-0.5,-0.5,-0.18,0.46,1.42},{0.22,-0.34,-0.58,-0.5,-0.1,0.62},{0.78,-0.02,-0.5,-0.66,-0.5,-0.02},{1.5,0.46,-0.26,-0.66,-0.74,-0.5},},mat_val_z[X_NUM][Y_NUM]={0},init_tuvw[4]={1,2,1,2},mat_crs[K][K]={0};FILE*p;intmain()//主程序入口{inti,j,x,y,kmin,tmp=0;doubletmpval;intgetval_z(double,double,int,int);intleasquare();voidresult_out(int);if(OUTPUTMODE){p=fopen("c:\\Result.txt","w+");fprintf(p,"计算结果：\n");fclose(p);}for(i=0;i<X_NUM;i++)for(j=0;j<Y_NUM;j++)tmp+=getval_z(0.08*i,0.5+0.05*j,i,j);if(!tmp)printf("迭代求解z=f(x,y)完毕。\n");elseprintf("迭代超过最大次数，计算结果可能不准确！\n");if(kmin=leasquare()){printf("近似表达式已求出！\n");for(i=0;i<8;i++)for(j=0;j<5;j++)tmp+=getval_z(0.1*(i+1),0.5+0.2*(j+1),i,j);if(!tmp){printf("迭代求解z=f(x*,y*)完毕。\n");for(x=0;x<8;x++){for(y=0;y<5;y++){tmpval=0;for(i=0;i<kmin;i++){for(j=0;j<kmin;j++){tmpval=tmpval+mat_crs[i][j]*pow(0.1*(x+1),i)*pow(0.5+0.2*(y+1),j);}}if(OUTPUTMODE){p=fopen("c:\\Result.txt","at+");fprintf(p,"x*[%d]=%f,y*[%d]=%f:f(x*[%d],