多项式函数对所给的坐标进行拟合:构造关于系数a j 的多元函数:n2s( ai 1，L，apq)g[f(xg，yg )Z g]g 1点（311，…，a pq ）是多元函数s （a 11,L ,a pq ）的极小点，其中 g 为权函数，默认为1,所以点（811，…，a pq ）必须满足方程组s3ijf(x,y)i 1 j 3j x yij 1,11ai i 1 j 1i 1 j 1 iX yf (x, y)a 11a 12y 2a 13yL q 1a 21x a 22xy2a 23xyLq 1a 2q xyM1i 1 i 12 Li 1 q 1a i1xy33Xya iqX y qMp 1p 1p 12Lp 1 qa p1X a p2X y a p3X ya pq X yp,qpq即1x 2x x M xp,yy2y M y q,A a 12La 1qa 22L a 2qM OMa p2La pqa iia 21Ma p1则函数又可表示为 f (x, y)x TAy ,拟合的目标就是求出系数矩阵 A 。
给定一组坐标(x g ’y g ’Z g ) , g 1,2,…,n ,表示有 n 个点。
要求用以下二元p q / i 1 j 1g ( a j X yi 1 j 1zg)2在g 1的情况下，有2[f (X g ,y g ) Z g ]g i2[f (X g ,y g )i 1 j 1Z g ]X g ygg in2g i因此可得nni 1 j 1 i 1 j 1 X g y gf(X g ,y g )X g y g Zgg 1g 1np qni X g1y g1 1 a X gy g 1i 1 jX gY g1Z gg 11 1g 1np,q ni X g1y g1 1 1aX g y gi 1 j 1X gy gz gg 11,1g 1p,qnnai 1 j(x gy g X g y g1、i 1Xy g1z g1,1g 1g 1p,qa u (i, j) v(i, j) (i, j)(1,1),…,(p,q)1,1上式实际共有p q 个等式，可将这比1(1,1) LU pq (1,1) anMO M MUn(p,q) L U pq (p,q) a pq也就是U*a=V 的形式，其中Un(1,1) L U pq (1,1)UM O MUn(P,q)LU pq (p,q)p q 个等式写成矩阵的形式有: v(1,1) M v( p, q)anv(1,1) a M ，V Ma pqv(p,q)2[f(X g ,y g ) 叩石If")]a ija ij x g 1y g 1f (xg ,y g ) x g+g zu (i, j)n(X g g 11yg1 i 1 X g y g 1), v(i,j)ni 1 j 1X g y gZ gg 1U为pq pq阶矩阵，实现函数为function A=leftmatrix（x,p,y,q）；V为长pq的列向量，实现函数为function B=rightmatrix（x,p,y,q,z）。
这样就可以算出列矩阵a, 然后转化成A。
某地区有一煤矿，为估计其储量以便于开采，先在该地区进行勘探。
假设该地区是一长方形区域，长为4公里，宽为5公里。
经勘探得到如下数据: 煤矿勘探数据表请你估计出此地区内（2x4, 1 y 5 ）煤的储量，单位用立方米表示,并用电脑画出该煤矿的三维图象。
如果直接画出三维曲面图形:clear;x=1:4;y=1:5;[X,Y]=meshgrid(x,y)Z=[13.72 25.80 8.47 25.27 22.32;15.47 21.33 14.49 24.83 26.19;23.28 26.48 29.14 12.04 14.58;19.95 23.73 15.35 18.01 16.29]'surf(X ,Y,Z);X =1 2 3 41 2 3 41 2 3 41 2 3 41 2 3 4Y =1 1 1 12 2 2 23 3 3 34 4 4 45 5 5 5Z =13.7200 15.4700 23.2800 19.950025.8000 21.3300 26.4800 23.73008.4700 14.4900 29.1400 15.350025.2700 24.8300 12.0400 18.010022.3200 26.1900 14.5800 16.2900粗略计算体积：底面积乘以平均高度。
p=sum(Z);q=P(：,[2,3,4]);h=sum(q')/15v=2000*4000*hh =20.0773v =1.6062e+008进行线性插值：xi=li nspace(1,4,31);yi=li nspace(1,5,41);[XI,YI]=meshgrid(xi,yi); ZI=i nterp2(X,Y,Z,XI,YI,'li near');surf(XI,YI,ZI);进行三次多项式插值：xi=li nspace(1,4, 31);yi=li nspace(1,5,Zl=interp2(X,Y,Z,XI,YI,' cubic ');surf(XI,YI,ZI);进行插值后计算体积：底面积乘以平均高度1 130、41);[XI,YI]=meshgrid(xi,yi);xi=li nspace(1,4,61);yi=linspace(1,5,81);[XI,YI]=meshgrid(xi,yi);Zl=i nterp2(X,Y,Z,XI,YI,'cubic');surf(XI,YI,ZI);H=0; n=0;for j=21:61for i=1:81H=H+ZI(i,j); n=n+1;endendnH=H/nS=2000*4000;V=S*H n =3321H =20.8222V =1.6658e+00830.上面是插值的方法解题，下面用拟合的方法解题。
为此编写了几个M函数：fun cti on U=leftmatrix(x,p,y,q)% U*a=V a 为系数列矩阵，长度为 p*q% U 为左边p*q 乘p*q 矩阵 % x,y 为长度一致的列矩阵，给定点的坐标 % p,q为拟合的函数中x,y 的幕的最高次数m=le ngth(x);if (nargin~=4) & (m~=le ngth(y))rightmatrix.mfun cti on V=rightmatrix(x,p,y,q,z)% U*a=V % V为一个列向量长为p*q% x y z 为点的坐标%p q 分别为xy 幕的最高次数 if n arg in~=5error( 'error check check! rightmatrix' endV=zeros(p*q,1); for i=1 : p*qx_z=quotie nt(i-1,q); y_z=mod(i-1,q);V(i,1)=qiuhe(x,x_z,y,y_z,z); endquuotie nt.mfun cti on sh=quotie nt(x,y) % sh 为x/y 的商sh=(x-mod(x,y))/y;error(end 'error check check!' );U_le ngth=p*q; U=zeros(U_le ngth,U_le ngth); for i=1 : p*q for j= 1 : p*qx_z=quotie nt(j-1,q)+quotie nt(i-1,q); y_z=mod(j_1,q)+mod(i_1,q); U(i,j)=qiuhe(x,x_z,y,y_z); end end% U 为p*q 阶方阵 %赋值0，目的是分配内存% x 的幕的次数，quotie nt %y 的幕的次数为求商qiuhe.mfun cti on he=qiuhe(x,p,y,q,z)% he x A p*y A q 从1 —>m 的和% x,y 向量长度相同% p,q分别为x,y的幕的次数m=le ngth(x);if (nargin<4 )&(m~=length(y))error( 'error check check!' );endif n arg in==4 %没有z , z=on es(m,1);end he=0;for i=1:mhe=he+x(i)Ap * y(i)Aq*z(i);% 1-->m 求和end先输入x=._、=••• z=. p=q=(注意x,y,z是向量); 拟合得到系数a,也就是得到了拟合的函数；根据拟合函数计算给定点(xx, yy)的函数值zz=f(xx, yy)并进行画图检验程序保存于M文件fit.m。
clear;[X,Y]=meshgrid(1:4,1:5);Z=[13.72 25.80 8.47 25.27 22.32;15.47 21.33 14.49 24.83 26.19;23.28 26.48 29.14 12.04 14.58;19.95 23.73 15.35 18.01 16.29]';x=reshape(X,20,1);y=reshape(Y,20,1);z=reshape(Z,20,1);p=4；q=5;U=leftmatrix(x,p,y,q); % U*a_ n=VV=rightmatrix(x,p,y,q,z);%a_n=i nv(U)*V;a_n=U\V;for i=1 : length(a_n) % 把长为p*q 的列向量a_n转换成p*q的矩阵aa ii=quotie nt(i-1,q)+1; % quotie nt %输入量至少为四，x,y行向量长度必需一样默认为元素全部为1的向量求商jj=mod(i-1,q)+1;aa(ii,jj)=a_n(i,1);endaa m=31;n=41; %m=4;n=5; [XI,YI]=meshgrid(linspace(1,4,m),linspace(1,5,n)); xx=reshape(XI,m*n,1); yy=reshape(YI,m*n,1);zz=zeros(m*n,1);xy=zeros(m*n,1);xt=zeros(m*n,1);yt=zeros(m*n,1);%zz=0; % zz for i=1 : p % for j=1 : q % aa xt=xx.A(i-1); yt=yy.A(j-1);xy=xt.*yt; zz=zz+aa(i,j).*xy;endendZI=reshape(zz,n,m); surf(XI,YI,ZI); %axis([1 4 1 5 0 30]) aa = 1.0e+003 *0.1465 -0.2678 0.2132 -0.0624 0.0058-0.7287 1.3972 -0.9275 0.2412 -0.02100.4416 -0.8415 0.5487 -0.1407 0.0122-0.0680 0.1295 -0.0839 0.0214 -0.0018是 xx,yy 代入所 拟 合的函数求出的函数 值函数为艺 aa(i,j)*x A i*y A j, (i=1...p,j=1...q) 为 pxq 的系数的矩 阵权函数在拟合的函数非常重要，不过她只能按你遇到的具体问题来取，我这里为1当p , q 越大时，拟合的函数与原数据的方差越小，但是有可能函数 本身抖动非常厉害，可以画图看出来。