应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » matlab循环计算转为矩阵 30金币求助
81/1返回列表
查看: 1562  |  回复: 7
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2012-11-22 13:33:23
已阅   回复此楼   关注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朵）

xhan
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 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭