北航数值分析大作业三
《数值分析》大作业（3）一、 算法设计方案1．算法的总体设计根据题目中给定的关于 z、t、u 的数表及方程组，可以得到关于 z、x、y 的数表，以此数表作为基础，可以得到 Z(x，y)的表达式要在区域 D = {( x, y ) 0 £ x £ 0.8, 0.5 £ y £ 1.5} 上作二元拟合函数 krsp ( x, y ) = å c rs x yr , s =010 20，并使 p ( x, y ) 在某一最小的 k 值达到如下精度，2s = åå ( f ( x , y ) - p ( x , y ))£ 10-7 ，其中 x = 0.08i ， y = 0.5 + 0.05 j i j i j i ji =0 j =0拟合节点的 xi ， y j 值已知，要求出相应的 z ij = g (t ij , u ij ) =f ( xi, y j ) 先解非线性方程组由一组 ( xi, y j ) 可得到一组 (t ij , u ij ) ，再根据所给的数表作二元分片插值得到函数 z = g (t, u ) ，最后将 (t ij , u ij ) 代入函数 z = g (t, u ) 得到z ij = g (t ij, u ij ) = f ( xi, y j ) 。
用到的算法有：求解非线性方程组，二元二次分片插值，二元函数拟合2．算法的具体实现（1）解非线性方程组题中的非线性方程组采用 Newton 法是收敛的，且有较快的收敛速度 在 x，y 已知的情况下，方程组是关于 t，u，v，w 的非线性方程组求解时近似解向量 x( k ) 满足条件x( k ) - x( k -1)x( k ) ¥ ¥£ 10-12 即对向量形式为ì f1 ( x1 , x2 ,…,xn ) = 0F ( x) = 0 ，一般形式为 …, = 的非线性方程组，ïï f2 ( x1 , x2 , xn ) 0íï Mïî fn ( x1 , x2 ,…,xn ) = 0令 ( k ) ( k ) ( k ) ( k ) ( k ) h = (h1, h2,L, hn), hj¹ 0, j = 1, 2,L n 1 1 1 1 f ( x( k ) + h ( k ) e ) - f ( x( k ) )h( k ) 1f ( x( k ) + h ( k ) e ) - f ( x( k ) )h 1 n 1 1 L ( k ) nJ ( x( k ) , h( k ) ) = M Mf ( x( k ) + h ( k ) e ) - f ( x( k ) )f ( x( k ) + h ( k ) e ) - f( x( k ) ) n 1 1 1 h( k ) 1 n n 1 n hL( k ) n其中 e j 是第 j 个 n 维基本单位向量。
迭代公式x( k +1) = x( k ) - J ( x( k ) , h( k ) )-1 F ( x( k ) ), k = 0,1,……（2）二元二次分片插值根据数表作二元二次分片插值得到函数 z = g (t, u ) 分片插值能保证收 敛性，且有较高的精度根据给定的数表，可将整个插值区域分成 16 个小 的区域分片插值最终要求 z ij = g (t ij, u ij ) = f ( xi, y j ) ，故先判断(t ij, u ij ) 所在的区域，再作此区域的插值，计算 z ij 3）二元函数拟合拟合节点为 ( x , y , z), i = 0,1L10, j = 0,1L20 基函数为乘积型的基函i j ijr s数{j ( x)y( y)(r = 0...k; s = 0...k )} ，其中j ( x) = xr ,y( y) = y s 拟合函数为r sk kp( x, y) = åå crsjr ( x)y s ( y) ，其中C = crs 为要求的系数矩阵实际上就是先s =0 r =0求解系数 C 的矩阵：C = (BT B)-1 BTUG(GT G)-1其中：k kæ1 x0L x0 ö æ1y0 Ly0 öB = ç1x1 Lx1 ÷,G = ç1y1 Ly1 ÷,U = ( f ( x , y ))ç M M O M÷ ç M M O M ÷i j n´mç ÷ ç ÷1 x L xk1 y L y kè n n ø èn n ø令 A = (BT B)-1 BTU， DT = G(GT G)-1，作变换得 BT BA = BTU ，GT GD = GT ，可通过列主元素 Gauss 法解得矩阵 A 和 D，再算得C = ADT 。
计算时的 k 值需要程序自动确定，要求最小的 k 值使精度达到：10 20i i i is =åå ( f ( x , y ) - p( x , y ))2 £ 10-7 i=0 j = 0二、计算结果1 ． 二 元 二 次 分 片 插 值 得 到 数 表 ：(i = 0,1,L,10; j = 0,1,L, 20) xi ， y j ，f ( xi, y j ) ，2．二元拟合选择过程的 k ，s 值3 ． 达 到 精 度 要 求 的 k ， s 的 值 及c rs ( r = 0,1,L, k; s = 0,1,L, k ) p ( x, y )中 的 系 数4．数表： * ， y * ， f ( x * , y * ) ， p ( x * , y * ) ，(i = 1, 2,L, 8; j = 1, 2,L, 5) xi j i j i j附录 全部源程序#include
/void main(){double u_0[6]={0,0.4,0.8,1.2,1.6,2}; double t_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}};int i,j,k;double u_1[11][21], t_1[11][21],z[11][21],x_0[4]={0.5,0.5,0.5,1};double delt[11],c[11][11]; double f_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.2f yi=%3.2f f(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=%d delt=%14.12e",i,delt[i-1]);printf("\n");}printf("\n");printf("\t"); printf("k=5 delt=%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 函数：求无穷范数-----------------------------------------------------------------------*/double fabsmax(double x[4]){double b=0; int i; b=fabs(x[0]); for(i=1;i<=3;i++){if(fabs(x[i])>b)b=fabs(x[i]);}return b;}/*----------------------------------------------------------------------- fabsmax_flag 函 数 ： 返 回 数 组 x[4] 中 最 大 元 素 所 在 的 位 置-----------------------------------------------------------------------*/int fabsmax_flag(double x[4],int k){double b=0; int i,flag=k; b=fabs(x[k]); for(i=k+1;i<=3;i++)if(fabs(x[i])>b){flag=i;b=fabs(x[i]);}return flag;}/*----------------------------------------------------------------------- fabsmax_flag_1 函数：比较数据绝对值的大小，并返回按模最大值所在的行号-----------------------------------------------------------------------*/int fabsmax_flag_1(double x[11],int k,int l){double b=0; int i,flag=k; b=fabs(x[k]); for(i=k+1;i<=l;i++)if(fabs(x[i])>b){flag=i;b=fabs(x[i]);}return flag;}/*----------------------------------------------------------------------- gauss 函 数 ： 列 主 元 素 Gauss 法 解 线 性 方 程 组-----------------------------------------------------------------------*/void gauss(double a_0[4][4],double b[4],double deltx[4]){double a[4][4],x[4],aa; int i,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=fabsmax_flag(x,k); if(flag!=k){for(j=k;j<=3;j++){aa=a[k][j]; a[k][j]=a[flag][j]; a[flag][j]=aa;}aa=b[k]; b[k]=b[flag]; b[flag]=aa;}for(i=k+1;i<=3;i++){aa=a[i][k]/a[k][k];for(j=k+1;j<=3;j++)a[i][j]-=aa*a[k][j];b[i]-=aa*b[k];}}deltx[3]=b[3]/a[3][3];for(k=2;k>=0;k--){for(j=k+1;j<=3;j++)b[k]-=a[k][j]*deltx[j];deltx[k]=b[k]/a[k][k];}}/*----------------------------------------------------------------------- function 函 数 ： Newton 法 解 非 线 性 方 程 组-----------------------------------------------------------------------*/void function(double t_1[11][21],double u_1[11][21],double x_0[4]){double a[4][4]={{0,1,1,1},{1,0,1,1},{0.5,1,0,1},{1,0.5,1,0}};double b[4],deltx[4],x[4]; int i,j,k; for(i=0;i<=10;i++)for(j=0;j<=20;j++){for(k=0;k<=3;k++)x[k]=x_0[k];for(;;){a[0][0]=-0.5*sin(x[0]); a[1][1]=0.5*cos(x[1]); a[2][2]=-sin(x[2]); a[3][3]=cos(x[3]);b[0]=-(0.5*cos(x[0])+x[1]+x[2]+x[3]-0.08*i-2.67);b[1]=-(x[0]+0.5*sin(x[1])+x[2]+x[3]-(0.5+0.05*j)-1.07);b[2]=-(0.5*x[0]+x[1]+cos(x[2])+x[3]-0.08*i-3.74);b[3]=-(x[0]+0.5*x[1]+x[2]+sin(x[3])-(0.5+0.05*j)-0.79);gauss(a,b,deltx);for(k=0;k<=3;k++)x[k]+=deltx[k];if((fabsmax(deltx)/fabsmax(x))<=1e-12){t_1[i][j]=x[0]; u_1[i][j]=x[1]; break;}}}}/*-----------------------------------------------------------------------f 函数：作二元分片插值时用到的函数，用于确定给定的（t，u）所在的插值区域-----------------------------------------------------------------------*/void f(double t_1,double u_1,int flag[2]){int flag1,flag2; flag1=(int)(t_1/0.1); flag2=(int)(u_1/0.2); if(flag1>=0&&flag1<3)flag1=0;else if(flag1>=3&&flag1<5)flag1=1;else if(flag1>=5&&flag1<7)flag1=2;elseflag1=3;if(flag2>=0&&flag2<3)flag2=0;else if(flag2>=3&&flag2<5)flag2=1;else if(flag2>=5&&flag2<7)flag2=2;elseflag2=3; flag[0]=flag1; flag[1]=flag2;}/*----------------------------------------------------------------------- chazhi 函 数 ： 分 片 二 元 插 值-----------------------------------------------------------------------*/ void chazhi(double t_1[11][21],double u_1[11][21],double z[11][21],double t_0[6],double u_0[6],double z_0[6][6]){double l_t[3],l_u[3]; int i,j,k,l,flag[2]; for(i=0;i<=10;i++)for(j=0;j<=20;j++){z[i][j]=0;f(t_1[i][j],u_1[i][j],flag);l_t[0]=(t_1[i][j]-t_0[flag[0]+1])*(t_1[i][j]-t_0[flag[0]+2])/0.08; l_t[1]=(t_1[i][j]-t_0[flag[0]])*(t_1[i][j]-t_0[flag[0]+2])/(-0.04); l_t[2]=(t_1[i][j]-t_0[flag[0]])*(t_1[i][j]-t_0[flag[0]+1])/0.08; l_u[0]=(u_1[i][j]-u_0[flag[1]+1])*(u_1[i][j]-u_0[flag[1]+2])/0.32; l_u[1]=(u_1[i][j]-u_0[flag[1]])*(u_1[i][j]-u_0[flag[1]+2])/(-0.16); l_u[2]=(u_1[i][j]-u_0[flag[1]])*(u_1[i][j]-u_0[flag[1]+1])/0.32;for(k=0;k<=2;k++)for(l=0;l<=2;l++)z[i][j]+=z_0[flag[0]+k][flag[1]+l]*l_t[k]*l_u[l];}}/*----------------------------------------------------------------------- gauss_1 函数：列主元素 Gauss 法解线性方程组-----------------------------------------------------------------------*/void gauss_1(double b[11][11],double u[11][21],double a[11][21],int l){double b_1[11][11],b_2[11][11],b_3[11][21],x[11],aa;int i,j,k,flag;for(i=0;i<=10;i++)for(j=0;j<=l;j++)b_2[j][i]=b[i][j];for(i=0;i<=l;i++)for(j=0;j<=20;j++){b_3[i][j]=0; for(k=0;k<=10;k++) b_3[i][j]+=b_2[i][k]*u[k][j];}for(i=0;i<=l;i++)for(j=0;j<=l;j++){b_1[i][j]=0; for(k=0;k<=10;k++) b_1[i][j]+=b_2[i][k]*b[k][j];}for(k=0;k
