arrayfun计算积分的问题

liuyalong008

这是在中文论坛看到的一个问题，原本觉得使用arrayfun解决应该会很好，但是还需要转来转去的挺麻烦，可好似发现计算要慢，不知论坛各位高人有没有更好的办法

%loop

1. clear
   
2. clc
   
3. tic;
   
4. i=0;
   
5. for a=-10:10
   
6. i=i+1;
   
7. j=0;
   
8. f1=inline('1./(1+x).*exp(-1*(a-x).^2./(1+x).^2)');
   
9. f11=vectorize(subs(f1));
   
10. z1=quad(f11,0,20);
    
11. for b=-10:10
    
12. j=j+1;
    
13. f2=inline('1./(1+x).*exp(-1*(b-x).^2./(1+x).^2)');
    
14. f22=vectorize(subs(f2));
    
15. z2=quad(f22,0,20);
    
16. z(i,j)=z1*z2;
    
17. end
    
18. end
    
19. a=-10:10;
    
20. b=-10:10;
    
21. [x y]=meshgrid(a,b);
    
22. surf(x,y,z)
    
23. 
    
24. toc;

复制代码

Elapsed time is 6.756710 seconds.

%arrayfun

1. clear
   
2. tic;
   
3. 
   
4. syms x x0 y y0
   
5. M=1/(x0+1)*exp(-(x-x0)^2/(1+x0)^2);
   
6. N=1/(y0+1)*exp(-(y-y0)^2/(1+y0)^2);
   
7. [X,Y]=meshgrid(-10:10,-10:10);
   
8. fstr1=char(vectorize(inline(M)));
   
9. fstr2=char(vectorize(inline(N)));
   
10. fstr=sprintf('quadl(inline(%s),0,20)*quadl(inline(%s),0,20)',fstr1,fstr2);
    
11. fun=inline(fstr,'x0','x','y0','y');
    
12. z=arrayfun(@(x,y)fun(x0,x,y0,y),X,Y);
    
13. 
    
14. surf(z)
    
15. toc;

复制代码

Elapsed time is 24.318110 seconds.

qibbxxt

这是我看到题目第一个写的程序

1. % 循环计算
   
2. clear;clc;close all
   
3. tic;
   
4. f=@(x)quad(@(x0)1./(1+x0).*exp(-(x-x0).^2./(1+x0).^2),0,20);
   
5. w=@(x,y)f(x).*f(y);
   
6. ezmesh(w,21);
   
7. toc;
   

复制代码

运行结果为：

1. Warning: Function failed to evaluate on array inputs; vectorizing the function may
   
2. speed up its evaluation and avoid the need to loop over array elements.
   
3. > In specgraph\private\ezplotfeval at 57
   
4.   In ezgraph3>ezeval at 611
   
5.   In ezgraph3>surfplot at 538
   
6.   In ezgraph3 at 49
   
7.   In ezmesh at 66
   
8.   In ex03 at 6
   
9. Warning: Function failed to evaluate on array inputs; vectorizing the function may
   
10. speed up its evaluation and avoid the need to loop over array elements.
    
11. > In specgraph\private\ezplotfeval at 57
    
12.   In ezgraph3>ezeval at 611
    
13.   In ezgraph3>gridfilt at 259
    
14.   In ezgraph3>surfplot at 543
    
15.   In ezgraph3 at 49
    
16.   In ezmesh at 66
    
17.   In ex03 at 6
    
18. Elapsed time is 1.892823 seconds.
    

复制代码

于是就按照警告，向量化

1. clear;clc;close all
   
2. tic;
   
3. f=@(x)arrayfun(@(i)quad(@(x0)1./(1+x0).*exp(-(x(i)-x0).^2./(1+x0).^2),0,20),1:numel(x));
   
4. w=@(x,y)f(x).*f(y);
   
5. ezmesh(w,21)
   
6. toc;
   

复制代码

运行结果：

1. Elapsed time is 2.310334 seconds.
   

复制代码

为了不同做法的效率，写了以下代码：

1. clear;clc;clear all
   
2. tic;
   
3. [a,b]=meshgrid(-10:10);
   
4. c=zeros(size(a));
   
5. f=@(x)quad(@(x0)1./(1+x0).*exp(-(x-x0).^2./(1+x0).^2),0,20);
   
6. w=@(x,y)f(x).*f(y);
   
7. for i=1:numel(a)
   
8.     c(i)=w(a(i),b(i));
   
9. end
   
10. surf(a,b,c);
    
11. toc;
    

复制代码

这个程序的运行结果为：

1. Elapsed time is 0.984482 seconds.
   

复制代码

1. clear;clc;close all
   
2. tic;
   
3. [a,b]=meshgrid(-10:10);
   
4. f=@(x)arrayfun(@(i)quad(@(x0)1./(1+x0).*exp(-(x(i)-x0).^2./(1+x0).^2),0,20),1:numel(x));
   
5. w=@(x,y)f(x).*f(y);
   
6. surf(a,b,reshape(w(a,b),size(a)));
   
7. toc;
   

复制代码

运行结果：

1. Elapsed time is 0.863089 seconds.
   

复制代码

1. 
   
2. clear;clc;close all
   
3. tic;
   
4. [a,b]=meshgrid(-10:10);
   
5. c=zeros(size(a));
   
6. for i=1:numel(a)
   
7.     c(i)=quad(@(x0)1./(1+x0).*exp(-(a(i)-x0).^2./(1+x0).^2),0,20)....
   
8.         *quad(@(x0)1./(1+x0).*exp(-(b(i)-x0).^2./(1+x0).^2),0,20);
   
9. end
   
10. surf(a,b,c);
    
11. toc;
    
12. 
    

复制代码

运行结果：

1. Elapsed time is 1.049432 seconds.
   

复制代码

然后看了你的程序，先说说关于你的程序，我的个人愚见
就是roc在上个帖子中所说，用syms导致你的速度降低，值得肯定的是你的char,sprintf,以及匿名函数来构造使用arrayfun还是值得大家借鉴的(inline就不推荐用了)
以上的例子大家似乎可以看出效率方面的差异，当然不大，可以调大网格数目，
测试得到更加有力的说法
arrayfun函数的效率似乎并不那么高，不像accumarray那样，
用arrayfun函数现在的目的更多只是简洁
以上只是个人看法，希望大家一起讨论

liuyalong008

qibbxxt 版主确实厉害，能列举出这么多种方案，学习了
另外，第一二种方案，我的明显要慢，要9秒多，其他的都差不多0.6秒左右
应该是matlab版本的问题吧，我的是2010b

qibbxxt

3# liuyalong008
我后来忘了改程序了，第1，2慢的原因是网格太密，我改一下，你再试一试

liuyalong008

3# liuyalong008  
我后来忘了改程序了，第1，2慢的原因是网格太密，我改一下，你再试一试
qibbxxt 发表于 2010-12-27 14:40

非常感谢，其中第三个是时间最快的

PS：绝对的高手

wanglu

用楼主的公式测试了一下Forcal的arrayfun函数的效率，发现比while循环慢些，可能是矢量化代码要分配内存的缘故吧？

1. 
   
2. !using["math","sys"];
   
3. fx0(x0::x)=1/(1+x0)*exp{-[(x-x0)^2/(1+x0)^2]};
   
4. fy0(y0::y)=1/(1+y0)*exp{-[(y-y0)^2/(1+y0)^2]};
   
5. w(xx,yy::x,y)= x=xx,y=yy, IMSL::QDAGS[HFor("fx0"),0,20,0,1e-6,0]+IMSL::QDAGS[HFor("fy0"),0,20,0,1e-6,0];
   
6. ArrayFunTime(:x,y,a,t0) =
   
7. {
   
8.   t0=clock[],
   
9.   oo{ndgrid[linspace(minus(10),10,21),linspace(minus(10),10,21),&x,&y], a=arrayfun[HFor("w"),x,y].Sum[]},
   
10.   printff{"\r\n结果＝{1,r}  arrayfun时间＝{2,r}毫秒",a,sys::clock()-t0}
    
11. };
    
12. ForTime(:i,j,a,t0)=
    
13. {
    
14.   t0=sys::clock(), a=0,
    
15.   i=-10,(i<=10).while{
    
16.     j=-10,(j<=10).while{
    
17.       a=a+w[i,j],
    
18.       j++
    
19.     },
    
20.     i++
    
21.   },
    
22.   printff{"\r\n结果＝{1,r}  while时间＝{2,r}毫秒",a,sys::clock()-t0}
    
23. };
    

复制代码

结果：
结果＝885.98422777071755  arrayfun时间＝125.毫秒
结果＝885.98422777071755  while时间＝78.毫秒