数学建模社区-数学中国
标题:
极限测试之Matlab与Forcal代码矢量化
[打印本页]
作者:
forcal
时间:
2011-8-2 15:02
标题:
极限测试之Matlab与Forcal代码矢量化
代码矢量化是matlab的特色,但这点似乎不难实现。代码矢量化的优势并不明显,通过一个例子说明。
//用C++代码描述为:
}. |, e4 O! C4 D7 y6 ]) t: T
s=0.0;
' X, ~- y4 l8 L
for(x=0.0;x<=1.0;x=x+0.0011)
" L# Q; L2 S' T! N' M
{
. T: J# {# q" c; b- s `& f2 s
for(y=1.0;y<=2.0;y=y+0.0009)
, t. w) C+ Y- r# p% D) E4 ~: L8 U1 r
{
- ^% \8 T0 ~* A& e2 T" o
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
$ ?7 I3 d9 e3 C8 k9 z+ C2 v
}
! V7 I0 { M) }1 A2 I
}
复制代码
Matlab代码:
tic
) Z3 \! f* p; t/ i5 h3 L" L
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
3 N, e+ p# I& q y! v) p+ Q* g5 F9 e
s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
7 m) a3 D) h1 p
toc
) i) l" t' P8 Y k! S
, p: d5 j5 G/ X$ v
s =
3 T4 }3 E% p6 [+ J0 b. k
6 f% q$ p4 j6 \' Q$ Q
1.0086e+006
' N! ?3 w. }4 i1 g
1 `: H% h% n: T
Elapsed time is 0.561108 seconds.
复制代码
Forcal代码:
!using["math","sys"];
0 d$ G3 h* n5 v8 Q' G
mvar:
2 j" U$ A2 q$ ^) P( f2 ~/ R' [7 F
t=clock(),
. E# @- G9 b4 [2 U
oo{
9 H) m# p! ^' c, v: T0 U, V) y' B
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
( x. ^+ `7 O, D& Q
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]
* W% E$ i9 _( p, M
};
: s! R$ y: p; g
[clock()-t]/1000;
复制代码
结果:
" K! N* H) R/ D9 [ R- S1 n/ c
1008606.64947441
$ F) @1 j$ q5 @9 a) Y
0.641
3 f0 J; `2 s6 k
5 g* Z2 b4 k& W& a% ^. E
Forcal比Matlab稍慢些。
% Q. P4 {% M% P' }+ k8 @
& W1 e' e* x5 E" n# T9 w
----------
" D# Z/ c; l5 k9 [8 V( `
. U) h& m! o: M- S5 F
再看循环效率。
- }# u- F4 \% i- l [' R" q3 s1 U% l
8 z4 e0 s+ j1 e! @
Matlab代码:
tic
. I* l$ t2 _7 p; o2 a4 X8 {/ \
s=0;
" n1 `. M1 P+ a: V
x = 0;
0 I3 ~% B; K; H D! h' l' T
y = 1;
( w* `. O, I6 ]0 t/ d6 N$ X: P
while x<1
2 h; W* ~* L5 y4 r" ~
while y<2;
, E/ t. i- o& R5 l) w
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
j# W( A6 {) ?) s* t8 R
y = y+0.0009;
9 G; ~3 j6 ?4 i; [. u" o: F
end
\$ J6 M6 p- `9 ~0 ?3 e
x = x+0.0011;
" T7 @) z O* A# e
y = 1;
( v S. R) ~& E
end
$ J4 g+ X6 ?- W3 U
s
~0 a: ^/ ^5 Z7 S/ I$ ~" g0 d' y; k
toc
8 O1 f4 u1 f" o: ^, a: }
6 N/ A6 U5 X+ M" w
s =
5 y9 A7 M X+ d+ W' A/ G
% e1 A% R4 b! {7 I
1.0086e+006
x% v/ |% F2 N2 G1 i* m, [
5 G9 N& F. D4 ]5 y7 j6 }% [
Elapsed time is 0.933513 seconds.
复制代码
Forcal代码:
mvar:
9 I2 X. c/ o. i2 T7 o, M' ?. W
t=sys::clock();
9 C, _: ?' ^3 h- t
s=0,x=0,
+ P4 L# c4 J! ~7 K; e( T
while{x<=1, //while循环算法;
3 Z0 E% O7 w3 _$ N6 R% k+ r% S
y=1,
5 H: Y+ x6 F& e) |" S4 ]' }- M; Z
while{y<=2,
+ ^* ^( v" | f/ i
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
; }4 d3 o# W: c; h$ `
y=y+0.0009
5 a9 H/ Y7 ]( l1 L/ C% l
},
3 z0 l7 t' q& q& K4 E: b* @0 X
x=x+0.0011
8 J5 G3 m& m4 r3 N5 L
},
9 o, J i" K5 _. M# z
s;
a' g% Y' [' F* M
[sys::clock()-t]/1000;
复制代码
结果:
. q: X4 n, s' ~6 c6 y4 l
1008606.64947441
7 J1 g9 e& Y9 L
0.734 //时间,秒
- G1 |; a: T) |8 {: j0 q: t
- G: C! [7 @/ L, F
我很奇怪,在这个例子中,matlab的JIT加速器为什么没有起作用?
$ N$ n+ W0 n- _$ V7 t6 s1 [* t( F
/ A4 z1 E5 t9 F
-------
* e( M+ r& U* ~
* V& e, A5 W: O2 q
Forcal中还有一个函数sum专门进行这种计算:
mvar:
' B! w" h o3 q5 i: e9 [# K
t=sys::clock();
1 e# c: l8 P2 N% `
f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
7 A1 u2 C( B- ]% W
sum["f",0,1,0.0011 : 1,2,0.0009];
- R/ M8 v) j) v( Q' {" s O
[sys::clock()-t]/1000;
复制代码
结果:
& C) X, f; T9 U3 c n
1008606.64947441
1 G7 p6 E# U. _8 l0 V. i+ v
0.719 //时间,秒
作者:
forcal
时间:
2011-8-2 15:24
说到代码矢量化,就不能不提到大名鼎鼎的arrayfun函数,现在再来看看它的表现。
6 Y$ d- d( T8 x$ y
$ A& ~4 g2 a; ]& y
matlab代码:
f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
: L: {; P) X5 \, c4 q+ t, z
tic
. _6 U/ z, Y0 t3 G. a, ^
[x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
0 M% O( N4 e3 |9 j
sum(sum(arrayfun(f,x,y)))
$ ~4 M+ ~0 a4 Z9 g( B9 r+ M; L
toc
- @7 J+ r* Q0 t: L( F' Z3 V" z* n: `
! f* @" O0 \- K3 ^" E
ans =
; e$ U7 K, c R0 I$ l; s* o% S
6 V& W4 w$ J$ m7 C! V" a
1.0086e+006
% G% Q2 M8 [8 C
6 E' W9 @* {( h0 _* U( B) j) w* @
Elapsed time is 7.339923 seconds.
复制代码
Forcal代码:
!using["math","sys"];
2 V" N A$ I, L* r L. R
mvar:
, u1 j! O+ S j4 _: w! D2 ^
f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
2 X- \3 N2 m! O) n
t=clock(),
6 s8 }+ g% s: D: V
oo{
g6 B8 Y% }+ W* S" K, j3 O% ^/ v& T
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
$ |4 x! E7 C; P! B9 m( O! G4 |
arrayfun[HFor("f"),x,y].Sum[]
' u, H% J- s2 j7 J; D! G/ C- }
};
/ n8 d/ `- j$ Y1 _3 z9 V; M
[clock()-t]/1000;
复制代码
结果:
' g5 [' z" ^7 u7 K
1008606.64947441
; a, w7 s1 y t$ R: q/ G1 M- r0 h
0.735 秒
, n( o2 `) K9 _# j
: h: _) I: p+ f, E) M; H
可以看出,Forcal的arrayfun函数远快于matlab的arrayfun函数。
7 h; }% q4 Y# G# I0 \. H8 A
$ V% l' p4 S/ q. d# r: ?
--------
! ^8 a! m: b2 D9 f! {0 d- }
" M0 r! w8 {9 j) `6 e* Q5 h
从这里似乎可以看出,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