一、算法设计方案
1、解非线性方程组
将各拟合节点（xi，yj）分别带入非线性方程组，求出与(xi,yi)相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。2、二元二次分偏插值
对数表z(t,u)进行分片二次代数插值，求得对应（tij，uij）处的值，即为f(xi,yj) 的值。根据给定的数表，可将整个插值区域分成 16 个小 的区域，故先判断?t ij, u ij ? 所在，的区域，再作此区域的插值，计算z ij，相应的Lagrange形式的插值多项式为：
p22(t,u)?k?m?1r?n?1??l(t)l?(u)f(t,u)krkrm?1n?1其中
lk(t)?t?tw (k=m-1, m, m+1)?t?tw?m?1kww?km?1l?r(u)?y?yw (r=n-1, n, n+1)?w?n?1yr?yww?rn?13、曲面拟合
从k=1开始逐渐增大k的值，使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，当??10?7时结束计算。拟合基函数φr(x)ψs(y)选择为φr(x)=xr，
?ψs(y)=ys。拟合系数矩阵c通过连续两次解线性方程组求得。C?[crs]，
C?(BTB)?1BTUG(GTG)?1 其中
1
kk?1x0?x0??1y0?y0???k?k?1x?x1y?y11?11?B?[?r(xi)]??，G?[?s(yj)]????????????????k?k????1x10?x10???1y10?y10??U?[f(xi,yj)]4、观察比较
**计算f(xi*,y*j),p(xi,yj)(i?1,2,???,8,j?1,2,???,5)的值并输出结果，以观察
p(x,y)逼近f(x,y)的效果。其中xi*?0.1i,y*j?0.5?0.2j。
二、全部源程序
// hean.cpp : 定义控制台应用程序的入口点。//
#include \#include
void Set_non_JacobiA(double* A,double* x);//求题中非线性方程组对应自变量向量x的雅克比//矩阵
void Set_non_B(double* B,double* x,double a,double b);//求非线性方程组Newton迭代法的右//端式:-F(x)
void Array_Mult_Array(double* A,double* B,int m,int s,int n,double* C);//矩阵相乘AB＝C
void Transpose(double *A,int m,int n,double* AT);//转置void Doolittle(double *A,int n,int *M);
void Solve_LU(double* A,int n,double* b,double* x);//解方程LUx＝b
2
void Solve_lin(double* A,int n,double* B,double* x,int m);//解线性方程组Ax＝B
double Vector_FanShu(double *V,int n);
void Solve_non_Equation(double a,double b,double* x);
//求解非线性方程组
int Near_Index(double* v,double a,int n);
//查找n维向量V中与常数a最接近的元素的下标
double ChaZhi(double a,double b);
//求数表z(t,u)在点（x，y）处的分片二次代数差值
void NiHe(double*U,double*x,double*y,int m,int n,int p,int q,double*C);
//对m×n数表U(x,y)进行二元多项式拟合
double Pxy(double *C,int p,int q,double x,double y);
//求多项式拟合函数P(x,y)在点（x，y）处的函数值
void main(){
double non[4]; double f[11][21];
//非线性方程组变量向量（t,u,v,w）//f(x,y)在拟合节点处的值
double x[11],y[21];//拟合节点double *C;double det=1e-7;double p,d;int i,j,k,t;//设置拟合节点for(i=0;i<11;i++)for(j=0;j<21;j++)
x[i]=0.08*i;y[j]=0.5+0.05*j;//拟合系数矩阵//拟合精度要求
//求拟合节点处的f(x,y)的值for(i=0;i<11;i++)
3
{
for(j=0;j<21;j++){
Solve_non_Equation(x[i],y[j],non);f[i][j]=ChaZhi(non[0],non[1]);}}
printf(\//打印数表f(xi,yi)
printf(\for(t=0;t<21/3;t++){
printf(\---\\n\);
printf(\);for(i=3*t;i<3*t+3;i++){
printf(\y=%4.2f |\,y[i]);}
printf(\);
for(j=0;j<11;j++){
printf(\,x[j]);for(i=3*t;i<3*t+3;i++)
printf(\%+.12e |\,f[j][i]);printf(\);}printf(\---\\n\);
printf(\);}
k=0;printf(\
不同k对应的精度
\);
printf(\);do{
C=(double*) calloc((k+1)*(k+1),sizeof(double));NiHe(*f,x,y,11,21,k+1,k+1,C);
4
数值分析第三次大作业\\n\);
f(xi,yi) \\n\);
//求在当前k值拟合的节点处误差d=0;
for(i=0;i<11;i++){
for(j=0;j<21;j++){
p=Pxy(C,k+1,k+1,x[i],y[j]);d+=((f[i][j]-p)*(f[i][j]-p));}}
printf(\if(d>=det)
free(C);else{
printf(\);printf(\break;
}
}while(++k<11);//打印拟合系数矩阵c
printf(\当k=%d时拟合函数p(x,y)的系数矩阵c：\\n\,k);printf(\);for(t=0;t<(k+1)/3;t++){
for(i=0;i<(k+1);i++){
for(j=3*t;j<3*t+3;j++)
printf(\%+.12e \,C[i*(k+1)+j]);printf(\);}printf(\);}
5
k=%d对应的σ= %.12e\,k,d);
故可知k=k=%d<%.0e，满足题设要求\,det);
北航数值分析A大作业3
