数学建模社区-数学中国
标题:
极限测试之Matlab与Forcal代码矢量化
[打印本页]
作者:
forcal
时间:
2011-8-2 15:02
标题:
极限测试之Matlab与Forcal代码矢量化
代码矢量化是matlab的特色,但这点似乎不难实现。代码矢量化的优势并不明显,通过一个例子说明。
//用C++代码描述为:
) \* O: o8 H2 H. e! ^% ]
s=0.0;
( Q% E2 r$ `4 {# N& H
for(x=0.0;x<=1.0;x=x+0.0011)
3 \3 V! K. d6 `5 |' D3 o
{
h9 L; Z& S: q2 R- W
for(y=1.0;y<=2.0;y=y+0.0009)
4 W# s1 ~7 K0 C& G
{
( y8 D h' a. S6 V% b
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
& j; r( U, }6 o, d
}
# m0 P8 v! w5 Z7 d5 d* w2 b
}
复制代码
Matlab代码:
tic
, K$ E! H2 T( b& f
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
1 A" w2 @) T1 q* o. f8 t7 |
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
T' S: o$ T" I/ }$ l+ O
toc
9 M h% _$ s5 H
9 @; e6 Y& T; x V# m
s =
! z; F/ w4 H, [* h) a
8 ~7 c) J; Z5 D( r3 D
1.0086e+006
) h0 a x* _9 U, n# ?
- ^2 |6 m% G, c9 a
Elapsed time is 0.561108 seconds.
复制代码
Forcal代码:
!using["math","sys"];
; ?3 M+ s$ T$ ~5 o
mvar:
0 H l+ I+ A) m* m. l3 ~2 b
t=clock(),
5 z* @. Z5 s. a$ b, h
oo{
u9 m7 w. O i* }6 P$ S- {
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
0 J* ?0 c$ ]# q- ]2 t
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]
( _+ S+ D/ n, ^/ d
};
2 |8 U G: v+ Y) W7 ~* x! q
[clock()-t]/1000;
复制代码
结果:
9 e8 q4 |& ~1 {/ G. {; Q
1008606.64947441
& P. N3 f5 ?8 f: _$ n, I8 x
0.641
# J% x3 j1 x" S- {+ k! H
) }( U9 H6 l" c* B6 j1 F8 R7 X) r
Forcal比Matlab稍慢些。
$ ]0 Q4 B Q# n$ e+ i4 r
! {; J/ z2 h3 i: p [
----------
, b8 T2 x8 w) u5 y' t2 w; d
, e- [) }6 S) L0 K
再看循环效率。
6 H* `% _6 y, g, T8 u, F& U5 `
1 _& }; X' Q x; x8 K4 X; C) {* a
Matlab代码:
tic
1 i' h8 I$ M" M7 V0 r$ Q; E
s=0;
% X8 h) g% a: u& \* T1 ~
x = 0;
: @) m6 Z- K: N5 l
y = 1;
: I9 q8 R( v! g. D" U9 a C3 ~
while x<1
$ b8 Y/ ?1 y5 m" O' p
while y<2;
, L' R0 ?# M! p5 A$ H. }3 I
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
; _9 i+ I& o! y# F z: v3 O4 A
y = y+0.0009;
8 q9 Z7 F- N* q% U, i& ^1 y
end
5 H. y: `0 j1 z, A
x = x+0.0011;
, X! X- ?$ b9 O& Y
y = 1;
* f$ m/ k4 k' f6 b9 m( B3 _
end
a& j1 k/ b6 v
s
( Z" g4 l# g5 j" D8 D$ o; ~
toc
1 @! {8 A( k: s9 W+ j7 m
- T- B3 s0 d- Z% F
s =
8 a( H: {' A$ c3 K5 `4 d3 t ?
* L% e1 j0 X4 h6 A
1.0086e+006
1 j2 @9 F' Q1 u% A
% ]) `" }7 i& _, i0 E1 L: z) |
Elapsed time is 0.933513 seconds.
复制代码
Forcal代码:
mvar:
' E: w, h% @, Y
t=sys::clock();
- V& Q3 I5 O% T) t' h/ u# Y1 n
s=0,x=0,
4 L5 c3 R( L- t: k9 j
while{x<=1, //while循环算法;
6 U) L+ G" J; R# n
y=1,
1 J1 j- D" V3 m: }6 ?
while{y<=2,
& f; n) O. S% b7 k, l+ @
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
: Y/ {" d) ~1 d. H2 y# _
y=y+0.0009
7 B8 Q E0 U5 j5 r1 I
},
% |6 _0 c" d! S& r9 k
x=x+0.0011
, ?5 \1 C+ G& S# P: k2 G
},
) B. Y: `3 d/ N- v; {
s;
0 i& m4 n5 Y6 B' R. l |; h4 l
[sys::clock()-t]/1000;
复制代码
结果:
4 m$ A4 e' f: n1 d0 l$ [, M3 ^
1008606.64947441
7 o) v3 E+ v. y: H9 w4 W
0.734 //时间,秒
& y& ~5 d4 C2 i/ e' f% f
H: @0 W" r- s5 Y0 Q" b- J
我很奇怪,在这个例子中,matlab的JIT加速器为什么没有起作用?
+ I1 {* N1 I+ o1 o* ?# S
/ w' x* D: A; ^ c, f4 r4 d
-------
) F0 t' e$ N D; ~' D. y
. k) L$ z9 H4 J- ^; w
Forcal中还有一个函数sum专门进行这种计算:
mvar:
L& E3 y3 G) y. `0 p( q
t=sys::clock();
* u, M2 S. j: l: j
f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
3 h# ]9 _# a+ V3 \* b9 D* F
sum["f",0,1,0.0011 : 1,2,0.0009];
) _' _3 Z' w& m/ R- `
[sys::clock()-t]/1000;
复制代码
结果:
" P9 x3 {+ O$ {: V0 T, _
1008606.64947441
I8 G" p1 [/ L; Y3 E
0.719 //时间,秒
作者:
forcal
时间:
2011-8-2 15:24
说到代码矢量化,就不能不提到大名鼎鼎的arrayfun函数,现在再来看看它的表现。
& z4 V) }1 ?1 m( v5 o; P# s( _# @, S
8 \0 u e& b3 @- F) \( v
matlab代码:
f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
% Q! l) |) }/ N& @. H) Z
tic
# N! ~- i: \$ o, g+ t/ ~
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
) M" T$ o8 c, E
sum(sum(arrayfun(f,x,y)))
; D5 }0 F0 O' _" u- {) M' z, [$ H
toc
4 z, H% O* M, Q0 Y" z
/ ]) k$ W$ y6 C8 f) }9 j% x" F
ans =
1 }2 t9 W* ]/ d- ]- U1 ]
+ N/ X2 \8 j4 G/ D/ t" j5 ?9 B1 X
1.0086e+006
+ F) F$ Z% q# d, i% r
/ ~$ e# M2 o# z& s- c
Elapsed time is 7.339923 seconds.
复制代码
Forcal代码:
!using["math","sys"];
) r/ `4 D. w+ d% |, g
mvar:
1 ~ ]1 d+ w/ c) f0 L, ]% t
f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
1 ~7 K( E. G) I
t=clock(),
: m: f/ f' T/ _3 _# q# m
oo{
+ w! c/ p2 L. c6 t4 X
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
2 e, v% s; x6 H5 n/ Y+ c- l
arrayfun[HFor("f"),x,y].Sum[]
$ }, }9 T1 \( Y% D5 {6 S
};
4 [5 \1 Z4 |6 w+ v$ S5 F
[clock()-t]/1000;
复制代码
结果:
7 a; ^, d/ l7 o0 i
1008606.64947441
h- I0 z# W Q5 \
0.735 秒
; Q% V/ {' Q! _3 u
4 A$ Q x; R, u# B
可以看出,Forcal的arrayfun函数远快于matlab的arrayfun函数。
9 p: J( O1 O' j- {4 g N" k
# [! J$ ~* S( o, b
--------
* i S2 z5 n3 O
: U' K% M4 Z, O
从这里似乎可以看出,matlab若借助于arrayfun函数计算三重及以上积分,其效率将远远落后于Forcal。当然,Forcal不使用arrayfun函数计算三重及以上积分。
作者:
我就在你背后
时间:
2011-8-2 16:10
学习中。。。。。。。。。
作者:
海水
时间:
2011-8-2 18:00
看看看。。。。。。。。。。。。。
作者:
alair005
时间:
2012-2-7 12:38
恩,参考一下。。
1449529676012957
作者:
sxjm567
时间:
2012-11-20 08:56
一起交流!楼主给咱们提供机会了
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5