北航数值分析作业三

《数值分析》第三次上机大作业姓名：学号：一、题目关于x，y，t，u，v，w的下列方程组0.5cost+u+v+w-x=2.67t+0.5sinu+v+w-y=1.070.5t+u+cosv+w-x=3.74t+0.5u+v+sinw-y=0.79以及关于z，t，u的下列二维数表zut00.40.81.21.620-0.5-0.340.140.942.063.50.2-0.42-0.5-0.260.31.182.380.4-0.18-0.5-0.5-0.180.461.420.60.22-0.34-0.58-0.5-0.10.620.80.78-0.02-0.5-0.66-0.5-0.021.01.50.46-0.26-0.66-0.74-0.5确定了一个二元函数z=f(x，y)。1．试用数值方法求出f(x，y)在区域D={(x,y)︱0≤x≤0.8,0.5≤y≤1.5}上的一个近似表达式,0(,)krsrsrspxycxy要求p(x,y)一最小的k值达到以下的精度10202700((,)(,))10ijijijfxypxy其中xi=0.08i，yj=0.5+0.05j。2．计算f(xi*,yj*),p(xi*,yj*)(i=1,2,…,8;j=1,2,…,5)的值，以观察p(x,y)逼近f(x,y)的效果，其中xi*=0.1i,yj*=0.5+0.2j。二、算法方案1．使用C++语言实现，使用牛顿迭代法求解非线性方程组，对ixi08.0，jyj05.05.0，(200,100ji)的共计11×21组jiyx,分别求出非线性方程组的解，即求出与,ijxy对应的,,,tuvw。,,,tuvw均为11×21的矩阵。2．由求出的jitu,，使用分片二次代数插值法对题中给出的utz,,数表进行插值得到ijz。即得到(,)Zfxy的11×21个数值解。3．k=0时的多项式拟合必然不符合要求，从k=1开始迭代，使用最小二乘法的曲面拟合法对(,)Zfxy进行拟合，计算在不符合要求的情况下增大k。当710时结束计算，输出结果。4．由3中得到的系数计算**(,)ijpxy的值，再次使用牛顿迭代法对**,ijxy进行求解得到,,,tuvw，再次进行二次插值得到结果**(,)ijfxy，以观察),(yxp逼近),(yxf的效果。其中*0.1ixi，*0.50.2jyj，(1,2,,8,1,2,,5)ij。三、源程序：#includeiostream#includevector#includecmath#includealgorithm#includeiomanip#defineN4//方程组未知个数#defineM6//z,t,u数表阶数#defineX_Num11#defineY_Num21//定义数表大小#defineEPSL1e-12//定义阶数，精度#defineEPSL21e-7usingnamespacestd;typedefvectorvectordoubleMat;//将二维数组简写为MatvectordoubleEquation(Matinput);//定义求解非线性方程的函数，同时供Inverse，Zxy函数调用MatInverse(intrank,Matinput2);//定义求解逆矩阵的函数doubleAccuracy(vectordoubleX_1,vectordoubleX_2);//定义求解近似向量精度的函数doubleInterpolation(doubleu_1,doublet_1);//定义分片代数二次插值函数MatCrs(vectordoubleX,vectordoubleY,MatU);//最小二乘法求解近似表达式系数MatZxy(vectordoubleX1,vectordoubleY1);//定义非线性方程组，调用Equation，Accuracy和Interpolation完成求解//所有的output应该调整，是否调整为输出到文件为好voidoutput(vectordoubleFinal1,vectordoubleFinal2,MatFinal3);//定义输出函数,输出矩阵voidoutput2(MatXi);doublevector_u[M]={0,0.4,0.8,1.2,1.6,2};doublevector_t[M]={0,0.2,0.4,0.6,0.8,1};doublemat_z[M][M]={{-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},};//定义初始数表z,t,u，此处使用数组，而其它矩阵和向量均使用的为vector以及二维vectorvoidmain(){inti,j;vectordoublex,y;x.resize(X_Num);y.resize(Y_Num);for(i=0;iX_Num;i++){x[i]=0.08*i;}for(i=0;iY_Num;i++){y[i]=0.5+0.05*i;}//定义插值节点MatZ=Zxy(x,y);//求出插值点函数值output(x,y,Z);MatCr=Crs(x,y,Z);//求出近似多项式output2(Cr);x.resize(8);y.resize(5);for(i=0;i8;i++){x[i]=0.1*(i+1);}for(j=0;j5;j++){y[j]=0.5+0.2*(j+1);}//新插值节点MatZ2=Zxy(x,y);MatP;P.resize(8);for(i=0;i8;i++){P[i].resize(5);}intk=Cr.size();//利用上一部的Cr获得P值即可doubletmp;intm,n;for(m=0;m8;m++){for(n=0;n5;n++){tmp=0;for(i=0;ik;i++){for(j=0;jk;j++){tmp=tmp+Cr[i][j]*pow(x[m],i)*pow(y[n],j);}}P[m][n]=tmp;}}//使用近似多项式得到的近似值for(i=0;i8;i++){for(j=0;j5;j++){coutsetprecision(2)scientificx[i];coutsetprecision(2)scientificy[j];cout插值setprecision(12)scientificZ2[i][j];cout拟合setprecision(12)scientificP[i][j]endl;}}system(pause);}MatZxy(vectordoubleX1,vectordoubleY1){inti,j,k;vectordoublex,y;x=X1;y=Y1;intX_No=x.size();intY_No=y.size();vectordoubleX_1,X_2;doublewrong;X_1.resize(N);//过渡用于判断误差X_2.resize(N);//此处调试发现x,y值略有差异MatA;//使用牛顿法迭代的带求非线性方程Matt,u,v,w;MatZ;//对应t,u的数表ZA.resize(N);for(i=0;iN;i++){A[i].resize(N+1);}t.resize(X_No);u.resize(X_No);v.resize(X_No);w.resize(X_No);Z.resize(X_No);for(i=0;iX_No;i++){t[i].resize(Y_No);u[i].resize(Y_No);v[i].resize(Y_No);w[i].resize(Y_No);Z[i].resize(Y_No);}for(i=0;iX_No;i++){for(j=0;jY_No;j++){t[i][j]=u[i][j]=w[i][j]=v[i][j]=1;}}//将待求量赋予初值for(i=0;iN;i++){for(j=0;jN;j++){A[i][j]=1;}}A[2][0]=0.5;A[3][1]=0.5;for(i=0;iX_No;i++)//求解对应x,y的t,u,v,w{j=0;while(jY_No){A[0][4]=2.67+x[i]-0.5*cos(t[i][j])-u[i][j]-v[i][j]-w[i][j];//此处应做修改A[1][4]=1.07+y[j]-0.5*sin(u[i][j])-t[i][j]-v[i][j]-w[i][j];A[2][4]=3.74+x[i]-cos(v[i][j])-0.5*t[i][j]-u[i][j]-w[i][j];A[3][4]=0.79+y[j]-sin(w[i][j])-0.5*u[i][j]-t[i][j]-v[i][j];A[0][0]=-0.5*sin(t[i][j]);A[1][1]=0.5*cos(u[i][j]);A[2][2]=-sin(v[i][j]);A[3][3]=cos(w[i][j]);vectordoubleChange=Equation(A);//调用求解方程得到第一组增量,此处需要注意Change赋值的问题，得到每组增量后怎么处理for(k=0;kN;k++){X_1[k]=Change[k];}X_2[0]=t[i][j];X_2[1]=u[i][j];X_2[2]=v[i][j];X_2[3]=w[i][j];wrong=Accuracy(X_1,X_2);if(wrong=EPSL){t[i][j]=X_2[0];u[i][j]=X_2[1];v[i][j]=X_2[2];w[i][j]=X_2[3];j++;}else{t[i][j]=X_1[0]+X_2[0];u[i][j]=X_1[1]+X_2[1];v[i][j]=X_1[2]+X_2[2];w[i][j]=X_1[3]+X_2[3];}}}for(i=0;iX_No;i++){for(j=0;jY_No;j++){Z[i][j]=Interpolation(u[i][j],t[i][j]);//在此处构建循环得到z的*21矩阵，即f(x,y)}}returnZ;}vectordoubleEquation(Matinput)//选主元的doolittle分解法，求非线性方程{Mata=input;//获得aintNum1=a.size();//a的阶数intNum2=a[0].size();intMax;vectordoubleMax2(Num1);vectordoubleX(Num1);vectordoubleS(Num1);doubleQ;//中转变量vectordoubleP(Num2);//中转变量inti,j,k,t;doubletmp;doubletmp1;//用于改善矩阵的条件数for(k=0;kNum1;k++){Max=k;for(i=k;iNum1;i++){tmp=0;for(t=0;tk;t++){tmp=tmp+a[i][t]*a[t][k];}S[i]=a[i][k]-tmp;//计算中间变量sif(abs(S[i])abs(S[Max]))//应为判断绝对值大小{Max=i;}}//寻找占优的行P=a[Max];a[Max]=a[k];a[k]=P;//