多元多项式的公式如图，现在有n个数据（x,y,f(x,y)）来拟合这个多项式，用最小二乘法，具体程序如下面的代码
Z=xlsread('Linton.xlsx');
x=Z(:,1);%x的数据，为风速，范围是0到40
y=Z(:,2);%y的数据，为风向，范围是0到360
z=Z(:,3);%频率数据
n = 7; % n为多项式的阶数，参数一共为(n+1)²个
p = length(x); q = length(y);
[X,Y,i1,i2] = ndgrid(x,y,0:n,0:n);
F = X.^i1.*Y.^i2;
F = reshape(F, p*q, (n+1)^2); %这一段是最小二乘法的左边矩阵，F矩阵具体形式如下图
A = F'*F;
b = F'*z;
a = A\b;%最小二乘法
然后现在就是得加入2个约束条件，一个是x在0到40，y在0到360，f(x,y)的积分等于1，然后f(x,y）永远大于0这两个约束条件。积分为1这个约束可以用拉格朗日乘子法来写。如果是下图的普通多项式，可以加下面一段代码解决积分为1
G = 1./(n+1:-1:1).*40.^(n+1:-1:1);% 这段是上图的普通多项式积分为1，直接写出了积分形式然后用拉格朗日乘子法
%现在需要对最上面的多元多项式积分为1，需要改成x在0到40，y在0到360的二重积分为1
aa = [A, F'; F, 0]\[b;1];%拉格朗日乘子法
现在是得把G矩阵写成二重积分，想问下这种多元多项式进行二重积分的代码应该怎么写。再就是f(x,y)永远大于0的约束要如何加上，以前只写过一维数据，然后普通多项式MATLAB也有代码，这个多元多项式得自己编写
