24小时热门版块排行榜    

查看: 1930  |  回复: 7

xhan

银虫 (小有名气)

[求助] matlab循环计算转为矩阵 30金币求助

程序是用来做菲涅尔积分的,由于程序运算时间太长,因此想把最耗时的两重for循环运算部分改为矩阵运算形式,但没找到合适的方法,循环中运行的语句本身就是一个矩阵运算,想请教各位怎么改成完全矩阵运算的程序。
X1=-22e-3;
X2=22e-3;
Y1=-22e-3;
Y2=22e-3;
Num0=128;
[x,y]=meshgrid(X1X2-X1)/(Num0-1):X2,Y1Y2-Y1)/(Num0-1):Y2);
for i=1:Num0
    x1=x(1,i);
    for j=1:Num0
        y1=y(j,1);
        Fij=Fun2.*exp(sqrt((x1-x).^2+(y1-y).^2)); %Fun2是一个矩阵函数
        Fun3(i,j)=sum(sum(Fij));
    end
end
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

meatball1982

铜虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
csgt0: 金币+3, 辛苦啦! 2012-11-23 09:45:47
csgt0: 金币+30, lz奖励金币 2012-12-17 09:15:17
你好。关于两个循环,实在是想不出什么好的方法,主要是求各个点到其它点的距离这一项,试了总是不行。
但如果你确定是用这一组X1,X2,Y1,Y2,就可以偷下懒。
我的程序暂时以Num0=8为例进行说明。只要X关于点对称,Y关于点对称,就好办。
你的F3是一个大的对称的矩阵。我们只要生成这个矩阵的1/8.通过转置,翻转等操作(matlab对矩阵的运算还是挺快的。),可以生成Fun3。

当Num0=8的时候,两层循环比较好。
Elapsed time is 0.000551 seconds.
Elapsed time is 0.012936 seconds.

当Num0=32的时候,两种方法差不多。
Elapsed time is 0.063168 seconds.
Elapsed time is 0.021089 seconds.

当Num0=128的时候,偷懒的方法要好些,可能是Num0越大,效果越好,但应该是有上限的,这要看你的内存值(用空间换时间就是)。
Elapsed time is 8.402686 seconds.
Elapsed time is 0.986670 seconds.


clear all
clc

% X1=-22e-3;
% X2=22e-3;
% Y1=-22e-3;
% Y2=22e-3;

% different X,Y
X1=-20e-3;
X2=20e-3;
Y1=-21e-3;
Y2=21e-3;

Num0=128;% u can change it for 8,32 64 and 128
tm_lu=zeros(Num0/2,Num0/2);
tm_lu=zeros(Num0/2,Num0);
tmp  =zeros(Num0,Num0);

Fun2=1;
x_l=X1X2-X1)/(Num0-1):X2;
y_l=Y1Y2-Y1)/(Num0-1):Y2;
[x,y]=meshgrid(x_l,y_l);
%% your way
tic
for i=1:Num0
    x1=x_l(i);
    for j=1:Num0
        y1=y_l(j);
        Fij=Fun2.*exp(sqrt((x1-x).^2+(y1-y).^2)); %Fun2是一个矩阵函数
        F3(i,j)=sum(Fij();
    end
end
toc

%%
tic
for i=1:Num0/2
    x1=x_l(i);
    for j=1:i
        y1=y_l(j);
        Fij=Fun2.*exp(sqrt((x1-x).^2+(y1-y).^2)); %Fun2是一个矩阵函数
        tm(i,j)=sum(Fij();
    end
end
tm_lu=tm+tm'-tril(tm');
tm_u=[tm_lu,flipud(tm_lu)'];
F4=[tm_u;flipud(tm_u)];
toc

% modified on : 22-Nov-2012 19:11:35
% typed by : m
唉。还是学吧。
2楼2012-11-22 19:23:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xhan

银虫 (小有名气)

引用回帖:
2楼: Originally posted by meatball1982 at 2012-11-22 19:23:09
你好。关于两个循环,实在是想不出什么好的方法,主要是求各个点到其它点的距离这一项,试了总是不行。
但如果你确定是用这一组X1,X2,Y1,Y2,就可以偷下懒。
我的程序暂时以Num0=8为例进行说明。只要X关 ...

非常感谢meatball,这样的处理确实能提高运行效率。不过对于我要做的程序还是不够,一个原因是里面用的积分最多的有二十多个,二是Num0的取值要到8000或者16000才能得到稳定解,这时候运行非常的慢,几倍的提高还是没办法解决问题。另外有个问题就是Fun2矩阵有时候不是对称的,因此得到的Fun3也不是对称的,不能用这个方法来算。不过还是很感谢你。要是后面没人能解决的话,这金币就给你了。
3楼2012-11-24 23:02:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

meatball1982

铜虫 (小有名气)

引用回帖:
3楼: Originally posted by xhan at 2012-11-24 23:02:47
非常感谢meatball,这样的处理确实能提高运行效率。不过对于我要做的程序还是不够,一个原因是里面用的积分最多的有二十多个,二是Num0的取值要到8000或者16000才能得到稳定解,这时候运行非常的慢,几倍的提高还是 ...

我再想想有没有办法,如果没能解决,金币还是算了吧。
唉。还是学吧。
4楼2012-11-25 21:58:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

himyou

银虫 (小有名气)


csgt0: 金币+1, 欢迎交流 2012-11-26 13:36:29
你用并行for试一下

parfor

» 本帖已获得的红花(最新10朵)

5楼2012-11-26 10:33:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

csgt0

荣誉版主 (著名写手)

彩色挂图

【答案】应助回帖


dbb627: 金币+1, 谢谢应助! 2012-11-27 14:36:27
菲涅尔积分用
c1=mfun('FresnelC',x)
s1=mfun('FresnelS',y)
计算,详细的你可以看看网上的一个算法,应该会快些的
http://wenku.baidu.com/view/76f633c29ec3d5bbfd0a74fd.html
showmethemoney
6楼2012-11-26 14:54:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xhan

银虫 (小有名气)

送鲜花一朵
引用回帖:
5楼: Originally posted by himyou at 2012-11-26 10:33:53
你用并行for试一下

parfor

谢谢himyou。第一次见到这个并行处理,学习了。我把程序用parfor改进了一下,虽然不能最终解决问题,但相比还是要好一些。
7楼2012-11-26 19:09:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xhan

银虫 (小有名气)

引用回帖:
6楼: Originally posted by csgt0 at 2012-11-26 14:54:08
菲涅尔积分用
c1=mfun('FresnelC',x)
s1=mfun('FresnelS',y)
计算,详细的你可以看看网上的一个算法,应该会快些的
http://wenku.baidu.com/view/76f633c29ec3d5bbfd0a74fd.html

谢谢csgt0.我仔细看了这个算法。但感觉其中的推导有些问题。第七式到第八式的推算部分直接把E0跟其它部分分别进行积分,即使不经过衍射孔这样积分也是不对的。
8楼2012-11-28 10:53:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 xhan 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[论文投稿] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 QTy3jDtz1uLt 2026-08-21 9/450 2026-08-23 03:04 by OEbVnUOu01ol
[博后之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +3 DpSrDtM079iu 2026-08-22 4/200 2026-08-23 01:04 by OEbVnUOu01ol
[基金申请] 93BebMhtakh前后11位开头都是大写 +7 且听虎啸 2026-08-17 8/400 2026-08-22 21:55 by 医学老男孩
[教师之家] 跳槽后在研项目怎么办? +3 简单化xn 2026-08-22 7/350 2026-08-22 19:46 by 简单化xn
[基金申请] 什么时候开奖? +7 CrisMessi 2026-08-18 8/400 2026-08-22 18:29 by 淀粉搬运工
[教师之家] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 2JOx3r2CYEgw 2026-08-21 4/200 2026-08-22 16:52 by sunzitan
[基金申请] 今天基金会出结果吗?20260819 +16 kkkl_v 2026-08-19 17/850 2026-08-22 16:12 by 阿布Abu
[基金申请] 放榜前的不淡定 40+4 snowwithsea 2026-08-19 14/700 2026-08-21 23:51 by cratir
[基金申请] 建议基金发布提前给出明确的时间点 +10 kulium 2026-08-21 13/650 2026-08-21 21:58 by alongwaytogo
[基金申请] 今日不放榜?网传国自然预计 8 月 27 日可查结果 +16 医学老男孩 2026-08-20 20/1000 2026-08-21 21:13 by Ldrop2023
[基金申请] 科研孤儿太难了 +17 我4大白菜 2026-08-20 18/900 2026-08-21 20:57 by zhangev
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 我面上完蛋了 +7 且听虎啸 2026-08-20 8/400 2026-08-21 12:31 by 酷酷墨镜
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
[基金申请] 今天系统多次维护,明天很可能放榜! +10 zju2000 2026-08-16 11/550 2026-08-20 20:02 by cl479861084
[基金申请] 重要消息,中午系统在维护 +11 yuleib84 2026-08-18 12/600 2026-08-20 11:09 by xskun
[基金申请] 今天放榜吗? +14 布布和一二 2026-08-19 15/750 2026-08-19 18:07 by gltch
[基金申请] 朋友圈看到的 +6 wangzilk 2026-08-18 8/400 2026-08-19 10:55 by Haru815
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
信息提示
请填处理意见