貴州大學(xué)2013級數(shù)值分析實(shí)驗報告

本次作業(yè)共七道題（第一道大題有三個小題，按三道題算），程序加載了StanfordCPPLib.lib庫，主要利用“grid.h”類似二維數(shù)組，“vector.h”類似一維數(shù)組，目的是使程序更通用，安全，方便編程。1.householder變換：程序：#include<iostream>#include<fstream>#include<cmath>#include<iomanip>#include"vector.h"#include"grid.h"#include"foreach.h"usingnamespacestd;boolreadFile(Grid<double>&Matrix);//從文件中讀入數(shù)據(jù)voidhouseholderFabrication(Grid<double>&Matrix);//對文件中的二維數(shù)組進(jìn)行householder變換boolcomputeW(Grid<double>&Matrix,intn,Vector<double>&temp);voiddoElimination(Grid<double>&Matrix,intn,Vector<double>&temp);doublexMultiplyW(Vector<double>&temp,Grid<double>&Matrix,introw,intcol);doublewTRansMultiplyX(Vector<double>&temp,Grid<double>&Matrix,introw,intcol);doublenorm(Vector<double>&temp);//計算向量的長度的平方squre(||A||)voidwriteFile(Grid<double>&Matrix);//變換后的三對角陣寫入文件，為以后使用intmain(){ Grid<double>Matrix; readFile(Matrix); householderFabrication(Matrix); writeFile(Matrix); return0;}boolreadFile(Grid<double>&Matrix){ intnumRows,numCols; stringfilename="abc.txt"; ifstreaminfile; infile.open(filename.c_str()); infile>>ws>>numRows>>ws>>numCols; Matrix.resize(numRows,numCols); for(inti=0;i<numRows;i++){ for(intj=0;j<numCols;j++){ infile>>ws>>Matrix[i][j]; } } infile.close(); returntrue;}voidhouseholderFabrication(Grid<double>&Matrix){//對矩陣A（非增廣矩陣）進(jìn)行變換 //Grid<double>Householder(Matrix.numRows(),Matrix.numCols()); Vector<double>temp; for(inti=0;i<Matrix.numRows()-2;i++){ boolelimination=computeW(Matrix,i,temp);//用i標(biāo)記消元的起始列，當(dāng)i=0時，對子矩陣Ai進(jìn)行消元，temp返回消元向量 if(elimination==true){ doElimination(Matrix,i,temp);//用i標(biāo)記消元的起始列，當(dāng)i=0時，對子矩陣Ai進(jìn)行消元， } }}boolcomputeW(Grid<double>&Matrix,intn,Vector<double>&temp){ temp.clear();//初始化temp for(inti=n+1;i<Matrix.numRows();i++){ temp.add(Matrix[i][n]); } doublelength=sqrt(norm(temp)); if(length==0)returnfalse;//若temp是零向量則不做任何操作，進(jìn)入下一個循環(huán)，否則會出錯 if(temp[0]<0){ temp[0]=length-temp[0]; }else{ temp[0]=-length-temp[0]; } //length=norm(temp); //temp[0]/=length; for(inti=1;i<temp.size();i++){ temp[i]=-temp[i]; } returntrue; }voiddoElimination(Grid<double>&Matrix,intn,Vector<double>&temp){ for(inti=n;i<Matrix.numCols();i++){ doubleintermediateValue=wTRansMultiplyX(temp,Matrix,n+1,i);//用n，i標(biāo)記HwX=X-2w（w'X）中的X,用i（行）控制，第一次消元即從第二行第一列開始； intermediateValue/=norm(temp); intk=0; for(intj=n+1;j<Matrix.numRows();j++){ Matrix[j][i]=Matrix[j][i]-2*temp[k++]*intermediateValue; } } for(inti=n;i<Matrix.numRows();i++){ doubleintermediateValue=xMultiplyW(temp,Matrix,i,n+1); intermediateValue/=norm(temp); intk=0; for(intj=n+1;j<Matrix.numCols();j++){ Matrix[i][j]=Matrix[i][j]-2*temp[k++]*intermediateValue; } }}doublexMultiplyW(Vector<double>&temp,Grid<double>&Matrix,introw,intcol){ intk=0; doubleintermidiate=0; for(inti=col;i<Matrix.numCols();i++){ intermidiate+=Matrix[row][i]*temp[k++]; } returnintermidiate;}doublewTRansMultiplyX(Vector<double>&temp,Grid<double>&Matrix,introw,intcol){ intk=0; doubleintermidiate=0; for(inti=row;i<Matrix.numRows();i++){ intermidiate+=Matrix[i][col]*temp[k++]; } returnintermidiate;}doublenorm(Vector<double>&temp){ doublenorm=0; foreach(doublenumintemp){ norm+=num*num; } returnnorm;}voidwriteFile(Grid<double>&Matrix){ ofstreamoutfile; outfile.open("123.txt"); for(inti=0;i<Matrix.numRows();i++){ for(intj=0;j<Matrix.numCols();j++){ outfile<<right<<fixed<<setprecision(7)<<setw(12)<<Matrix[i][j]<<""; } outfile<<'\n'; } outfile.close();}輸入數(shù)據(jù)（從文件）：計算結(jié)果（文件顯示）：算法分析：輸入矩陣為A，householder矩陣為Hw=I-2*w*w’,w=(v-u)/(||v-u||2),Hw*u=v,A為9*9矩陣，第一次Hw*A*Hw后，計算Am，Am為8*8矩陣，一次下去，實(shí)際并不生成Hw，生成矩陣既耗費(fèi)內(nèi)存，有需要多次的乘法操作，而是Hw*a=（I-2*w*w’）*x=I*x-2*w*w’*x=x-2*w*(w’*x),使用向量w就可以計算Hw*x，不需要形成大型矩陣，減少運(yùn)算次數(shù)，x為A的列向量，Hw左乘A后，對A右乘。2．超松弛法：程序：#include<iostream>#include<fstream>#include<iomanip>#include"grid.h"#include"vector.h"#include"foreach.h"usingnamespacestd;constdoubleSOR_W=1.4;constintITERATION_TIMES=9;boolreadFile(Grid<double>&Matrix);voidSOR_Iteration(Grid<double>&Matrix,Vector<double>&Solution,Vector<double>&SolutionNew);voidreplaceOldSolution(Vector<double>&Solution,Vector<double>&SolutionNew);voidprintSolution(Vector<double>&Solution);intmain(){ Grid<double>Matrix; readFile(Matrix); Vector<double>Solution(Matrix.numRows(),0); Vector<double>SolutionNew(Matrix.numRows(),0); SOR_Iteration(Matrix,Solution,SolutionNew); printSolution(Solution); return0;}boolreadFile(Grid<double>&Matrix){ intnumRows,numCols; stringfilename="abc.txt"; ifstreaminfile; infile.open(filename.c_str()); infile>>ws>>numRows>>ws>>numCols; Matrix.resize(numRows,numCols); for(inti=0;i<numRows;i++){ for(intj=0;j<numCols;j++){ infile>>ws>>Matrix[i][j]; } } returntrue;}voidSOR_Iteration(Grid<double>&Matrix,Vector<double>&Solution,Vector<double>&SolutionNew){ for(inti=0;i<ITERATION_TIMES;i++){ for(intj=0;j<Matrix.numRows();j++){ doublesum=0; for(intk=0;k<j;k++){ sum+=Matrix[j][k]*SolutionNew[k]; } for(intk=j+1;k<Matrix.numRows();k++){ sum+=Matrix[j][k]*Solution[k]; } SolutionNew[j]=(1-SOR_W)*Solution[j]-SOR_W/Matrix[j][j]*(sum-Matrix[j][Matrix.numCols()-1]); } replaceOldSolution(Solution,SolutionNew); }}voidreplaceOldSolution(Vector<double>&Solution,Vector<double>&SolutionNew){ for(inti=0;i<Solution.size();i++){ Solution[i]=SolutionNew[i]; }}voidprintSolution(Vector<double>&Solution){ cout<<"Thesolution:"; foreach(doublexinSolution){ cout<<fixed<<setprecision(10)<<x<<""; } cout<<'\n';}輸入（從文件）即householder程序的計算結(jié)果（123.txt）輸出結(jié)果：SOR迭代適合大型稀疏矩陣，在情況好時，矩陣收斂，則計算速度快，迭代速度由迭代矩陣的最大特征值決定。3.列主元素消去法：#include<iostream>#include<fstream>#include<iomanip>#include"grid.h"#include"vector.h"#include"foreach.h"#include<cmath>#include"error.h"usingnamespacestd;boolreadFile(Grid<double>&Matrix);boolelimination(Grid<double>&Matrix);voidexchangeRows(Grid<double>&Matrix,introw);voidbackSubstitution(Grid<double>&Matrix,Vector<double>&Solution);voidprintSolution(Vector<double>&Solution);intmain(){ Grid<double>Matrix; readFile(Matrix); Vector<double>Solution(Matrix.numRows()); elimination(Matrix); backSubstitution(Matrix,Solution); printSolution(Solution); return0;}boolreadFile(Grid<double>&Matrix){ intnumRows,numCols; stringfilename="abc.txt"; ifstreaminfile; infile.open(filename.c_str()); infile>>ws>>numRows>>ws>>numCols; Matrix.resize(numRows,numCols); for(inti=0;i<numRows;i++){ for(intj=0;j<numCols;j++){ infile>>ws>>Matrix[i][j]; } } returntrue;}boolelimination(Grid<double>&Matrix){ for(inti=0;i<Matrix.numRows()-1;i++){ exchangeRows(Matrix,i)； for(intj=i+1;j<Matrix.numRows();j++){ doubleci=Matrix[j][i]/Matrix[i][i]; for(intk=i;k<Matrix.numCols();k++){ Matrix[j][k]=Matrix[j][k]-Matrix[i][k]*ci; } } for(intl=0;l<Matrix.numRows();l++){ for(intm=0;m<Matrix.numCols();m++){ cout<<setw(10)<<Matrix[l][m]<<""; } cout<<'\n'; } cout<<'\n'; } returntrue;}voidexchangeRows(Grid<double>&Matrix,introw){ doublemax=fabs(Matrix[row][row]); intlabel=row; doubletemp; for(inti=row;i<Matrix.numRows();i++){ if(fabs(Matrix[i][row])>max){ max=fabs(Matrix[i][row]); label=i; } } if(max==0)error("Thereisnosolution!Thematrixisnotinvertible"); if(label!=row){ for(intj=row;j<Matrix.numCols();j++){ temp=Matrix[label][j]; Matrix[label][j]=Matrix[row][j]; Matrix[row][j]=temp; } }}voidbackSubstitution(Grid<double>&Matrix,Vector<double>&Solution){ intkNumRows=Matrix.numRows(); intkNumCols=Matrix.numCols(); intk=0; for(inti=kNumRows-1;i>=0;i--){ doubletemp=0; k++; for(intj=kNumCols-2;j>kNumCols-k-1;j--){ temp+=Matrix[i][j]*Solution[j]; } temp=Matrix[i][kNumCols-1]-temp; Solution[i]=temp/Matrix[i][i]; }}voidprintSolution(Vector<double>&Solution){ cout<<"Thesolution:"; foreach(doublexinSolution){ cout<<setprecision(10)<<x<<""; } cout<<'\n';}輸入（從文件）即householder程序的計算結(jié)果（123.txt）輸出結(jié)果：高斯列主元消去法：高斯消去法是一個自然的想法，通過選主元，可以減少舍入誤差，對于1000列以下的良好方陣是可以的，對于大型矩陣，舍入誤差是不穩(wěn)定的，若矩陣可以分塊，則可以先對其分塊進(jìn)行運(yùn)算。程序的不足：沒有使用指針，效率不高，選主元時使用指針，可避免數(shù)據(jù)在內(nèi)存中的移動，對于大型矩陣，可顯著提高效率。4．固定樣條插值 程序：#include<iostream>#include<fstream>#include<iomanip>usingnamespacestd;#include"vector.h"constdoubleSENTINEL=-1;voidreadFile(Vector<double>&XNodes,Vector<double>&YNodes,double&S1,double&Sn);voidspanLinearSystem(Vector<double>&XNodes,Vector<double>&YNodes,double&S1,double&Sn,Vector<double>&Hi,Vector<double>&upper,Vector<double>&diag,Vector<double>&lower,Vector<double>&b);voidsolveTridiagnalMatrix(Vector<double>&upper,Vector<double>&diag,Vector<double>&lower,Vector<double>&b,Vector<double>&Solution);voidcalculatePP(Vector<double>&Solution,Vector<double>&Hi,Vector<double>&YNodes,Vector<double>&pp);voidcalculateInterpolateValue(Vector<double>&pp,Vector<double>&XNodes,Vector<double>&Hi);intmain(){ Vector<double>XNodes;//x節(jié)點(diǎn) Vector<double>YNodes;//y節(jié)點(diǎn) Vector<double>Hi;//distanceofXNodes Vector<double>upper;//生成系數(shù)矩陣對角線上部元素 Vector<double>lower;//系數(shù)矩陣對角線下部元素 Vector<double>b;//Ax=b Vector<double>pp;//樣條系數(shù)向量，p1（x）=ax^3+bx^2+cx+d;則pp=[abcd] doubleS1,Sn; readFile(XNodes,YNodes,S1,Sn); Vector<double>diag(XNodes.size(),2);//系數(shù)矩陣對角線元素 spanLinearSystem(XNodes,YNodes,S1,Sn,Hi,upper,diag,lower,b);//生成樣條曲線的方程，并按序?qū)懭胛募? Vector<double>Solution(diag.size());//Ax=b的解 solveTridiagnalMatrix(upper,diag,lower,b,Solution); calculatePP(Solution,Hi,YNodes,pp); calculateInterpolateValue(pp,XNodes,Hi); return0;}voidreadFile(Vector<double>&XNodes,Vector<double>&YNodes,double&S1,double&Sn){ ifstreaminfile; infile.open("InterpolationNodes.txt"); doublenodes; infile>>nodes; doublenumX,numY; for(inti=0;i<nodes;i++){ infile>>numX>>numY; XNodes.add(numX); YNodes.add(numY); } infile>>S1>>Sn; infile.close(); }voidspanLinearSystem(Vector<double>&XNodes,Vector<double>&YNodes,double&S1,double&Sn,Vector<double>&Hi,Vector<double>&upper,Vector<double>&diag,Vector<double>&lower,Vector<double>&b){ intN=XNodes.size(); for(inti=0;i<N-1;i++){ Hi.add(XNodes[i+1]-XNodes[i]); } upper.add(1); for(inti=0;i<N-2;i++){ lower.add(Hi[i]/(Hi[i]+Hi[i+1])); upper.add(1-lower[i]); } lower.add(1); b.add(6*(YNodes[1]-YNodes[0])/(Hi[0]*Hi[0])-6*S1/Hi[0]); for(inti=1;i<N-1;i++){ b.add(6.0/(Hi[i]+Hi[i-1])*((YNodes[i+1]-YNodes[i])/Hi[i]-(YNodes[i]-YNodes[i-1])/Hi[i-1])); } b.add(-6*(YNodes[N-1]-YNodes[N-2])/(Hi[N-2]*Hi[N-2])+6*Sn/Hi[N-2]); /*for(inti=0;i<upper.size();i++){ cout<<upper[i]<<""; } cout<<endl; for(inti=0;i<diag.size();i++){ cout<<diag[i]<<""; } cout<<endl; for(inti=0;i<lower.size();i++){ cout<<lower[i]<<""; } cout<<endl; for(inti=0;i<b.size();i++){ cout<<b[i]<<""; } cout<<endl;*/}voidsolveTridiagnalMatrix(Vector<double>&upper,Vector<double>&diag,Vector<double>&lower,Vector<double>&b,Vector<double>&Solution){ intN=diag.size(); for(inti=0;i<N-1;i++){ lower[i]=lower[i]/diag[i]; diag[i+1]=diag[i+1]-lower[i]*upper[i]; b[i+1]=b[i+1]-b[i]*lower[i]; } Solution[N-1]=b[N-1]/diag[N-1]; for(inti=N-2;i>=0;i--){ Solution[i]=(b[i]-Solution[i+1]*upper[i])/diag[i]; } /*for(inti=0;i<Solution.size();i++){ cout<<Solution[i]<<""; } cout<<endl;*/ }voidcalculatePP(Vector<double>&Solution,Vector<double>&Hi,Vector<double>&YNodes,Vector<double>&pp){ intN=YNodes.size()-1; for(inti=0;i<N;i++){ pp.add((Solution[i+1]-Solution[i])/(6*Hi[i])); pp.add(Solution[i]/2); pp.add((YNodes[i+1]-YNodes[i])/Hi[i]-Hi[i]*(Solution[i+1]+2*Solution[i])/6); pp.add(YNodes[i]); } /*for(inti=0;i<N;i++){ cout<<setprecision(12)<<fixed<<right<<pp[i]<<""; cout<<pp[i+1]<<""; cout<<pp[i+2]<<""; cout<<pp[i+3]<<""; cout<<endl; } cout<<endl;*/ }voidcalculateInterpolateValue(Vector<double>&pp,Vector<double>&XNodes,Vector<double>&Hi){ cout<<"whenyouenter-1,theprogrammewillstop"<<endl; while(true){ cout<<"Pleaseenterx:"; doublex; cin>>x; if(x==SENTINEL)break; intn=(x-XNodes[0])/Hi[0];//等距節(jié)點(diǎn)時才可以用此方式 if(x==10)n--; x=(x-XNodes[n]); intN=4*n; doubley=pp[N]*x*x*x+pp[N+1]*x*x+pp[N+2]*x+pp[N+3]; doubleyprime=3*pp[N]*x*x+2*pp[N+1]*x+pp[N+2]; cout<<"Thecorrespondyis"<<y<<endl; cout<<"Thecorrespondy'is"<<yprime<<endl; cout<<'\n'; }輸入數(shù)據(jù)：計算結(jié)果：程序有子程序和主程序組成，便于后續(xù)的修改、擴(kuò)展，此題為固定結(jié)點(diǎn)樣條插值，若為自由結(jié)點(diǎn)或非結(jié)點(diǎn)樣條插值只需將自程序spanLinearSystem略微修改即可。程序通用，對任意固定結(jié)點(diǎn)樣條插值，只要將已知數(shù)據(jù)按序?qū)懭胛募纯勺詣舆M(jìn)行計算。5.newton’smethod:程序：#include<iostream>#include<cmath>usingnamespacestd;constdoubleERROR=1e-5;doubley(doublex);doubleyprime(doublex);doublepower(doublex,inta);voidnewtonIteration(double&x,double&error);intmain(){ doublex,error; cout<<"Entertheinitialvalue:"; cin>>x; while(true){ newtonIteration(x,error); if(fabs(error)<ERROR){ newtonIteration(x,error); if(fabs(error)<ERROR){ cout<<"Thesolutionis"<<x<<endl; break; } } }}voidnewtonIteration(double&x,double&error){ doublex1; x1=x-y(x)/yprime(x); error=x1-x; x=x1;}doubley(doublex){ returnpower(x,7)-28*power(x,4)+14;}doubleyprime(doublex){ return7*power(x,6)-28*4*power(x,3);}doublepower(doublex,inta){ doubles=1; for(inti=0;i<a;i++){ s*=x; } returns;}計算結(jié)果：連續(xù)兩次迭代達(dá)到誤差限，則停止迭代，保證收斂。6.Romberg程序：#include<iostream>#include<math.h>#include"vector.h"constdoubleERROR=1e-5;constdoubleBEGGIN=1;constdoubleEND=3;constdoubleH=END-BEGGIN;usingnamespacestd;doubleTh_k(int&k);voidspan_T_table(double&Th,int&k,Vector<double>&T_table_column1,Vector<double>&T_table_column2);voidspanThisColumn(Vector<double>&T_table_column1,Vector<double>&T_table_column2,intk);voidspanThisColumn(Vector<double>&T_table_column1,Vector<double>&T_table_column2,intk);booltestTolerance(int&k,Vector<double>&T_table_column1,Vector<double>&T_table_column2,double&quad);booltest(Vector<double>&T_table_column);doublefunction(double&x);intmain(){ intk=1; Vector<double>T_table_column1; Vector<double>T_table_column2; doublequad=0; boolT=false; while(true){ doubleTh=Th_k(k); //cout<<Th<<endl; span_T_table(Th,k,T_table_column1,T_table_column2); if(k>3)T=testTolerance(k,T_table_column1,T_table_column2,quad); if(T)break; k++; } cout<<quad<<endl; return0;}doubleTh_k(int&k){ doubleh=H/pow(2.0,k-1); doublesum=0; for(doublei=BEGGIN+h;i<END;i+=h){ sum+=2*function(i); } doublebeggin=BEGGIN,end=END; sum+=(function(beggin)+function(end)); sum*=(h/2); returnsum;}voidspan_T_table(d