北航数值分析大作业第三题2014秋

数值分析大作业3参考前人程序基础上自己编写的，大家参考(copy)的时候又多了一个版本了，不要问我是谁，请叫我红领巾SY1405***1.算法设计方案1.1.定义矩阵、向量的类为了能将矩阵、复数、向量作为函数的返回值，也方便计算，定义矩阵类Matrix、向量类Vector，上次大作业二也采用类编写矩阵、向量，本次作业重用了上次作业的大部分代码。如Matrix类、Vector类、列主元高斯消去法求解方程𝑨𝒙=𝒃的函数VectorMainRowGauss(MatrixA,Vectorb)等。1.1.1.矩阵类变量：doubledata[Mmax][Nmax];用于存放数据。[Mmax][Nmax]为最多的行数列数，此存储方法并没有达到最优化。目的是为简化编程，由此造成空间占用量较大，如4*4阶矩阵，却申请了Mmax*Nmax的空间，其他无用空间存放零。由于没有使用全局变量，申请空间在子函数调用完后就自动被系统删除。故总的空间占用量可以接受。intm,n;//m行n列阶数，默认为Mmax*Nmax阶程序中：#defineMmax21//定义Matrix最多允许的行数#defineNmax21//定义Matrix最多允许的列数函数：Matrix(inthang=Mmax,intlie=Nmax)构造hang*lie阶零矩阵，不输入hang与lie时，hang与lie默认为Mmax，Nmaxvoidshow(ofstream&out)以精确到小数点后3位、固定小数点位置的形式在屏幕输出矩阵的结果；并以精确到小数点后11位、科学计数法的形式将结果输出到文件result.txt中Matrixoperator*(MatrixB)矩阵的乘法运算，当前矩阵this为m*n阶矩阵，B为B.m*B.n阶矩阵，this.n==B.m才可运算Matrix&operator=(Matrixa)在程序中可直接写成this=a实现矩阵a向矩阵this赋值1.1.2.向量类变量：doubledata[VectorMmax];用于存放数据。intm;m阶的列向量函数：Vector(intx=VectorMmax)构造x阶的零向量，默认x=VectorMmax程序中：#defineVectorMmax21//定义Vector最多允许的阶数Vector(MatrixA,intk)用矩阵A的第k列构造向量，注意列号k从0开始Vector&operator=(Vectora)赋值，this=a意为将向量a的数据和阶数赋给当前向量thisVectoroperator+(Vectorb)返回当前向量this+向量b的和1.2.除主要函数外方便编程的小函数VectorMainRowGauss(MatrixA,Vectorb)列主元高斯消去法求解方程𝑨𝒙=𝒃的函数，返回值为解向量x。doubleMax(Vectorx)返回向量x的无穷范数MatrixT(MatrixA)计算矩阵A的转置。Matrixinv(MatrixA)返回方阵A的逆，求逆思路为：若AB=I，B=[b1,b2,…,bn]，I=[e1,e2,…,en]，利用列主元高斯消去法求解Abi=ei，i=1,2,3,…,n即可。得出的B即是A-11.3.Newton迭代法求解非线性方程组的解VectorNewtonSolve(doublexx,doubleyy)函数参数为xx，yy的值：0.08(0,1,,10)ixii和0.50.05(0,1,,20)jyjj，返回值为tuvw的解向量Vectorx。采用函数Vector_F(Vectorx,doublexx,doubleyy)求解−𝑭(𝒙(𝑘))，采用函数MatrixJ(Vectorx)求解Jaccobi矩阵𝑭′(𝒙(𝑘))，最后用VectorMainRowGauss(MatrixA,Vectorb)求解线性方程组𝑭′(𝒙(𝑘))∆𝑥(𝑘)=−𝑭(𝒙(𝑘))1.4.分片二次插值doublez(doublet,doubleu)由上一步求得的t,u的值在题中所给z-tu数表中分片二次插值得到z。具体做法为：采取分片二次插值，选择（m，n）满足,2322mmhhtttm,2322nnuuun如果12htt或42htt，则m=1或4；如果12uu或42uu，则n=1或n=4。选择(,)(1,,1;1,,1)krtukmmmrnnn为插值节点，相应的Lagrange形式的插值多项式为112211(,)()()(,)mnkrkrkmrnptultluftu其中11()mwkwmkwwkttlttt(k=m-1,m,m+1)11()nwrwnrwwryyluyy(r=n-1,n,n+1)综合以上两步，得到函数doublef(doublex,doubley)，即描述z与x,y关系的函数z=f(x，y)。于是得到了数表,,(,)ijijxyfxy。1.5.曲面拟合按照教材P142-143的方法拟合数表,,(,)ijijxyfxy。分为以下两个子函数：1）求拟合为k阶多项式时的矩阵C：MatrixC(intk,Matrix&U)2）得到系数矩阵C后得到拟合曲面,0(,)krsrsrspxycxyk从2递增，每步验证是否满足精度条件10202700((,)(,))10ijijijfxypxy2.源代码以下C++代码在VS2010下编译通过#includeiostream#includefstream#includeiomanipusingnamespacestd;#defineMmax21//定义Matrix最多允许的行数#defineNmax21//定义Matrix最多允许的列数#defineVectorMmax21//定义Vector最多允许的阶数#defineepsilon1e-12//题中要求的精度classMatrix{public:doubledata[Mmax][Nmax];intm,n;//m行n列Matrix(inthang=Mmax,intlie=Nmax){m=hang;n=lie;for(inti=0;iMmax;i++){for(intj=0;jNmax;j++){data[i][j]=0;}}}voidshow(ofstream&out){for(inti=0;im;i++){for(intj=0;jn;j++){//setiosflags(ios::fixed)表示不用科学计数法表示，即0.00001记为0.000而不是1.000e-5coutsetprecision(3)setiosflags(ios::fixed)data[i][j]\t;}coutresetiosflags(ios::fixed)endl;//coutendl;//若不resetiosflags(ios::fixed)则出现错误}cout-----------------endl;for(inti=0;im;i++){for(intj=0;jn;j++){outsetprecision(11)setw(20)setiosflags(ios::scientific)data[i][j]\t;}outendl;}out-----------------endl;}Matrixoperator*(MatrixB){MatrixResult(m,B.n);if(n!=B.m){cout当前m*n阶矩阵与运算符*后的矩阵不能相乘！\n;returnResult;}for(inti=0;im;i++)for(intj=0;jB.n;j++){doubletemp=0;for(intt=0;tn;t++)temp+=data[i][t]*B.data[t][j];Result.data[i][j]=temp;}returnResult;}Matrix&operator=(Matrixa)//this的无用单元（即超过行号列号的元素）不改变值{m=a.m;n=a.n;for(inti=0;im;i++)for(intj=0;jn;j++)data[i][j]=a.data[i][j];return*this;}};classVector{public:doubledata[VectorMmax];intm;Vector(intx=VectorMmax){m=x;for(inti=0;iVectorMmax;i++)data[i]=0;}Vector(MatrixA,intk)//用A的第k列构造向量，注意列号k从0开始{m=A.m;for(inti=0;iA.m;i++)data[i]=A.data[i][k];}Vector&operator=(Vectora)//this的无用单元的值没有改变{m=a.m;for(inti=0;im;i++)data[i]=a.data[i];return*this;}Vectoroperator+(Vectorb){VectorResult(m);if(m!=b.m){cout阶数不同，向量与向量不能相减！\n;returnResult;}for(inti=0;im;i++)Result.data[i]=data[i]+b.data[i];returnResult;}};VectorMainRowGauss(MatrixA,Vectorb)//解Ax=b，A须方阵，b.m==A.m==A.n，返回值为x向量//完成后A、b都没有改变（改变的只是形参，实参不便）{Vectorx(b.m);if((A.m!=A.n)||(A.n!=b.m)){coutA不是方阵或执行MainRowGauss解Ax=b时方阵A与b的阶数不同，无法求解\n;returnx;}//消元过程for(intk=0;kA.n-1;k++){intwheremax=k;for(inti=k+1;iA.n;i++)if(A.data[i][k]A.data[wheremax][k])wheremax=i;doubletemp;for(intj=k;jA.n;j++){temp=A.data[wheremax][j];A.data[wheremax][j]=A.data[k][j];A.data[k][j]=temp;}temp=b.data[wheremax];b.data[wheremax]=b.data[k];b.data[k]=temp;for(inti=k+1;iA.n;i++){doublemik=A.data[i][k]/A.data[k][k];for(intj=k+1;jA.n;j++)A.data[i][j]=A.data[i][j]-mik*A.data[k][j];b.data[i]=b.data[i]-mik*b.data[k];}}//回代过程x.data[b.m-1]=b.data[b.m-1]/A.data[A.n-1][A.n-1];//b.m==A.m==A.nfor(intk=A.n-2;k=0;k--){doublesum=0;for(intj=k+1;jA.n;j++)sum+=A.data[k][j]*x.data[j];x.data[k]=(b.data[k]-sum)/A.data[k][k];}returnx;}Vector_F(Vectorx,doublexx,doubleyy)/*返回-F(x)*/{VectorF(4);F.data[0]=-(0.5*cos(x.data[0])+x.data[1]+x.data[2]+x.data[3]-xx-2.67);F.data[1]=-(x.data[0]+0.5*sin