北航数值分析大作业(三)Word文件下载.docx

1、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+语言实现，使用牛顿迭代法求解非线性方程组，对，（）的共计1121组分别求出非线性方程组的解，即求出与对应的。均为1121的矩阵。2由求出的，使用分片二次代数插值法对题中给出的数表进行插值得到。即得到的1121个数值解。3k=0时的多项式拟合必然不符合要求，从k=1开始迭代，使用最小二乘法的曲面拟合法对进行拟合，计算在不符合要求的情况下增大。当时结束计算，输出结果。4由3中得到的系数计算的值，再次使

2、用牛顿迭代法对进行求解得到，再次进行二次插值得到结果，以观察逼近的效果。其中，。三、源程序：#includevectorcmathalgorithmiomanip#define N 4/方程组未知个数#define M 6/z,t,u数表阶数#define X_Num 11#define Y_Num 21/定义数表大小#define EPSL 1e-12/定义阶数，精度#define EPSL2 1e-7using namespace std;typedef vectorvector Mat;/将二维数组简写为Mat Equation（ Mat input ）;/定义求解非线性方程的函数，同时

3、供Inverse，Zxy函数调用Mat Inverse（int rank, Mat input2）;/定义求解逆矩阵的函数double Accuracy （ vector X_1,vector X_2）;/定义求解近似向量精度的函数double Interpolation （ double u_1,double t_1）;/定义分片代数二次插值函数Mat Crs（vectorX,vectorY,Mat U）;/最小二乘法求解近似表达式系数Mat Zxy（vectorX1,vectorY1）;/定义非线性方程组，调用Equation，Accuracy和Interpolation完成求解/所有的o

4、utput应该调整，是否调整为输出到文件为好void output（vectorFinal1,vectorFinal2,Mat Final3）;/定义输出函数,输出矩阵void output2（Mat Xi）;double vector_uM = 0,0.4,0.8,1.2,1.6,2;double vector_tM = 0,0.2,0.4,0.6,0.8,1;double mat_zMM = -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

5、,-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以及二维vectorvoid main（）int i,j;vector x,y;x.resize（X_Num）;y.resize（Y_Num）;for（i=0;iX_Num;i+）xi=0.08*i;Y_Num;yi=0.5+0.05*i;/定义插值节点Mat Z = Zxy（x,y）;/求出插值点函数值output（x,y,

6、Z）;Mat Cr = Crs（x,y,Z）;/求出近似多项式output2（Cr）;x.resize（8）;y.resize（5）;8;xi=0.1*（i+1）;for（j=0;j5;j+）yj=0.5+0.2*（j+1）;/新插值节点Mat Z2 = Zxy（x,y）;Mat P;P.resize（8）;Pi.resize（5）;int k = Cr.size（）;/利用上一部的Cr获得P值即可double tmp;int m,n;for（m=0;mm+）for（n=0;nn+）tmp = 0;for（i=0;k;for（j=0;tmp = tmp + Crij*pow（xm,i）*pow

7、（yn,j）;Pmn = tmp;/使用近似多项式得到的近似值for（j=0;coutsetprecision（2）scientificxi ;yj插值setprecision（12）Z2ij拟合Pijendl;system（pause）;Y1）int i,j,k;x,y;x=X1;y=Y1;int X_No=x.size（）;int Y_No=y.size（）; X_1,X_2;double wrong;X_1.resize（N）;/过渡用于判断误差X_2.resize（N）;/此处调试发现x,y值略有差异Mat A;/使用牛顿法迭代的带求非线性方程Mat t,u,v,w;Mat Z;/对应

8、t,u的数表ZA.resize（N）;N;Ai.resize（N+1）;t.resize（X_No）;u.resize（X_No）;v.resize（X_No）;w.resize（X_No）;Z.resize（X_No）;X_No;ti.resize（Y_No）;ui.resize（Y_No）;vi.resize（Y_No）;wi.resize（Y_No）;Zi.resize（Y_No）;Y_No;tij=uij=wij=vij=1;/将待求量赋予初值Aij=1;A20=0.5;A31=0.5;i+）/求解对应x,y的t,u,v,wj=0;while（jY_No）A04=2.67 + xi -

9、 0.5 * cos （tij） - uij - vij - wij;/此处应做修改A14=1.07 + yj - 0.5 * sin （uij） - tij - vij - wij;A24=3.74 + xi - cos （vij） - 0.5 * tij - uij - wij;A34=0.79 + yj - sin （wij） - 0.5 * uij - tij - vij;A00=-0.5 * sin （tij）;A11=0.5 * cos （uij）;A22=-sin （vij）;A33=cos （wij）;vector Change = Equation （ A ）; /调用求解方

10、程得到第一组增量,此处需要注意Change赋值的问题，得到每组增量后怎么处理for（k=0;kk+）X_1k = Changek;X_20=tij;X_21=uij;X_22=vij;X_23=wij;wrong = Accuracy（X_1,X_2）;if（wrong=EPSL）tij=X_20;uij=X_21;vij=X_22;wij=X_23;j+;elsetij=X_10+X_20;uij=X_11+X_21;vij=X_12+X_22;wij=X_13+X_23;Zij = Interpolation（uij,tij）;/在此处构建循环得到z的*21矩阵，即f（x,y）return

11、 Z; Equation （Mat input）/选主元的doolittle分解法，求非线性方程 Mat a = input;/获得a int Num1 = a.size（）;/a的阶数 int Num2 = a0.size（）; int Max; vector Max2（Num1）; X（Num1）; S（Num1）; double Q;/中转变量 P（Num2）; int i,j,k,t; double tmp; double tmp1;/用于改善矩阵的条件数 for（ k=0 ; kNum1 ; k+） Max=k; for（ i=k ; i i+） tmp=0; for （ t=0 ;

12、 tabs（SMax）/应为判断绝对值大小 Max=i; /寻找占优的行 P=aMax; aMax=ak; ak=P;/交换行 Q=SMax; SMax=Sk; Sk=Q;/交换S akk=Sk; for（ j=k+1 ; j j+） tmp = 0; for（ t=0 ; k ;t+） tmp = tmp + akt * atj; akj = akj - tmp; for （ i=k+1 ; aik = Si / akk; for（ i=1 ; double tmp = 0; for （ t=0 ;t=0 ; i-） tmp = 0; for（ j=i+1 ; tmp = tmp + aij

13、 * Xj; Xi = （aiNum2-1 - tmp） / aii; return X; X_2）/迭代精度 A = X_1; B = X_2;int Num1 = A.size（）;double tmp1=0;double tmp2=0;tmp1=abs（A0）;for（i=1;tmp1=max（tmp1,abs（Ai）;/增量的无穷范数tmp2=abs（B0）;tmp2=max（tmp2,abs（Bi）;/迭代向量的无穷范数double Final;Final = tmp1/tmp2;return Final;double Interpolation （ double u_1,doubl

14、e t_1）/分片插值函数double u = u_1;double t = t_1;double z=0; tmp1; tmp2;tmp1.resize（3）;tmp2.resize（3）;int k,r;if（t=0.3）k=0;if（0.3t）&（t=0.5）k=1;if（0.50.7）k=3;/判断t的范围if（u=0.6）r=0;if（0.6u）&（u=1.0）r=1;if（1.01.4）r=3;tmp10=（t-vector_tk+1）*（t-vector_tk+2）/（vector_tk-vector_tk+1）*（vector_tk-vector_tk+2）;tmp11=（t-

15、vector_tk）*（t-vector_tk+2）/（vector_tk+1-vector_tk）*（vector_tk+1-vector_tk+2）;tmp12=（t-vector_tk）*（t-vector_tk+1）/（vector_tk+2-vector_tk+1）*（vector_tk+2-vector_tk）;tmp20=（u-vector_ur+1）*（u-vector_ur+2）/（vector_ur-vector_ur+1）*（vector_ur-vector_ur+2）;tmp21=（u-vector_ur）*（u-vector_ur+2）/（vector_ur+1-ve

16、ctor_ur）*（vector_ur+1-vector_ur+2）;tmp22=（u-vector_ur）*（u-vector_ur+1）/（vector_ur+2-vector_ur+1）*（vector_ur+2-vector_ur）;3;tmp=1;tmp=tmp*mat_zk+ir+j*tmp1i*tmp2j;z=z+tmp;return z;Mat Inverse（ int rank, Mat input2 ）/用于矩阵求逆，并返回逆矩阵Mat Copy = input2;int R = rank;Mat X;X.resize（R）;R;Xi.resize（R）;for（k=0;i

17、f（k=i）Copyk.push_back（1）;Copyk.push_back（0）;/该法求解逆矩阵的原理，为在input2方阵后，添加列向量，全体列向量组成单位阵，则各解向量组成逆矩阵vector x = Equation（Copy）;Xji=xj;Copyk.pop_back（）;return X;Y,Mat U）/此处可能重定义,XY为x,y的坐标x = X;y = Y;int X_No = x.size（）;int Y_No = y.size（）;int i,j,k,l,m,n;Mat B,G;B.resize（X_No）;G.resize（Y_No）;Mat U_0 = U;for（k=1;k+）/增加迭代次数的限制，并设置误差的分析Mat C,C_1;Mat U_1;/U的过渡矩阵,B的转制称U_0,阶数为k+1乘Y_NoU_1.resize（k+1）;Mat B_1,G_1;C.resize（k+1）;C_1.resize（k+1）;B_1.resize（k+1）;G_1.resize（k+1）;for（i=0;k+1;Ci.resize（k+1）;C_1i.resize（k+1）;U_1i.resize（Y_No）;B_1i