求教 Matlab 三维数值积分

[cross]

求教数值积分

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

各位如果有时间, 也可以帮忙回答到这个邮箱里面 [email protected]
不胜感激.

[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

[cross]

恩, 我还是没有想出来错在哪里, 把坐标上的函数值设定为1, 这样用结果与积分空间的体积作比较, 网格增多, 结果保持等于体积, 说明我的算法没有额外加入点使结果偏大. 但是仍然不理解为何收敛会出现回升.

[laosam280]

请把你的完整的问题发给我，我用并行算法给你试试。50万个数据点的操作对于普通的电脑是根本承受不了的。

[laosam280]

我的邮箱[email protected]