求教 Matlab 三维数值积分

[cross]

目前我采用如下方法计算三维积分, 但是出现积分结果随着剖分单元数量增加先下降收敛, 然后反弹上升的曲线. 大家能否帮我指出思路中的错误, 先谢谢了.
1. 数据结构: 50万个离散数据点, 在一圆柱体空间非均匀分布, 对于每个点, 有坐标值(x,y,z)和其上对应函数值, 例如 psi(x,y,z).
2. 对全部积分空间按照坐标划分立方体单元, 寻找包含在每个单元内的数据点, 及其单元中心坐标.
3. 若单元内点数少于10, 以此单元外壁为起点, 步长(长,或宽,或高的0.1倍)固定, 在相邻单元内搜索临近点, 找到的邻点大于10,停止.
4. 若单元内点数(#)大于100, 以ceil(#/100)为步长, 抽取少于等于100个数据点.
5. 用每个单元内数据点线性回归(regress) 表达式f=a0+a1*x+a2*y+a3*z+a4*x^2+a5*y^2+a6*z^2+a7*x*y+a8*x*z+a9*y*z 的系数a0到a9.
6. 在每单元内对f三次积分(triplequad). 并对全部单元的积分值求和.
下面是我使用的代码:
clear
load data_2s
tic
data=ones(474382,1); %%%%% for test, must get only the volume
x1=coord(:,1); y1=coord(:,2); z1=coord(:,3);
nx= 12 ; %element number along x-axis
ny= 10 ; %element number along y-axis
nz= 12 ; %element number along z-axis
step_x=(max(x1)-min(x1))/nx ; % set the scale of each rectangular element
step_y=(max(y1)-min(y1))/ny ;
step_z=(max(z1)-min(z1))/nz ;
for i=1:nx
for j=1:ny
for k=1:nz
disp([i,j,k])
search=find(((x1>=(min(x1)+(i-1)*step_x))&(x1<(min(x1)+i*step_x)))&((y1>=(min(y1)+(j-1)*step_y))&(y1<(min(y1)+j*step_y)))&((z1>=(min(z1)+(k-1)*step_z))&(z1<(min(z1)+k*step_z))));
element_center=[(min(x1)+(i-1/2)*step_x),(min(y1)+(j-1/2)*step_y),(min(z1)+(k-1)*step_z)];
kk=5; % initial value
while (length(search)<10)
kk=kk+1;
search=find((abs(x1-element_center(1))<(step_x/10)*kk)& abs(y1-element_center(2))<(step_y/10)*kk & (abs(z1-element_center(3))<(step_z/10)*kk));
end
r_size2(i,j,k)=length(search);
while (length(search)>100)
pick=(1:ceil(length(search)/100):length(search)); % shrink # of search to <100
search=search(pick);
end
clear pick element_center
if (length(search)>=10)
if length(find(data(search)==0))==length(search)
Q(i,j,k)=0;
else
n=[ones(length(search),1),x1(search),y1(search),z1(search),x1(search).*x1(search),y1(search).*y1(search),z1(search).*z1(search),x1(search).*y1(search),x1(search).*z1(search),y1(search).*z1(search)];
[b,bint,r,rint,s]=regress(data(search),n);
clear n
fun=strcat(num2str(b(1)),'+',num2str(b(2)),'*x+',num2str(b(3)),'*y+',num2str(b(4)),'*z+',num2str(b(5)),'*x.*x+',num2str(b(6)),'*y.*y+',num2str(b(7)),'*z.*z+',num2str(b(8)),'*x.*y+',num2str(b(9)),'*x.*z+',num2str(b(10)),'*y.*z'); % modified 08.11.03
Q(i,j,k)=triplequad(inline(fun),min(x1)+(i-1)*step_x,min(x1)+i*step_x,min(y1)+(j-1)*step_y,min(y1)+j*step_y,min(z1)+(k-1)*step_z,min(z1)+k*step_z);
end
end

end
end
end
int=sum(reshape(Q,1,[]))
time=toc