北航数值分析大作业三

《数值分析》大作业（3）一、算法设计方案1．算法的总体设计根据题目中给定的关于z、t、u的数表及方程组，可以得到关于z、x、y的数表，以此数表作为基础，可以得到Z(x，y)的表达式。要在区域,00.8,0.51.5Dxyxy上作二元拟合函数,0,ksrrsrspxyycx，并使,pxy在某一最小的k值达到如下精度，21020700,,10iijjijfpyyxx，其中0.08iix，0.50.05jjy。拟合节点的ix，jy值已知，要求出相应的,,ijijijijgfytuxz。先解非线性方程组由一组,ijyx可得到一组,ijijtu，再根据所给的数表作二元分片插值得到函数,zgtu，最后将,ijijtu代入函数,zgtu得到,,ijijijijgfytuxz。用到的算法有：求解非线性方程组，二元二次分片插值，二元函数拟合。2．算法的具体实现（1）解非线性方程组题中的非线性方程组采用Newton法是收敛的，且有较快的收敛速度。在x，y已知的情况下，方程组是关于t，u，v，w的非线性方程组。求解时近似解向量kx满足条件11210kkkxxx。即对向量形式为()0Fx，一般形式为11221212(,,)0(,,)0(,,)0nnnnfxxxfxxxfxxx…,…,…,的非线性方程组，令()()()()()12(,,,),0,1,2,kkkkknjhhhhhjn()()()()()()1111111()()1()()()()()()()()1111()()1()()()()(,)()()()()kkkkkknkknkkkkkkkknnnnkknfxhefxfxhefxhhJxhfxhefxfxhefxhh其中je是第j个n维基本单位向量。迭代公式(1)()()()1()(,)(),0,1,kkkkkxxJxhFxk……（2）二元二次分片插值根据数表作二元二次分片插值得到函数,zgtu。分片插值能保证收敛性，且有较高的精度。根据给定的数表，可将整个插值区域分成16个小的区域。分片插值最终要求,,ijijijijgfytuxz，故先判断,ijijtu所在的区域，再作此区域的插值，计算ijz。（3）二元函数拟合拟合节点为j(,,),0,110,0,120iijxyzij。基函数为乘积型的基函数{()()(0...;0...)}rsxyrksk，其中(),()rsrsxxyy。拟合函数为00(,)()()kkrsrssrpxycxy，其中rsCc为要求的系数矩阵。实际上就是先求解系数C的矩阵：11()()TTTCBBBUGGG其中：000011111111,((,))11kkkkijnmkknnnnxxyyxxyyBGUfxyxxyy　,　　令1()TTABBBU，1()TTDGGG，作变换得TTBBABU，TTGGDG，可通过列主元素Gauss法解得矩阵A和D，再算得TCAD。计算时的k值需要程序自动确定，要求最小的k值使精度达到：102027i=00=((,)(,))10iiiijfxypxy。二、计算结果1．二元二次分片插值得到数表：ix，jy，,ijfyx，0,1,,10;0,1,,20ij。2．二元拟合选择过程的k，值。3．达到精度要求的k，的值及,pxy中的系数0,1,,;0,1,,rsrkskc。4．数表：*ix，*jy，**,ijfyx，**,ijpyx，1,2,,8;1,2,,5ij。附录全部源程序#includeiostream.h#includestdio.h#includemath.hdoublefabsmax(doublex[4]);intfabsmax_flag(doublex[4],intk);voidgauss(doublea_0[4][4],doubleb[4],doubledeltx[4]);voidfunction(doublet_1[11][21],doubleu_1[11][21],doublex_0[4]);voidchazhi(doublet_1[11][21],doubleu_1[11][21],doublez[11][21],doublet_0[6],doubleu_0[6],doublez_0[6][6]);voidf(doublet_1,doubleu_1,intflag[2]);intnihe(doubleu[11][21],doublec[11][11],doubledelt[11]);voidgauss_1(doubleb[11][11],doubleu[11][21],doublea[11][21],intl);voidgauss_2(doubleg[21][21],doubled[21][21],intl);intfabsmax_flag_1(doublex[11],intk,intl);voidf_p(doublef_xy[8][5],doubleg_xy[8][5],doublex_0[4],doubleu_0[6],doublet_0[6],doublez_0[6][6],doublec[11][11]);/*-----------------------------------------------------------------------main函数：在给定区域上作二元拟合得到满足精度要求的P(X,Y)，为确定拟合节点(x,y,z)需要解非线性方程组和作二元分片插值。-----------------------------------------------------------------------*/voidmain(){doubleu_0[6]={0,0.4,0.8,1.2,1.6,2};doublet_0[6]={0,0.2,0.4,0.6,0.8,1.0};doublez_0[6][6]={{-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}};inti,j,k;doubleu_1[11][21],t_1[11][21],z[11][21],x_0[4]={0.5,0.5,0.5,1};doubledelt[11],c[11][11];doublef_xy[8][5],g_xy[8][5];function(t_1,u_1,x_0);chazhi(t_1,u_1,z,t_0,u_0,z_0);k=nihe(z,c,delt);for(i=0;i=10;i++)for(j=0;j=20;j++){printf(\n);printf(\t);printf(xi=%3.2fyi=%3.2ff(xi,yi)=%14.12e,0.08*i,(0.5+0.05*j),z[i][j]);}printf(\n);for(i=1;i=5;i++){printf(\t);printf(k=%ddelt=%14.12e,i,delt[i-1]);printf(\n);}printf(\n);printf(\t);printf(k=5delt=%14.12e,delt[4]);printf(\n);printf(\n);for(i=0;i=k;i++)for(j=0;j=k;j++){printf(\t);printf(c%d%d=%14.12e,i,j,c[i][j]);printf(\n);}printf(\n);f_p(f_xy,g_xy,x_0,u_0,t_0,z_0,c);printf(\n);for(i=1;i=8;i++)for(j=1;j=5;j++){printf(\n);printf(f(%2.1f,%2.1f)=%12.11ep(%2.1f,%2.1f)=%12.11e,0.1*i,(0.5+0.2*j),f_xy[i-1][j-1],0.1*i,(0.5+0.2*j),g_xy[i-1][j-1]);}printf(\n);printf(\n);}/*-----------------------------------------------------------------------fabsmax函数：求无穷范数-----------------------------------------------------------------------*/doublefabsmax(doublex[4]){doubleb=0;inti;b=fabs(x[0]);for(i=1;i=3;i++){if(fabs(x[i])b)b=fabs(x[i]);}returnb;}/*-----------------------------------------------------------------------fabsmax_flag函数：返回数组x[4]中最大元素所在的位置-----------------------------------------------------------------------*/intfabsmax_flag(doublex[4],intk){doubleb=0;inti,flag=k;b=fabs(x[k]);for(i=k+1;i=3;i++)if(fabs(x[i])b){flag=i;b=fabs(x[i]);}returnflag;}/*-----------------------------------------------------------------------fabsmax_flag_1函数：比较数据绝对值的大小，并返回按模最大值所在的行号-----------------------------------------------------------------------*/intfabsmax_flag_1(doublex[11],intk,intl){doubleb=0;inti,flag=k;b=fabs(x[k]);for(i=k+1;i=l;i++)if(fabs(x[i])b){flag=i;b=fabs(x[i]);}returnflag;}/*-----------------------------------------------------------------------gauss函数：列主元素Gauss法解线性方程组-----------------------------------------------------------------------*/voidgauss(doublea_0[4][4],doubleb[4],doubledeltx[4]){doublea[4][4],x[4],aa;inti,j,k,flag;for(i=0;i=3;i++)for(j=0;j=3;j++)a[i][j]=a_0[i][j];for(k=0;k=2;k++){for(i=k;i=3;i++)x[i]=a[i][k];flag=