24小时热门版块排行榜    

查看: 1921  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 今天基金会出结果吗?20260819 +14 kkkl_v 2026-08-19 15/750 2026-08-20 08:21 by sunzitan
[基金申请] 欢迎发来filecode的Mz6后的代码验证其规律 +55 医学老男孩 2026-08-13 123/6150 2026-08-19 23:44 by stan_et
[基金申请] 时间戳变了,能看出什么问题? +12 基诺咪客 2026-08-17 16/800 2026-08-19 22:04 by kk 的小木虫
[基金申请] filecode,4个jtjc了 +11 ziyangfang 2026-08-19 13/650 2026-08-19 21:28 by aasahr
[教师之家] 为什么余额宝的年化利率越来越低?主要原因有哪些? +5 瞬息宇宙 2026-08-15 5/250 2026-08-19 20:23 by super2002521
[基金申请] 今天放榜没戏了吧 +4 yuleib84 2026-08-19 5/250 2026-08-19 19:52 by hhs666
[基金申请] 今天放榜吗? +14 布布和一二 2026-08-19 15/750 2026-08-19 18:07 by gltch
[基金申请] 2027广东省杰青 +4 奶牛小黑 2026-08-15 10/500 2026-08-19 11:02 by wanfengnew
[基金申请] 快农历七夕节了,轻松一下,男人悄悄话,女施主请不要进来。 +6 Tide man 2026-08-14 7/350 2026-08-19 09:56 by ZJTJZ
[基金申请] 什么时候开奖? +6 CrisMessi 2026-08-18 6/300 2026-08-19 08:06 by Equinoxhua
[基金申请] 重要消息,中午系统在维护 +10 yuleib84 2026-08-18 11/550 2026-08-19 07:30 by 羽毛枫f
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[论文投稿] 投稿咨询 +4 wwm09 2026-08-17 6/300 2026-08-18 15:36 by wwm09
[基金申请] 时间戳又变了8-15 +14 archvillain 2026-08-15 26/1300 2026-08-18 13:37 by phantomgost
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
[基金申请] 哪位老哥知道今年的国自然具体哪一天放榜? +13 Ldrop2023 2026-08-13 16/800 2026-08-18 12:25 by 淀粉搬运工
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
[基金申请] 感觉是下周放榜了 +6 angus9576 2026-08-17 11/550 2026-08-17 23:57 by angus9576
[基金申请] 今天系统多次维护,明天很可能放榜! +8 zju2000 2026-08-16 9/450 2026-08-17 12:20 by lmz0216
[基金申请] 各位道友,我要去昆明玩几天,回来见。 +7 Tide man 2026-08-14 8/400 2026-08-15 01:11 by arzu_hma
信息提示
请填处理意见