北航数值分析A大作业3

}while(k

printf(\法在该初值不收敛\\n\);}

//查找n维向量V中与常数a最接近的元素的下标int Near_Index(double* v,double a,int n){

int i,Index;double min;

min=abs(v[0]-a);Index=0;

for(i=1;i

if(min>abs(v[i]-a)){

min=abs(v[i]-a);Index=i;}}

return Index;}

//求数表z(t,u)在点（x，y）处的分片二次代数差值

double ChaZhi(double a,double b){

double z[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};

double x[6]={0,0.2,0.4,0.6,0.8,1};double y[6]={0,0.4,0.8,1.2,1.6,2};

double L[3],L_[3],Z;int i,j,k,r,t;

i=Near_Index(x,a,6);j=Near_Index(y,b,6);

//插值区域边界处插值节点的选取

11

if(i==0)i=1;if(i==5)i=4;if(j==0)j=1;if(j==5)j=4;

for(k=i-1;k<=i+1;k++){

L[k-i+1]=1;

for(t=i-1;t<=i+1;t++){

if(t!=k)L[k-i+1]*=((a-x[t])/(x[k]-x[t]));}}

for(r=j-1;r<=j+1;r++){

L_[r-j+1]=1;

for(t=j-1;t<=j+1;t++){

if(t!=r)L_[r-j+1]*=((b-y[t])/(y[r]-y[t]));}}

Z=0;

for(k=i-1;k<=i+1;k++){

for(r=j-1;r<=j+1;r++){

Z+=(L[k-i+1]*L_[r-j+1]*z[k][r]);}}

return Z;}

//对m×n数表U(x,y)进行二元多项式拟合

void NiHe(double*U,double*x,double*y,int m,int n,int p,int q,double*C)//U为拟合数据点处的函数值{

int i,j;

double *B,*BT,*BTB,*G,*GT,*GTG,*BTUG,*temp,*tempT,*CT;B=(double*) calloc(m*p,sizeof(double));BT=(double*) calloc(p*m,sizeof(double));BTB=(double*) calloc(p*p,sizeof(double));

12

G=(double*) calloc(n*q,sizeof(double));GT=(double*) calloc(q*n,sizeof(double));GTG=(double*) calloc(q*q,sizeof(double));BTUG=(double*) calloc(p*q,sizeof(double));temp=(double*) calloc(p*n,sizeof(double));for(i=0;i

for(j=0;j

}

for(i=0;i

for(j=0;j

}

Transpose(B,m,p,BT);Transpose(G,n,q,GT);

Array_Mult_Array(BT,B,p,m,p,BTB);Array_Mult_Array(GT,G,q,n,q,GTG);Array_Mult_Array(BT,U,p,m,n,temp);Array_Mult_Array(temp,G,p,n,q,BTUG);free(B);free(BT);free(G);free(GT);free(temp);

temp=(double*) calloc(p*q,sizeof(double));Solve_lin(BTB,p,BTUG,temp,q);free(BTB);free(BTUG);

tempT=(double*) calloc(q*p,sizeof(double));CT=(double*) calloc(q*p,sizeof(double));Transpose(temp,p,q,tempT);Solve_lin(GTG,q,tempT,CT,p);Transpose(CT,q,p,C);free(temp);

13

free(CT);free(tempT);free(GTG);}

//求多项式拟合函数P(x,y)在点（x，y）处的函数值double Pxy(double *C,int p,int q,double x,double y){

int i,j;

double sum=0;for(i=0;i

for(j=0;j

sum+=(C[i*q+j]*pow(x,i)*pow(y,j));

}

return sum;}

三、输出结果

14

15