北航數(shù)值分析全部三次大作業(yè)

北航《數(shù)值分析》計算實習題目一一、算法設(shè)計：要計算矩陣的最大最小特征值，可先通過冪法求得模最大的特征值λa，再根據(jù)λa進行原點平移求出下一特征值λb，比較兩值的大小，循環(huán)下去，最終可以確定最小特征值λ1和最大特征值λ501。要計算矩陣按模最小的特征值λs，可以用反冪法求得，此處需要求解中間過渡向量，采用LU分解法解帶狀方程組，故要設(shè)計Doolite分解。要求解與給定數(shù)值接近的特征值，可以對該數(shù)做漂移量，新數(shù)組特征值倒數(shù)的絕對值滿足反冪法的要求，故通過反冪法即可求得。要計算cond(A)，可由公式cond(A)=│λmax/λmin│求得，其中λmax和λmin分別矩陣A的模為最大和模為最小的特征值，λmax可以通過冪法求得，而λmin可已通過反冪法求得。要計算detA，可對A作LU分解，即A=LU，則│A│=│L││U│，其中L為單位下三角陣，U為上三角陣，則detA即為U矩陣主對角元素之積。由于矩陣A中的零元素較多，為了節(jié)省存儲空間，可以將矩陣A中的元素存儲在5*501的二維數(shù)組中。編譯環(huán)境：vc++6.0編譯函數(shù)：冪法，反冪法，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個判斷大小的函數(shù)*/intMax(intv1,intv2){ return((v1>v2)?v1:v2);}intMin(intv1,intv2){ return((v1<v2)?v1:v2);}/*將矩陣值轉(zhuǎn)存在一個數(shù)組里*/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));//設(shè)置計算精度 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;}三、運行結(jié)果:λ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四、實驗分析：實驗發(fā)現(xiàn)，當初始向量u0中0元素的個數(shù)較多時，計算結(jié)果與準確值相差較大；而當增加初始向量u0中1元素的個數(shù)，發(fā)現(xiàn)其結(jié)果更加靠近準確值。所以，初始向量的選取會對計算結(jié)果產(chǎn)生一定的影響。另外，初始向量選擇不當時將導致迭代中x1的系數(shù)等于零。但是，由于舍入誤差的影響，經(jīng)若干步迭代后，按照基向量展開時，x1的系數(shù)可能不等于零，把這一向量看作初始向量，用冪法繼續(xù)求向量序列，仍然會得出預期的結(jié)果，不過收斂速度較慢。如果收斂很慢，可改換初始向量，所以初始向量的選取還會影響到迭代的次數(shù)。北航《數(shù)值分析》計算實習題目二一．算法設(shè)計方案:1.矩陣的QR分解。把矩陣A分解為一個正交矩陣Q與一個上三角矩陣R的乘積，稱為矩陣A的正三角分解，簡稱QR分解。QR分解的算法如下：記，并記，令（n階單位矩陣）對于r=1,2,…,n-1執(zhí)行若全為零，則令，轉(zhuǎn)(5);否則轉(zhuǎn)(2)計算(若,則取)令計算繼續(xù)當此算法執(zhí)行完后就得到正交矩陣和上三角矩陣且有。2.矩陣的

擬上三角化。對實矩陣A的擬上三角化具體算法如下：記，并記的第r列到第n列的元素為。對于執(zhí)行若全為零，則令,轉(zhuǎn)(5);否則轉(zhuǎn)(2)。計算 令計算繼續(xù)算法執(zhí)行完后，就得到與原矩陣A相似的擬上三角矩陣。3.帶雙步位移的QR方法求實矩陣全部特征值的具體算法如下：使用矩陣的擬上三角化的算法把矩陣化為擬上三角矩陣；給定精度水平和迭化最大次數(shù)L。記，令。如果，則得到A的一個特征值，置（降階），轉(zhuǎn)（4）；否則轉(zhuǎn)（5）。如果,則得到A的一個特征值，轉(zhuǎn)(11)；如果,則直接轉(zhuǎn)(11);如果轉(zhuǎn)（3）。求二階子陣的兩個特征值和，即計算二次方程的兩個根和如果,則得到A的兩個特征值和,轉(zhuǎn)(11)；否則轉(zhuǎn)(7)。如果,則得到A的兩個特征值和，置(降階),轉(zhuǎn)(4);否則轉(zhuǎn)(8)。如果k=L,則計算終止，未得到A的全部特征值；否則轉(zhuǎn)9。記,計算置，轉(zhuǎn)（3）。A的全部特征值已計算完畢，停止計算。4.算法的主過程:根據(jù)公式生成原始矩陣對矩陣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();}三．實驗運行結(jié)果：（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北航《數(shù)值分析》計算實習題目三算法設(shè)計方案:先利用二維數(shù)表，通過分片二次插值，得到。然后再由數(shù)列和，用牛頓法迭代求出非線性方程組，得到共組值（這些均在二維數(shù)表的范圍內(nèi)）。由插值函數(shù)，求得關(guān)于、的函數(shù)。曲面擬合的基底采用冪函數(shù)。系數(shù)矩陣的形式為。這其中有矩陣求逆的運算。求逆可通過解線性方程組的方法實現(xiàn)。，求n個線性方程組，得到即。具體來說是：（1）使用牛頓迭代法(簡單迭代法不收斂)，對原題中給出的，，()的11*21組分別求出原題中方程組的一組解，于是得到一組和對應的。（2）對于已求出的，使用分片二次代數(shù)插值法對原題中關(guān)于的數(shù)表進行插值得到。于是產(chǎn)生了z=f(x,y)的11*21個數(shù)值解。（3）從k=1開始逐漸增大k的值，并使用最小二乘法曲面擬合法對z=f(x,y)進行擬合，得到每次的。當時結(jié)束計算，輸出擬合結(jié)果。（4）計算的值并輸出結(jié)果，以觀察逼近的效果。其中。編譯環(huán)境：Windows7系統(tǒng)，VC++6.02.全部源程序:#defineN6//矩陣的階；#defineM4//方程組未知數(shù)個數(shù)；#defineEPSL1.0e-12//迭代的精度水平；#defineMAXDIGIT1.0e-219#defineSIGSUM1.0e-7#defineL500//迭代最大次數(shù)；#defineK10//最小二乘法擬合時的次數(shù)上限；#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,"計算結(jié)果：\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("迭代超過最大次數(shù)，計算結(jié)果可能不準確！\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],