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

1、数值分析大作业3参考前人程序基础上自己编写的，大家参考(copy)的时候又多了一个版本了，不要问我是谁，请叫我红领巾SY1405*1. 算法设计方案1.1. 定义矩阵、向量的类为了能将矩阵、复数、向量作为函数的返回值，也方便计算，定义矩阵类Matrix、向量类Vector，上次大作业二也采用类编写矩阵、向量，本次作业重用了上次作业的大部分代码。如Matrix类、Vector类、列主元高斯消去法求解方程Ax=b的函数Vector MainRowGauss(Matrix A,Vector b)等。1.1.1. 矩阵类变量：double dataMmaxNmax;用于存放数据。MmaxNmax为最多

2、的行数列数，此存储方法并没有达到最优化。目的是为简化编程，由此造成空间占用量较大，如4*4阶矩阵，却申请了Mmax*Nmax的空间，其他无用空间存放零。由于没有使用全局变量，申请空间在子函数调用完后就自动被系统删除。故总的空间占用量可以接受。int m,n;/m行n列阶数，默认为Mmax*Nmax阶程序中：#define Mmax 21/定义Matrix最多允许的行数#define Nmax 21/定义Matrix最多允许的列数函数：Matrix(int hang=Mmax,int lie=Nmax)构造hang*lie阶零矩阵，不输入hang与lie时，hang与lie默认为Mmax，Nma

3、xvoid show(ofstream &out)以精确到小数点后3位、固定小数点位置的形式在屏幕输出矩阵的结果；并以精确到小数点后11位、科学计数法的形式将结果输出到文件result.txt中Matrix operator*(Matrix B)矩阵的乘法运算，当前矩阵this为m*n阶矩阵，B为B.m*B.n阶矩阵，this.n=B.m才可运算Matrix& operator=(Matrix a)在程序中可直接写成this=a实现矩阵a向矩阵this赋值1.1.2. 向量类变量：double dataVectorMmax;用于存放数据。int m;m阶的列向量函数：Vecto

4、r(int x=VectorMmax)构造x阶的零向量，默认x=VectorMmax程序中：#define VectorMmax 21/定义Vector最多允许的阶数Vector(Matrix A,int k) 用矩阵A的第k列构造向量，注意列号k从0开始Vector& operator=(Vector a) 赋值，this=a意为将向量a的数据和阶数赋给当前向量thisVector operator+(Vector b)返回当前向量this+向量b的和1.2. 除主要函数外方便编程的小函数Vector MainRowGauss(Matrix A,Vector b)列主元高斯消去法求解

5、方程Ax=b的函数，返回值为解向量x。double Max(Vector x)返回向量x的无穷范数Matrix T(Matrix A)计算矩阵A的转置。Matrix inv(Matrix A)返回方阵A的逆，求逆思路为：若AB=I，B=b1,b2,bn，I=e1,e2,en，利用列主元高斯消去法求解Abi=ei，i=1,2,3,n即可。得出的B即是A-11.3. Newton迭代法求解非线性方程组的解Vector NewtonSolve(double xx,double yy)函数参数为xx，yy的值：和，返回值为tuvw的解向量Vector x。采用函数Vector _F(Vector x,

6、double xx,double yy)求解-F(xk)，采用函数Matrix J(Vector x)求解Jaccobi矩阵F'(x(k)，最后用Vector MainRowGauss(Matrix A,Vector b)求解线性方程组F'xkxk=-F(xk)1.4. 分片二次插值double z(double t,double u)由上一步求得的t,u的值在题中所给z-tu数表中分片二次插值得到z。具体做法为：采取分片二次插值，选择（m，n）满足如果或，则m=1或4；如果或，则n=1或n=4。选择为插值节点，相应的Lagrange形式的插值多项式为其中 (k=m-1, m,

7、 m+1) (r=n-1, n, n+1)综合以上两步，得到函数double f(double x,double y)，即描述z与x, y关系的函数z=f(x，y)。于是得到了数表。1.5. 曲面拟合按照教材P142-143的方法拟合数表。分为以下两个子函数：1） 求拟合为k阶多项式时的矩阵C：Matrix C(int k,Matrix &U)2） 得到系数矩阵C后得到拟合曲面k从2递增，每步验证是否满足精度条件2. 源代码以下C+代码在VS2010下编译通过#include<iostream>#include<fstream>#include<ioman

8、ip>using namespace std;#define Mmax 21/定义Matrix最多允许的行数#define Nmax 21/定义Matrix最多允许的列数#define VectorMmax 21/定义Vector最多允许的阶数#define epsilon 1e-12/题中要求的精度class Matrixpublic:double dataMmaxNmax;int m,n;/m行n列Matrix(int hang=Mmax,int lie=Nmax)m=hang;n=lie;for(int i=0;i<Mmax;i+)for(int j=0;j<Nmax;

9、j+)dataij=0;void show(ofstream &out)for(int i=0;i<m;i+)for(int j=0;j<n;j+)/setiosflags(ios:fixed)表示不用科学计数法表示，即0.00001记为0.000而不是1.000e-5cout<<setprecision(3)<<setiosflags(ios:fixed)<<dataij<<"t"cout<<resetiosflags(ios:fixed)<<endl;/cout<<e

10、ndl;/若不resetiosflags(ios:fixed)则出现错误cout<<"-"<<endl;for(int i=0;i<m;i+)for(int j=0;j<n;j+)out<<setprecision(11)<<setw(20)<<setiosflags(ios:scientific)<<dataij<<"t"out<<endl;out<<"-"<<endl;Matrix operator*

11、(Matrix B)Matrix Result(m,B.n);if(n!=B.m)cout<<"当前m*n阶矩阵与运算符*后的矩阵不能相乘！n"return Result;for(int i=0;i<m;i+)for(int j=0;j<B.n;j+)double temp=0;for(int t=0;t<n;t+)temp+=datait*B.datatj;Result.dataij=temp;return Result;Matrix& operator=(Matrix a)/this的无用单元（即超过行号列号的元素）不改变值m=a.

12、m;n=a.n;for(int i=0;i<m;i+)for(int j=0;j<n;j+)dataij=a.dataij;return *this;class Vectorpublic:double dataVectorMmax;int m;Vector(int x=VectorMmax)m=x;for(int i=0;i<VectorMmax;i+)datai=0;Vector(Matrix A,int k)/用A的第k列构造向量，注意列号k从0开始m=A.m;for(int i=0;i<A.m;i+)datai=A.dataik;Vector& opera

13、tor=(Vector a)/this的无用单元的值没有改变m=a.m;for(int i=0;i<m;i+)datai=a.datai;return *this;Vector operator+(Vector b)Vector Result(m);if(m!=b.m)cout<<"阶数不同，向量与向量不能相减！n"return Result;for(int i=0;i<m;i+)Result.datai=datai+b.datai;return Result;Vector MainRowGauss(Matrix A,Vector b)/解Ax=b，

14、A须方阵，b.m=A.m=A.n，返回值为x向量/完成后A、b都没有改变（改变的只是形参，实参不便）Vector x(b.m);if(A.m!=A.n)|(A.n!=b.m)cout<<"A不是方阵 或 执行MainRowGauss解Ax=b时方阵A与b的阶数不同，无法求解n"return x;/消元过程for(int k=0;k<A.n-1;k+)int wheremax=k;for(int i=k+1;i<A.n;i+)if(A.dataik>A.datawheremaxk)wheremax=i;double temp;for(int j=

15、k;j<A.n;j+)temp=A.datawheremaxj;A.datawheremaxj=A.datakj;A.datakj=temp;temp=b.datawheremax;b.datawheremax=b.datak;b.datak=temp;for(int i=k+1;i<A.n;i+)double mik=A.dataik/A.datakk;for(int j=k+1;j<A.n;j+)A.dataij=A.dataij-mik*A.datakj;b.datai=b.datai-mik*b.datak;/回代过程x.datab.m-1=b.datab.m-1/A

16、.dataA.n-1A.n-1;/b.m=A.m=A.nfor(int k=A.n-2;k>=0;k-)double sum=0;for(int j=k+1;j<A.n;j+)sum+=A.datakj*x.dataj;x.datak=(b.datak-sum)/A.datakk;return x;Vector _F(Vector x,double xx,double yy)/*返回-F(x)*/Vector F(4);F.data0=-(0.5*cos(x.data0)+x.data1+x.data2+x.data3-xx-2.67);F.data1=-(x.data0+0.5*

17、sin(x.data1)+x.data2+x.data3-yy-1.07);F.data2=-(0.5*x.data0+x.data1+cos(x.data2)+x.data3-xx-3.74);F.data3=-(x.data0+0.5*x.data1+x.data2+sin(x.data3)-yy-0.79); return F;Matrix J(Vector x)/*求解Jaccobi矩阵F'(x)*/Matrix Result(4,4);Result.data00=-0.5*sin(x.data0); Result.data01=1;Result.data02=1;Result

18、.data03=1;Result.data10=1;Result.data11=0.5*cos(x.data1);Result.data12=1;Result.data13=1;Result.data20=0.5;Result.data21=1;Result.data22=-1*sin(x.data2);Result.data23=1;Result.data30=1;Result.data31=0.5;Result.data32=1;Result.data33=cos(x.data3);return Result;double Max(Vector x)/返回x的无穷范数double max=

19、0;for(int i=0;i<x.m;i+)if(fabs(x.datai)>max)max=fabs(x.datai);return max;Vector NewtonSolve(double xx,double yy)/返回值为tuvw的解向量Vector xVector x(4);/t,u,v,w，初值为1,1,1,1for(int i=0;i<4;i+)x.datai=1;Vector dx(4);/x，初值为零for(int k=0;k<1000;k+)dx=MainRowGauss(J(x),_F(x,xx,yy);if(Max(dx)/Max(x)<

20、;=epsilon)/|x|/|x|<=epsilonreturn x;else x=x+dx;double z(double t,double u)/由t,u分片二次插值得到z的值int m,n;double u06=0,0.4,0.8,1.2,1.6,2;double t06=0,0.2,0.4,0.6,0.8,1;double lt3,lu3;double ztu66=-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.

21、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,;if(t<=t01+0.1) m=1;else if(t>t04-0.1) m=4;else for(m=2;m<=3;m+)if(t0m-0.1<t)&&(t<=t0m+0.1) break;if(t<=u01+0.2) n=1;else if(t>u04-0.2) n=4;else for(n=2;n<=3;n+)if(u0n-0.2<t)&&am

22、p;(t<=u0n+0.2) break;/*计算拉格朗日插值基函数*/lt0=(t-t0m)*(t-t0m+1)/(t0m-1-t0m)/(t0m-1-t0m+1);lt1=(t-t0m-1)*(t-t0m+1)/(t0m-t0m-1)/(t0m-t0m+1);lt2=(t-t0m-1)*(t-t0m)/(t0m+1-t0m-1)/(t0m+1-t0m);lu0=(u-u0n)*(u-u0n+1)/(u0n-1-u0n)/(u0n-1-u0n+1); lu1=(u-u0n-1)*(u-u0n+1)/(u0n-u0n-1)/(u0n-u0n+1);lu2=(u-u0n-1)*(u-u0n

23、)/(u0n+1-u0n-1)/(u0n+1-u0n);double result=0;for(int k=m-1;k<=m+1;k+)for(int r=n-1;r<=n+1;r+)result+=ltk-m+1*lur-n+1*ztukr; /*计算插值函数值*/return result;double f(double x,double y)Vector tuvw(4);tuvw=NewtonSolve(x,y);return z(tuvw.data0,tuvw.data1);Matrix T(Matrix A)Matrix Result(A.n,A.m);for(int i

24、=0;i<A.m;i+)for(int j=0;j<A.n;j+)Result.dataji=A.dataij;return Result;Matrix inv(Matrix A)/A须为方阵,A.m=A.nMatrix B(A.m,A.n);if(A.m!=A.n) return B;for(int k=0;k<A.m;k+)Vector e(A.m);/定义单位向量ek=(00.00100.0)e.datak=1;Vector x=MainRowGauss(A,e);for(int t=0;t<A.m;t+)B.datatk=x.datat;return B;Mat

25、rix C(int k,Matrix &U)/拟合为k次多项式时的C,U(11,21)为避免重复计算而引入k+;/k次多项式时，B、G的列数为k+1Matrix B(11,k);Matrix G(21,k);for(int j=0;j<k;j+)for(int i=0;i<11;i+)B.dataij=pow(0.08*i,j);for(int j=0;j<k;j+)for(int i=0;i<21;i+)G.dataij=pow(0.5+0.05*i,j);return inv(T(B)*B)*T(B)*U*G*inv(T(G)*G);/返回值为(k+1)*(

26、k+1)阶的Cdouble p(double x,double y,Matrix C)/C为(k+1)*(k+1)阶方阵double result=0;int k=C.m-1;for(int r=0;r<=k;r+)for(int s=0;s<=k;s+)result+=C.datars*pow(x,r)*pow(y,s);return result;void main()ofstream out("result.txt");cout<<"xi,yj,f(xi,yj)数表为n"Matrix U(11,21);/用于储存f(xi,y

27、j)的值，避免重复计算Vector xi(11),yj(21);for(int i=0;i<=10;i+)xi.datai=0.08*i;for(int j=0;j<=20;j+)yj.dataj=0.5+0.05*j;for(int i=0;i<=10;i+)for(int j=0;j<=20;j+)U.dataij=f(xi.datai,yj.dataj);for(int i=0;i<=10;i+)for(int j=0;j<=20;j+)cout<<"x"<<i<<"=t"&l

28、t;<xi.datai<<"t"<<"y"<<j<<"=t"<<yj.dataj<<"t"<<setiosflags(ios:scientific)<<setprecision(12)<<U.dataij<<resetiosflags(ios:scientific)<<endl;int k;Matrix Crs;for(k=2;k<10;k+)Crs=C(k,U);doub

29、le sigma=0;for(int i=0;i<=10;i+)for(int j=0;j<=20;j+)sigma+=pow(U.dataij-p(xi.datai,yj.dataj,Crs),2);cout<<"k="<<k<<"t"<<"="<<setiosflags(ios:scientific)<<setprecision(12)<<sigma<<endl;if(sigma<=1e-7) break;cout&l

30、t;<resetiosflags(ios:scientific);Crs.show(out);for(int i=1;i<=8;i+)for(int j=1;j<=5;j+)double xxing=0.1*i,yxing=0.5+0.2*j;cout<<"x"<<i<<"*=t"<<xxing<<"t"<<"y"<<j<<"*=t"<<yxing<<&quo

31、t;t"<<setiosflags(ios:scientific)<<setprecision(12)<<f(xxing,yxing)<<"t"<<p(xxing,yxing,Crs)<<resetiosflags(ios:scientific)<<endl;3. 运行结果注：以下带框的段落为从控制台直接复制出来的结果。数表如下：xi,yj,f(xi,yj)数表为x0= 0 y0= 0.5 4.465040184807e-001x0= 0 y1= 0.55 3.2468326292

32、77e-001x0= 0 y2= 0.6 2.101596866827e-001x0= 0 y3= 0.65 1.030436083160e-001x0= 0 y4= 0.7 3.401895562675e-003x0= 0 y5= 0.75 -8.873581363800e-002x0= 0 y6= 0.8 -1.733716327497e-001x0= 0 y7= 0.85 -2.505346114666e-001x0= 0 y8= 0.9 -3.202765063876e-001x0= 0 y9= 0.95 -3.826680661097e-001x0= 0 y10= 1 -4.3779

33、57667384e-001x0= 0 y11= 1.05 -4.857589414438e-001x0= 0 y12= 1.1 -5.266672548835e-001x0= 0 y13= 1.15 -5.606384797965e-001x0= 0 y14= 1.2 -5.877965387677e-001x0= 0 y15= 1.25 -6.082697790899e-001x0= 0 y16= 1.3 -6.221894528764e-001x0= 0 y17= 1.35 -6.296883781856e-001x0= 0 y18= 1.4 -6.308997600028e-001x0=

34、 0 y19= 1.45 -6.259561525454e-001x0= 0 y20= 1.5 -6.149885466094e-001x1= 0.08 y0= 0.5 6.380152265113e-001x1= 0.08 y1= 0.55 5.066117551467e-001x1= 0.08 y2= 0.6 3.821763692774e-001x1= 0.08 y3= 0.65 2.648634911537e-001x1= 0.08 y4= 0.7 1.547802002848e-001x1= 0.08 y5= 0.75 5.199268349094e-002x1= 0.08 y6=

35、0.8 -4.346804020490e-002x1= 0.08 y7= 0.85 -1.316010567885e-001x1= 0.08 y8= 0.9 -2.124310883088e-001x1= 0.08 y9= 0.95 -2.860045510580e-001x1= 0.08 y10= 1 -3.523860789794e-001x1= 0.08 y11= 1.05 -4.116554565222e-001x1= 0.08 y12= 1.1 -4.639049115188e-001x1= 0.08 y13= 1.15 -5.092367247005e-001x1= 0.08 y1

36、4= 1.2 -5.477611179623e-001x1= 0.08 y15= 1.25 -5.795943883391e-001x1= 0.08 y16= 1.3 -6.048572588895e-001x1= 0.08 y17= 1.35 -6.236734213318e-001x1= 0.08 y18= 1.4 -6.361682484133e-001x1= 0.08 y19= 1.45 -6.424676566901e-001x1= 0.08 y20= 1.5 -6.426971026996e-001x2= 0.16 y0= 0.5 8.400813957666e-001x2= 0.

37、16 y1= 0.55 6.997641656732e-001x2= 0.16 y2= 0.6 5.660614423517e-001x2= 0.16 y3= 0.65 4.391716081176e-001x2= 0.16 y4= 0.7 3.192421380408e-001x2= 0.16 y5= 0.75 2.063761923874e-001x2= 0.16 y6= 0.8 1.006385238914e-001x2= 0.16 y7= 0.85 2.060740067837e-003x2= 0.16 y8= 0.9 -8.935402476698e-002x2= 0.16 y9=

38、0.95 -1.736269688648e-001x2= 0.16 y10= 1 -2.507999561599e-001x2= 0.16 y11= 1.05 -3.209322694446e-001x2= 0.16 y12= 1.1 -3.840977350046e-001x2= 0.16 y13= 1.15 -4.403821754175e-001x2= 0.16 y14= 1.2 -4.898811523126e-001x2= 0.16 y15= 1.25 -5.326979655338e-001x2= 0.16 y16= 1.3 -5.689418792921e-001x2= 0.16

39、 y17= 1.35 -5.987265495151e-001x2= 0.16 y18= 1.4 -6.221686297503e-001x2= 0.16 y19= 1.45 -6.393865356972e-001x2= 0.16 y20= 1.5 -6.504993507878e-001x3= 0.24 y0= 0.5 1.051515091803e+000x3= 0.24 y1= 0.55 9.029274308310e-001x3= 0.24 y2= 0.6 7.605802668596e-001x3= 0.24 y3= 0.65 6.247151981456e-001x3= 0.24

40、 y4= 0.7 4.955197560009e-001x3= 0.24 y5= 0.75 3.731340427746e-001x3= 0.24 y6= 0.8 2.576567488723e-001x3= 0.24 y7= 0.85 1.491505594102e-001x3= 0.24 y8= 0.9 4.764698677337e-002x3= 0.24 y9= 0.95 -4.684932320146e-002x3= 0.24 y10= 1 -1.343567603849e-001x3= 0.24 y11= 1.05 -2.149133449274e-001x3= 0.24 y12=

41、 1.1 -2.885737006348e-001x3= 0.24 y13= 1.15 -3.554063647857e-001x3= 0.24 y14= 1.2 -4.154913964886e-001x3= 0.24 y15= 1.25 -4.689182499695e-001x3= 0.24 y16= 1.3 -5.157838831247e-001x3= 0.24 y17= 1.35 -5.561910752001e-001x3= 0.24 y18= 1.4 -5.902469305629e-001x3= 0.24 y19= 1.45 -6.180615482412e-001x3= 0

42、.24 y20= 1.5 -6.397468392579e-001x4= 0.32 y0= 0.5 1.271246751483e+000x4= 0.32 y1= 0.55 1.115002018147e+000x4= 0.32 y2= 0.6 9.646077272157e-001x4= 0.32 y3= 0.65 8.203473694751e-001x4= 0.32 y4= 0.7 6.824476781795e-001x4= 0.32 y5= 0.75 5.510852085975e-001x4= 0.32 y6= 0.8 4.263923859018e-001x4= 0.32 y7=

43、 0.85 3.084629956332e-001x4= 0.32 y8= 0.9 1.973571296919e-001x4= 0.32 y9= 0.95 9.310562085941e-002x4= 0.32 y10= 1 -4.285992234035e-003x4= 0.32 y11= 1.05 -9.483392529689e-002x4= 0.32 y12= 1.1 -1.785729903640e-001x4= 0.32 y13= 1.15 -2.555537790546e-001x4= 0.32 y14= 1.2 -3.258401501575e-001x4= 0.32 y15

44、= 1.25 -3.895069883634e-001x4= 0.32 y16= 1.3 -4.466382045995e-001x4= 0.32 y17= 1.35 -4.973249517677e-001x4= 0.32 y18= 1.4 -5.416640326994e-001x4= 0.32 y19= 1.45 -5.797564797951e-001x4= 0.32 y20= 1.5 -6.117062881476e-001x5= 0.4 y0= 0.5 1.498321052482e+000x5= 0.4 y1= 0.55 1.334998632066e+000x5= 0.4 y2

45、= 0.6 1.177125123739e+000x5= 0.4 y3= 0.65 1.025024055020e+000x5= 0.4 y4= 0.7 8.789600231744e-001x5= 0.4 y5= 0.75 7.391451087035e-001x5= 0.4 y6= 0.8 6.057448714871e-001x5= 0.4 y7= 0.85 4.788838610666e-001x5= 0.4 y8= 0.9 3.586506258818e-001x5= 0.4 y9= 0.95 2.451022361964e-001x5= 0.4 y10= 1 1.382683509

46、285e-001x5= 0.4 y11= 1.05 3.815486540699e-002x5= 0.4 y12= 1.1 -5.525282116814e-002x5= 0.4 y13= 1.15 -1.419868808137e-001x5= 0.4 y14= 1.2 -2.220944390959e-001x5= 0.4 y15= 1.25 -2.956352324598e-001x5= 0.4 y16= 1.3 -3.626795115028e-001x5= 0.4 y17= 1.35 -4.233061642240e-001x5= 0.4 y18= 1.4 -4.7760103613

47、25e-001x5= 0.4 y19= 1.45 -5.256554266672e-001x5= 0.4 y20= 1.5 -5.675647436551e-001x6= 0.48 y0= 0.5 1.731892740383e+000x6= 0.48 y1= 0.55 1.562034577209e+000x6= 0.48 y2= 0.6 1.397216918208e+000x6= 0.48 y3= 0.65 1.237801006739e+000x6= 0.48 y4= 0.7 1.084087532678e+000x6= 0.48 y5= 0.75 9.363227723149e-001x6= 0.48 y6= 0.8 7.947044490537e-001x6= 0.48 y7= 0.85 6.593871980282e-001x6= 0.48 y8= 0.9 5.304875868400e-001x6= 0.48 y9= 0.95 4.08088685