1、，若，则选；若，等上述方法一样。其中：4.二元函数拟合已通过求解非线性方程组和插值得到的1121的（x,y）,z）二元矩形数表为基准，进行给定不同次数的二元函数最小二乘法拟合。拟合过程选取正交函数系作为基函数，即拟合函数可表示为：，其中和分别为次数为r和s的正交多项式系。将其写为矩阵形式为： 同理：所以：其中利用x和y的正交函数系同其幂函数系数阵以及（2）中所求系数矩阵，最终得到所求系数矩阵，即：该子程序中的矩阵乘积运算调用了通用的矩阵相乘的子程序。二、全部源程序#include stdio.hmath.hint p1,q1;/全局变量/二元函数插值/x为插值表的x坐标点，y为插值表的y坐标点
2、，z为插值表的函数值/u为待插值点的x坐标，v为为待插值点的y坐标double cz（double x,double y,double z,double u,double v） double h,f,hh,w; double bb10; int i,j,ii,jj,k,k1,k2,kk; h=0.4;f=0.2;hh=0.0; if（ux4-h/2） k1=4; else for（i=2;i（xi-h/2）&u=（xi+h/2） k1=i; if（vy4-f/2） k2=4; for（j=2;j（yj-f/2）&v=（yj+f/2） k2=j; for（ii=k1-1;iik1+2;ii+）
3、bbii=0.0; for（jj=k2-1;jjk2+2;jj+） k=ii*6+jj;hh=zk; for（kk=k2-1;kkfan） fan=fabs（arrayi）; return fan;/矩阵转置/a为待转置矩阵，维数为m*n，转置后存放在矩阵b中zhuanzhi（double a,double b,int m,int n） int i,j;m; for（j=0; bm*j+i=ai*n+j; /矩阵相乘/a矩阵维数为m*n，b矩阵维数为n*k，乘积结果存放在c矩阵中 void brmul（a,b,m,n,k,c） int m,n,k; double a,b,c; int i,j,
4、l,u; for （i=0; i=m-1; i+） for （j=0; j=k-1; j+） u=i*k+j; cu=0.0; for （l=0; l=n-1; l+） cu=cu+ai*n+l*bl*k+j;/对矩阵部分求逆/a为待求逆矩阵，维数为n1*n1，对其中n*n进行操作int brinv（double A,int n,int n1） double *a,d,p; int *is,*js,i,j,k,l,u,v; a=malloc（n*n*sizeof（double）; is=malloc（n*sizeof（int）; js=malloc（n*sizeof（int）; *（a+i*n
5、+j）=*（A+i*n1+j）; for （k=0; kd） d=p; isk=i; jsk=j; if （d+1.0=1.0） free（is）; free（js）; free（a）; printf（err*not invn）; return（0）; if （isk!=k） u=k*n+j; v=isk*n+j; p=au; au=av; av=p; if （jsk! u=i*n+k; v=i*n+jsk; l=k*n+k; al=1.0/al; if （j! au=au*al; if （i! u=i*n+j; au=au-ai*n+k*ak*n+j; au=-au*al; for （k=n
6、-1; k=0; k-） if （jsk! v=jsk*n+j; v=i*n+isk; *（A+i*n1+j）=*（a+i*n+j）; free（is）; free（a）; return（1）;/曲面拟合/m为正交多项式系最高次值，array为拟合系数矩阵/（u,v）为待拟和点double qmlh（int m,double u,double v,double array） int r,s; double sum=0.0; for（s=0;ss+） for（r=0;rr+） sum+=arrayr*m+s*pow（u,r）*pow（v,s）; return sum;/Gauss列主元素消去法解
7、线性方程组/a,系数矩阵; b，值向量，返回时存放解向量 int agaus（a,b,n） int n; double a,b; int *js,l,k,i,j,is,p,q; double d,t; l=1;kd） d=t; is=i; if （d+1.0=1.0） l=0; p=i*n+k; q=i*n+jsk; t=ap; ap=aq; aq=t; if （is! for （j=k; p=k*n+j; q=is*n+j; t=bk; bk=bis; bis=t; if （l=0） free（js）; printf（failn d=ak*n+k; for （j=k+1; ap=ap/d;
8、bk=bk/d; for （i=k+1; for （j=k+1; p=i*n+j; ap=ap-ai*n+k*ak*n+j; bi=bi-ai*n+k*bk; d=a（n-1）*n+n-1; if （fabs（d）+1.0=1.0） bn-1=bn-1/d; for （i=n-2;ii-） t=0.0; for （j=i+1; t=t+ai*n+j*bj; bi=bi-t; jsn-1=n-1;kk-） t=bk; bk=bjsk; bjsk=t;main（） int i,j,r,s,b1,b2,k,ss,kk; double xx11=0.,yy21=0.,x6=0.,y6=0.,xin8=
9、0.,yin5=0.; double a44=0.,b4=0.,h4=0.,H4=0.,g4=0.; double C1111=0.,B1111=0.,BZ1111=0.,BT1111=0.;double BTT1111=0.,BU1121=0.; double G2111=0.,GZ1121=0.,GT1111=0.;double GTT2111=0.,U1121=0.; double bk,eps,u,v,uu,vv,w,delta,t1,t2,sigama,lh,llh; double X4=0.5,1.5,-0.5,1.5;/初始迭代向量 double z66=-0.5,-0.42,-
10、0.18,0.22,0.78,1.5,-0.34,-0.5,-0.5,-0.34,-0.02,0.46,0.14,-0.26,-0.5,-0.58,-0.5,-0.26,0.94,0.3,-0.18,-0.5,-0.66,-0.66,2.06,1.18,0.46,-0.1,-0.5,-0.74,3.5,2.38,1.42,0.62,-0.02,-0.5;/插值表 FILE *fp; if（fp=fopen（SY0604229.txt,w+）=NULL）/建立一个新文件存放输出结果non bk=1.0;eps=1e-12;delta=10e-7;t1=0.;t2=0.;sigama=0.;lh=
11、0.; for（p1=0;p111;p1+） xxp1=0.08*p1; for（q1=0;q1=eps）;/达到精度退出循环 u=X1;v=X0;/得出插值坐标 w=cz（x,y,z,u,v）;/对给定点进行插值计算 Up1q1=w; fprintf（fp,x=%13.12e y=%13.12e f（x,y）=%13.12en,xxp1,yyq1,w）;/将插值结果输出到文件SY0604229.txtnn fprintf（fp,nnn do for（k=1; for（j=0;=k; for（i=0; xi=0.08*i;/待拟合点的x坐标向量 Bij=pow（xi,j）;/x的各阶正交函数的
12、幂次矩阵 for（b2=0;b2b2+） for（b1=0;b1=1e-7）;%13.12e ,Cij）；/将拟合系数矩阵输出到文件SY0604229.txtnLH=%13.12e sigama=%13.12en,lh,sigama）;/将误差输出到文件SY0604229.txt6;C%d%d=%13.12en,i,j,Cij）; bk=0.0;8; xinp1=0.1*p1+0.1;/初始化x*向量5;q1+） yinq1=0.7+0.2*q1;/初始化y*向量i+） xi=0.4*i;/初始化插值节点的x坐标j+） yj=0.2*j;/初始化插值节点的y坐标 dnewta2（X,b,xin,yin）;v=X
