| 查看: 504 | 回复: 3 | |||
| 当前主题已经存档。 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[交流]
【求助】 Matlab 数值积分
|
|||
|
求教数值积分 对于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 不胜感激. [search]数值积分[/search] [ Last edited by sunxiao on 2009-3-9 at 08:48 ] |
回复此楼» 猜你喜欢
有没有人能给点建议
已经有5人回复
-
假如你的研究生提出不合理要求
已经有12人回复
-
实验室接单子
已经有7人回复
-
全日制（定向）博士
已经有5人回复
-
萌生出自己或许不适合搞科研的想法，现在跑or等等看？
已经有4人回复
-
Materials Today Chemistry审稿周期
已经有4人回复
-
参与限项
已经有3人回复
-
对氯苯硼酸纯化
已经有3人回复
-
所感
已经有4人回复
-
要不要辞职读博？
已经有7人回复
|
目前我采用如下方法计算三维积分, 但是出现积分结果随着剖分单元数量增加先下降收敛, 然后反弹上升的曲线. 大家能否帮我指出思路中的错误, 先谢谢了. 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 |
3楼2008-11-14 13:51:12
2楼2008-11-08 12:02:34
4楼2008-11-18 11:38:04