对于3D空间中100k个坐标点及其上面数据值使用matlab进行数值积分, 有没有成熟的程序段哪? 我一直没有找到. 自己采用类似复化梯形法划分坐标网格,对网格内数据求平均计算可以较快解决, 但是精度很低. 如果采用插值, 用interp3, 又会对3维矩阵的大小产生限制, 无法引入足够多的原始数据点, 精度依然无法提高. 对于采用单元内平均值的方法,用一维验证,类似于梯形法,积分sin(x)在-pi到pi之间,结果在0.01和0.001之间,很难说趋近于0.  第三种我采用的办法是划分坐标网格(其实我的模型是柱体)后,对每个单元单元内部点上的数据拟合出 f=a0+a1*x+a2*y+a3*z 在每个单元内求积分解析解,然后对全部单元求和, 这样仍然大量消耗内存. (刚刚经过每次将近一个小时的运算, 发现随着剖分单元增加, 积分结果并不收敛,附上程序草稿)

请教大家有没有什么好办法. 急!!!!!!!!!!! 感谢.

各位如果有时间, 也可以帮忙回答到这个邮箱里面 yaoyuanzhao@gmail.com不胜感激.

[code]

clear

load data_1s(E)R13

tic

total_number=474382;  % total data point number

%%%%%%%% extract %%%%%%%%%%%%%

% with import data structure coord_s=(x,y,z),data_s=(d),coord_p=(x,y,z),data_p=(d) for s and p-state with same size

extract=60;   % extract one data point from each point group with point number=extract

for ii=1:(floor(length(coord)/extract))

coord_new_s(ii,1:3)=coord(1+(ii-1)*extract,1:3);

data_new_s(ii)=data(1+(ii-1)*extract);

end

x1=coord_new_s(:,1); y1=coord_new_s(:,2); z1=coord_new_s(:,3); d1=data_new_s;

clear coord_new_s data_new_s coord_new_p data_new_p

%%%%%%%%%%%%%%% find the data points in each mesh element

nx= 40  ;  %element number along x-axis

ny= 10  ;  %element number along y-axis

nz= 35  ;  %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(y1)+(j-1)*step_y))&(y1=(min(z1)+(k-1)*step_z))&(z1

element{i,j,k}=search;     % store the index of search result into a cell array

element_center=[(min(x1)+(i-1/2)*step_x),(min(y1)+(j-1/2)*step_y),(min(z1)+(k-1)*step_z)];

record_search1{i,j,k}=search;

kk=1;

while (length(search)<4)

kk=kk+1;

search=find((abs(x1-element_center(1))

if (kk==ceil((extract)^(1/3)))&(length(search)==0)   % for empty element far from the near pts, set the up limit

b=[0,0,0,0];

break

end

end

record_search2{i,j,k}=search;  % record for test

while (length(search)>20)

pick=(1:floor(length(search)/2)).*2; % shrink # of search to 1/2

search=search(pick);

end

clear pick element_center

if (length(search)>=4) %  & length(search)<10)

if length(find(d1(search)==0))==length(search)

b=[0;0;0;0];

else

x_1=x1(search); y_1=y1(search); z_1=z1(search);

n=[ones(length(search),1),x_1,y_1,z_1];

[b,bint,r,rint,s]=regress(d1(search)',n);

clear n x_1 y_1 z_1

end

end

if length(find(b==0))==4

Q(i,j,k)=0;

else

fun=strcat(num2str(b(1)),'+',num2str(b(2)),'*x+',num2str(b(3)),'*y+',num2str(b(4)),'*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

int=sum(reshape(Q,1,[]))

time=toc