【原创】：新年礼物---戏说arrayfun系列函数

bainhome

                                   戏说Arrayfun系列函数
摘要：本文介绍了MATLAB软件中自7.0版本之后的一个新函数系列，包括：“arrayfun”、“cellfun”以及“strfun”等，对该函数系列的基本使用方法、应用范围选择、与MATLAB传统循环流程的比较以及在具体数学或者工程问题中的应用方法做简单介绍与分析。



§1 “*fun”系列函数简介
1.1 关于“*fun”系列函数的名称
这里的““*fun”系列函数”是指实际用法和原理基本相同的几个函数如：arrayfun、cellfun和structfun的统称，它们均用于对某个或者某些函数句柄进行批量调用。当然MATLAB中还有针对稀疏矩阵函数句柄进行调用的函数“spfun”，但一般只在特定场合使用，本文介绍从略。
“*fun”函数中，从基本原理上讲，可看作一种特殊循环形式，只是循环体内依循规律改变的对象多为函数，往往不是一个或者一系列具体的“数”。不过，对于具体数字的循环也同样可以调用这些函数实现。以最简单的多点求值为例，代码与运行结果如下：

1. x=[1 2 3];                                   %给定数据
   
2. yNormal=x.^2                                 %普通求值
   
3. yArFun=arrayfun(@(i1) x(i1)^2,1:length(x))   %调用arrayfun函数求值
   
4. xCldata={pi,pi/2,pi/6,pi/4};                 %定义cell数据
   
5. yClFun=cellfun(@sin,xCldata)                 %以cellfun调用"sin"方法求值
   
6. strData=struct('a1',5,'a2',6,'a3',7);        %定义结构数组
   
7. ystrFun=[structfun(@(x1) sin(2*x1),strData)]'%以structfun函数调用"sin"求值

复制代码

执行结果：

1. yNormal =
   
2.      1     4     9
   
3. yArFun =
   
4.      1     4     9
   
5. yClFun =
   
6.     0.0000    1.0000    0.5000    0.7071
   
7. ystrFun =
   
8.    -0.5440   -0.5366    0.9906  
   

复制代码

从结果可以看出：“*fun”系列函数均可对数据实现批量求值。不过，这不是本文的重点，更远远不是这些函数真正应当发挥作用的场合，因为从函数名称中的“fun”字段，我们能得出初步判断：这些函数最拿手的本领，应该是在特殊的结构体内操纵函数句柄实现循环！

bainhome

在进行具体应用“*fun”系列函数进行实例计算之前，我们有必要先看看MATLAB的help是怎么描述这些函数的，因为从这些内容中，我们可以从轮廓上看出，编写这些函数的程序员希望函数如何被用户使用。不妨以平时使用频率最高的arrayfun函数为例：
================================================
Arrayfun函数的定义：将确定的函数方法作用于数组中的每个元素
================================================
解释：
确定的函数方法：相当于处理不同材料的机床刀具；
数组中的每个元素：相当于待加工原料
arrayfun：相当于机床

既然是“作用于每个元素”，这当然暗示着一次循环遍历，所产生的结果当然就是arrayfun"流水线"上的最终产品，因此不难在脑海中生成如图 1.1所示的图景：

                                                         
                                                 图 1.1 Arrayfun函数处理数据的方式图解

bainhome

为了更详细地说明arrayfun函数的使用方式，下面给出类似“hello world！”的示例，由于国家体育总局“注意感谢对象”的友情提示，使我们深刻认识到传统的hello world已经不符合中国国情，因此全面修改，以期配合当前形式的发展需要：

1. 
   
2. clear;clc;close all
   
3. strThanks={'城管','发改委','人民公安', '国税局','和谐社会'};
   
4. str='感谢';
   
5. arrayfun(@(x)strcat(str,strThanks{x}),1:length(strThanks),'UniformOutput',false)'  

复制代码

执行结果：

1. ans =
   
2.     '感谢城管'
   
3.     '感谢发改委'
   
4.     '感谢人民公安'
   
5.     '感谢国税局'
   
6.     '感谢和谐社会'

复制代码

从上述示例可以看出：“确定的函数方法”是“@strcat ”，待加工的原料当然就是感谢对象：strThanks字符串数组中的每个数组元素相当于等待“被”感谢的伟大对象，在arrayfun这台机床上，我们首先换上“@strcat”刀具，命令它对五个元素进行同样方式的加工---让“感谢”字符和每个字符串数组元素实现“并”组合。于是，就实现了和谐而有条不紊地感谢。于是，世界就安静了。代码中每个元素的基本功用如图 1.2所示：

                                             
                                            图 1.2 arrayfun函数对strcat的操作原理图示
但是，arrayfun显然并不是实现这一功能最合理的方法，在MATLAB中有矢量化的字符串操作函数：strcat，实现代码如下：

1. str=strcat({'感谢'},strThanks')

复制代码

怎么样，是不是更加简单，关于这个问题后面还有几个例子，所以把总结放在后面，现在看官不妨跟着笔者的思路先走一走。

bainhome

我想，当你认认真真按上述方法虔诚地“感谢”过之后，关于句柄和匿名函数的基本原理就应该已经有一个比较明确的“直觉”了。之所以说仅仅是直觉而非一个明确理论，是因为如果把句柄和匿名函数真正的定义拿出来，还要按术语分析解释其内部的数据结构、重载方式等，这实在是杀死脑细胞的优良方法，另外也小声说一句：“我们和你一样，对这些也不懂…”。
但至少，这些看似深奥晦涩的理论不能阻止我们把各种函数如“@sin”、“@plot”、“@strcat”等等，想象成为一些基本的刀具库，它们就是机床上的铣刀、车刀…(当然还有它们下属的各种型号刀具)，而arrayfun这台万能机床，则不管你是什么刀，只要这十八般武器你敢拿来放在原料(待处理数组)的适当位置上，我就敢一样一样按你给定的规则，去一样一样耍那些待处理的数组原料。同样：cellfun，structfun也和arrayfun一样，可以看成是另外两台基本原理相同的万能机床，只是处理的原料不同，一个是处理cell数组、一个处理结构数组而已，从调用的“函数刀具库”角度来看，三者相同。
至于函数句柄就更简单了，这是让arrayfun和你知道到底抽出来的是弓还是剑、是刀还是枪的标识，有人说这个句柄官方解释太抽象，让人摸不着头脑，很不喜欢！我只能举个也许并不恰当的例子：“你家菜刀能不能没有刀把”？“信用卡能不能不分名字”？必要的标识是使用何种工具的重要信号，而这个信号，在计算机中，就可以看成是索引，最终给出一个名字“句柄”，这是向计算机索要某种确定工具使用权的消费存根。所以，有了句柄，就使得我们可以明确地告诉arrayfun，我们到底要干什么，甚至利用简单工具，实现多个函数方法操作的组合。而这个组合，就叫做“匿名函数”。
这个名字听起来有点儿鬼鬼祟祟的。废话不多说，先用一个简单例子说明一下，以下列函数为例：
从MATLAB角度来讲，使用了种不同的函数方法，包括：正弦、开根号、÷和余弦，此时就不好麻烦arrayfun它一个一个方法去调用，哲人说过：“在简洁面前，神马都是浮云”，因此，祭出匿名函数法宝，把这些处理方法都组合起来，形成一个“函数刀具”的组合，再把这个组合的句柄信息扔给arrayfun，告诉它：“原料来了之后，就按这个我所给的这个信息执行吧”。基本上就是这个意思，具体代码如下：

1. f=@(x) sin(x)+sqrt(x./cos(x))                %定义匿名函数
   
2. x=pi*[1/12:1/20:1/3]                         %给出数组原料
   
3. Results=arrayfun(@(ii) f(x(ii)),1:length(x)) %arrayfun按数组长度往复加工
   

复制代码

执行结果：

1. f =
   
2.     @(x)sin(x)+sqrt(x./cos(x))
   
3. x =
   
4.     0.2618    0.4189    0.5760    0.7330    0.8901    1.0472
   
5. Results =
   
6.     0.7794    1.0839    1.3733    1.6623    1.9664    2.3132
   

复制代码

这个例子提示我们：复杂的方法组合，不妨就用这个鬼头鬼脑的匿名函数先封装起来，一并扔给arrayfun去处理，零售改批发。

bainhome

即使是从未接触过*fun系列函数，恐怕现在也琢磨出一点儿味道了：这一系列花里胡哨的函数，好像根本就是个循环嘛？
好问题！尤其是在关于循环的效率已经被现今的MATLAB很好解决了的背景下，提出这个问题是值得深思的，例如计算的积分：

这种非规则区域的二重积分，可以利用两次quadl函数，注意到quadl函数要求积分表达式必须写成向量化形式，所以我们构造的函数必须能接受向量输入，见如下代码：

1. function IntDemo
   
2. function f1 = myfun1(x)
   
3.      f1 = zeros(size(x));
   
4.      for k = 1:length(x)
   
5.          f1(k) = quadl(@(y) x(k)*y,sin(x(k)),cos(x(k)));
   
6.      end
   
7.      end
   
8. y = quadl(@myfun1,1,2)
   
9. end
   

复制代码

上述myfun1函数就是原始被积函数对y积分后关于x的函数，由于quadl要求其要能接受向量形式的x输入，所以写这个函数的时候考虑到x是向量的情况——里面用for语句嵌套了一个循环。如果我们利用arrayfun函数，可以将IntDemo函数精简成一句：

1. quadl(@(x) arrayfun(@(x) quadl(@(y) x*y,sin(x),cos(x)),x),1,2)

复制代码

代码变得紧凑简洁的多，更重要的是，按照这个套路可以写出一般区域的三重甚至n重循环代码。例如：

1. quadl(@(x)arrayfun(@(x)quadl(@(y)arrayfun(@(y) quadl(@ (z)x.*y.*z,x*y,2*x*y),y),x,2*x),x),1,2)

复制代码

计算的是在非矩形区域内的积分：






如果不利用arrayfun函数，还是利用子函数+通常的循环体的话，这个写起来回来会非常繁琐。上述写法具有明显的规律，可以利用递归函数求出非规则区域多重积分的表达式，甚至还可以通过先算二重，进而计算最外层积分的“2+1”的办法，代码如下：

1. quad2d(@(x,y) x.*y,1,2,@(x) sin(x),@(x) cos(x))

复制代码

继续算三重积分：

1. quadl(@(x) arrayfun(@(x) quad2d(@(y,z) x.*y.*z,x,2*x,@(y) x.*y,@(y) 2*x.*y),x),1,2)

复制代码

而这种情况下还是用子函数+循环方式写的话，繁琐程度会超乎想象。


综上，之所以有的时候不用循环而用arrayfun是为了简写代码，使代码在“形式上”可以向量化从而适应MATLAB一些向量化函数的要求。在这些情形下，用arrayfun比用循环要来的简洁明快。

bainhome

1.1 简单应用1:关于help中的example1

看完前面的例子，可能会让人觉得：“好像也蛮简单的，不就是给循环戴了副墨镜，耍耍酷什么的”。嗯…好吧，算你狠，被你看穿了。不过也别太急着下结论，其实函数学习过程中还是有些东西比较好玩儿，因为现在我们有了一些初步的“操纵句柄序列进行循环”的模糊思想(谢天谢地)，现在反过来去看help中的那个几乎令人在初识arrayfun时完全摸不着头脑的例子example1，就比较容易看出点什么了。实例如下：

1. %构造一个维数为1×15的结构数组序列
   
2. for k=1:15
   
3.    s(k).f1 = rand(k+3,k+7) * 10;            
   
4.    s(k).f2 = rand(k+3,k+7) * 10;
   
5. end

复制代码

这种构造方式很新鲜，而且可以负责任地说，几乎没有几个使用MATLAB的人会本能地采用这种构造数组的方式，它通过一重循环，对s结构变量实现扩维至1×15，最终生成一结构数组列，如图 2‑1所示：
                     

                                                   图2‑1 example1循环产生的结构数组图解
可看出：这套循环形成15个小结构数组s(i)，每个小结构数组下有两个域“f1”和“f2”，这部分内容参照struct数据类型方面的内容，略过不提。每个f1域下面有一个随机数矩阵，维数随着i的增加而变化，最后所有的小结构数组排成一行，形成一个大的结构数组s。接下来对这个生成的大结构数组进行简单处理，令这15个小结构数组中的第3、9、12个，其f2域内的随机数矩阵等于其f1内的随机数矩阵，整体排列如图 2‑2所示：

1. s(3).f2 = s(3).f1;s(9).f2 = s(9).f1;
   
2. s(12).f2 = s(12).f1;

复制代码

         

                      图2‑2 结构数组的调整
显然，数组源构建完毕，等于原料齐备，下面就对这个复杂结构数组实施加工。函数“isequal”适时登场，它对第s(1,i)(i=1,2,…,15)个子结构数组的域“f1”和“f2”进行比较，如果相等，向屏幕报告“1”，如果不等，向屏幕报告“0”，显然，我们心里有数，这是一个知道答案的过程，最后如果报告了15个数字并组成一个由0和1构成的结果数组Result，则必然是第3、9、12个数字为1，结果如图 2‑3所示。
                          
               图2‑3 isequal对每个结构数组数据进行比较的图示
代码与结果显示：

1. Result = arrayfun(@(x)isequal(x.f1, x.f2), s)  

复制代码

1. 
   
2. 
   
3. Result =
   
4.   Columns 1 through 14
   
5.      0     0     1     0     0     0     0     0     1     0     0     1     0     0
   
6.   Column 15
   
7.      0  

复制代码

结果并不出乎预料，因此不重要，重要的是：arrayfun对结构数组进行比较时，都是怎么做手脚的？例如同样的问题，如果用循环，会以如下方式实现：

1. 
   
2. for x=1:15
   
3.        Result(x)=isequal(s(1,x).f1,s(1,x).f2);
   
4. end
   
5. Result  
   

复制代码

1. 
   
2. Result =
   
3.   Columns 1 through 14
   
4.      0     0     1     0     0     0     0     0     1     0     0     1     0     0
   
5.   Column 15
   
6.      0  

复制代码

从循环体内看，s数组内部索引一直在以步长“1”的速度递增，这容易理解，就好比“f1”和“f2”是小马仔，包含他们的元素s(1,i)是小马仔他爹(比如李×)，但他爹也不是最大，外面罩着的s才是抗把子。但无论如何，用for时，到底循环谁，能看得一清二楚，每一次都必须清楚告诉for，是s部所属s(1,i)，其再下属的f1和f2实现比较；反观arrayfun内部，此时则出现了一些新鲜玩意儿：它不像我们在前面所举的任何一个实例，其循环数目“x”从字面上看似与s毫无关系！因为从help的字面意思上看arrayfun的第二个变量只是一个数组，仅此而已。这充分说明，arrayfun的手可以直接伸进结构数组s的内部，把内部的s(x).f1和s(x).f2直接从里面拽出来扔到isequal这个函数方法面前随意处置。哇！这个能量就很大了！这总能令我联想到某些事情，却模模糊糊不知道到底该联想些什么…反正不管了，总之，这个实例说明arrayfun和for还是有少许不同的。

bainhome

简单应用2:格式“不统一”的多输出

再来看一个帮助中的实例example3，不妨重新定义一个子结构数组比较少的变量s，代码如下：

1. rand('state', 0);s(1).f1 = rand(7,4) * 10;s(2).f1 = rand(3,7) * 10;
   
2. s(3).f1 = rand(5,5) * 10;

复制代码

这回s的机构被大大精简，中型部门s(i)从原来的15个变成3个，同时最小的f2部门被撤销，只剩下f1，同时f1数组维数发生变化，如图 2‑4所示，小黑旗代表该部分不存在：
                           
                           图2-4 构造新的结构数组s
下面对刚刚构成的结构数组进行LU分解[1]，代码：

1. [l u] = arrayfun(@(x)lu(x.f1), s, 'UniformOutput', false)  

复制代码

1. 
   
2. l =
   
3.     [7x4 double]    [3x3 double]    [5x5 double]
   
4. u =
   
5.     [4x4 double]    [3x7 double]    [5x5 double]  

复制代码

观察workspace发现l、u、p三个输出结果均为cell数组，而之前我们却并没有做任何初始化的工作，这说明当多输出的格式有所不同时，多个输出数据被arrayfun自动变成cell分别存储，以本次LU分解为例，每个矩阵s(i).f1做LU分解产生三个输出矩阵，依次存储在三个cell数组的第i个位置，验证代码：

1. 
   
2. ii=1:3;
   
3. q=strcat({'f'},num2str(ii','%d'));
   
4. arrayfun(@(x) assignin('base',q{x},l{x}*u{x}),1:3)
   
5. f1,f2,f3  

复制代码

1. 
   
2. f1 =
   
3.     9.5013    0.1850    1.7627    3.5287
   
4.     2.3114    8.2141    4.0571    8.1317
   
5.     6.0684    4.4470    9.3547    0.0986
   
6.     4.8598    6.1543    9.1690    1.3889
   
7.     8.9130    7.9194    4.1027    2.0277
   
8.     7.6210    9.2181    8.9365    1.9872
   
9.     4.5647    7.3821    0.5789    6.0379
   
10. f2 =
    
11.     2.7219    7.4679    4.6599    5.2515    8.3812    3.7948    7.0947
    
12.     1.9881    4.4510    4.1865    2.0265    0.1964    8.3180    4.2889
    
13.     0.1527    9.3181    8.4622    6.7214    6.8128    5.0281    3.0462
    
14. f3 =
    
15.     1.8965    1.5087    5.9356    8.1797    5.3408
    
16.     1.9343    6.9790    4.9655    6.6023    7.2711
    
17.     6.8222    3.7837    8.9977    3.4197    3.0929
    
18.     3.0276    8.6001    8.2163    2.8973    8.3850
    
19.     5.4167    8.5366    6.4491    3.4119    5.6807

复制代码

再看看变换前的矩阵s(i).f1：

1. 
   
2. s(1).f1,s(2).f1,s(3).f1  

复制代码

1. 
   
2. ans =
   
3.     9.5013    0.1850    1.7627    3.5287
   
4.     2.3114    8.2141    4.0571    8.1317
   
5.     6.0684    4.4470    9.3547    0.0986
   
6.     4.8598    6.1543    9.1690    1.3889
   
7.     8.9130    7.9194    4.1027    2.0277
   
8.     7.6210    9.2181    8.9365    1.9872
   
9.     4.5647    7.3821    0.5789    6.0379
   
10. ans =
    
11.     2.7219    7.4679    4.6599    5.2515    8.3812    3.7948    7.0947
    
12.     1.9881    4.4510    4.1865    2.0265    0.1964    8.3180    4.2889
    
13.     0.1527    9.3181    8.4622    6.7214    6.8128    5.0281    3.0462
    
14. ans =
    
15.     1.8965    1.5087    5.9356    8.1797    5.3408
    
16.     1.9343    6.9790    4.9655    6.6023    7.2711
    
17.     6.8222    3.7837    8.9977    3.4197    3.0929
    
18.     3.0276    8.6001    8.2163    2.8973    8.3850
    
19.     5.4167    8.5366    6.4491    3.4119    5.6807

复制代码

结果完全相同。
因为LU分解得到的结果全部是矩阵，而且维数可能不同，俗话说：“牛人有孤胆，怂人有群胆”，这帮货单个的时候老实得很，凑在一起就很嚣张，可不能像单个数字那样，整齐排成一队接受检阅，为了迁就它们，完整地把它们保存起来和平共处，就分门别类地把相应内容成块保存在单独的cell中。所以，一方面需要把"Uniformoutput"属性设定为'False'来达到这个目的，同时arrayfun的“多输出”，实际上是指在arrayfun的一个循环流程中输出多个变量，在下一次循环流程中，为了不把上次输出的值给冲掉，就另开一间cell把它们分别保管，这样大家都哈皮，如图 2‑6所示。现在看来，cell这名字起得还真形象！
               
                         图2-6 多个输出变量的存储方式图解




[1]
LU分解的理论内容可自行查阅数值分析教材。

bainhome

近来很多人在问：“MATLAB能否自动生成变量并相应存储动态计算得到的数据”？(这话听起来实在是别扭)，不妨用一个简单的例子说明一下：

1. for i=1:3
   
2.         str=['x',num2str(i),'=2*',num2str(i),'^2;'];
   
3.         disp(['第',num2str(i),'个执行字串:"', 'x',num2str(i),'=2*',num2str(i),'^2;"'])
   
4.         eval(str)
   
5. end
   
6. ResultArray=[x1,x2,x3]

复制代码

1. 
   
2. 第1个执行字串:"x1=2*1^2;"
   
3. 第2个执行字串:"x2=2*2^2;"
   
4. 第3个执行字串:"x3=2*3^2;"
   
5. ResultArray =
   
6.      2     8    18

复制代码

这个例子很说明问题，表现在：
-循环体内没有直接出现所希望生成的变量"x1～x3"，它们都藏在字符串str里；

-由eval把str送进MATLAB的命令执行机构，完成对每个xi的赋值；

因此，eval是把字符串变成执行命令的中转站，它不管那些字符串包含什么内容，只是负责将这些字符串扔进MATLAB和用户交流时打开的那个黑黢黢的洞口，然后自另一端就加工出各种字符串里希望产生的内容，在本例中，当然就是那一堆自行改变脚标的"x-档案"。字符串中循环中每次都是变化的，它通过num2str这个“把数字变为字串”的命令实现了和循环值i的沟通，导致每次循环时eval扔进MATLAB去处理的str都不一样，所以上述执行结果的变化也就很自然了。用数学方法来讲就是：

当然，MATLAB处理这种字符串的方法还不止一种，比如还可以这样：

1. stry=strcat('y',num2str([1:3]'));
   
2. for i=1:3
   
3.         assignin('base',stry(i,:),2*i^2)
   
4. end
   
5. ResultAssignin=[y1,y2,y3]

复制代码

1. 
   
2. ResultAssignin =
   
3.      2     8    18

复制代码

原理和效果都与前面使用eval的代码相同，只是多了strcat和assignin两个函数，strcat用于生成三个变量的名称，assignin用于在workspace中向生成的三个变量赋值。注意strcat函数构造多个变量名称的方法非常有意思，它构造字符串的思路如图 2?7所示：

图2-7 strcat合并字符串原理图解


这等于给字符串的“合并方法”做了一个类似矩阵数乘的运算，字符串矩阵的每行元素都与公共字符元素“x”合并，这是构造变量名称相对聪明的做法。而assignin函数用于给刚才构造的workspace变量进行赋值，这个函数在simulink搭建模块是常常用于simulink输入变量以及运算结果与workspace空间变量之间的眉来眼去。此方面的应用后面还会提及，此处暂时略过不提。下面主要说明*fun函数与eval、assignin函数结合使用时需要注意的问题。
显然，明眼人早已看出：这小节前面两个循环都属暗送秋波、偷笑脸红的暧昧前戏，真刀真枪的高潮肯定要落在arrayfun包装改造两个循环的狗血桥段上。
…
好吧，又被你看穿了，我们先按这个思路往下走走，试试能否弄个人人哈皮的结局：
首先对第一个循环用arrayfun改造：

1. 
   
2. arrayfun(@(i) eval(['x',num2str(i),'=2*',num2str(i),'^2;']),1:3)  

复制代码

??? Attempt to add "x1" to a static workspace.

See <a href="matlab: helpview([docroot '/techdoc/matlab_prog/matlab_prog.map'],'adding_variables_dynamically')">MATLAB Programming, Restrictions on Assigning to Variables</a> for details.

Error in ==> @(i)eval(['x',num2str(i),'=2*',num2str(i),'^2;'])

第一个例子就出问题！
完全正确的语法，却得到一长串的错误提示的红叉，那么再看看assignin的表现如何：

1. arrayfun(@(i) assignin('base',['x',num2str(i)],2*i^2),1:3)
   
2. ResultAssignin=[x1,x2,x3]

复制代码

1. 
   
2. ResultAssignin =
   
3.      2     8    18

复制代码

原理相同的代码，又得出正确结果，问题出在哪？我们先分析执行eval得到的红色错误:

??? Attempt to add "x1" to a static workspace.

这个错误的大致意思是：“静态的workspace中不能增加一个名为"x1"的变量”。但奇怪的是为什么for循环可以，而同样在使用arrayfun的前提下，循环体内用assignin函数时也能够顺利执行，独独就是eval在arrayfun中不可以呢？为回答这个问题，首先在MATLAB的command window中执行如下代码

1. arrayfun(@(i) eval(['x',num2str(i),'=2*',num2str(i),'^2;']),1:3)

复制代码

出现提示：
??? Attempt to add "x1" to a static workspace.

See MATLAB Programming, Restrictions on Assigning to Variables for details.
Error in ==> @(i)eval(['x',num2str(i),'=2*',num2str(i),'^2;'])
请注意下划线部分是关于此错误在help中的超链接，打开之后发现，使用匿名函数时，MATLAB规定了一些禁忌，其中就包括：当使用Nested function以及某些时候使用匿名函数时，所有变量只允许在函数体内部或者当前workspace中使用，禁止在运行程序时动态增加变量，这就解释了上面两个问题：
问1：为什么for循环中可以使用eval：
答1：因为for循环中我们并没有使用匿名函数或者nested function。
问2：为什么用arrayfun中用eval出问题，assignin却正常运行？
答2：因为assignin中每次给一个变量赋值，类似于在workspace中执行"a=rand(2)"这样的语句，而eval是执行字符串，而这个字符串向MATLAB提出凭空增加一个变量x1，这属于“禁忌”的非份请求，就遭到MATLAB蛮横拒绝。
对这种“eval与匿名函数齐飞”的混乱情况，help中有关"How to Avoid Using Dynamic Assignment"部分给出建议：" As a general suggestion, it is best to avoid using the eval, evalin, and assignin functions altogether."。翻译过来大致意思是：
“好汉，我说你没事儿能不能别老把eval,evalin,assignin这种函数搅和在一起乱用”？
有些“一根筋”的同学会感到奇怪：有时写函数我没用这些函数啊？怎么还会有这种错误？此时，我就会十分诚恳地告诉你：“不使用nested function时，函数调用基本是一个树状结构，也就是A调用B、B又调用C，…，这些繁琐的脏活累活，MATLAB都不会告诉您；此外，鉴于您动不动就图省事，用符号工具箱命令，syms满天飞，这些命令底层多使用eval、evalin、assignin这些行走在黑暗中的函数，自然难免问题”。
坦白讲，此类拿字符串当执行函数用的方法无论从效率还是使用规则上讲，都充满禁忌，个人意见是除非你知道你在做什么，否则没那本事能不碰最好别碰。
当然，也有"天赋异禀"的同学，一听见"禁忌"两个字，就一下"high"起来，弄得浑身上上下下躁动起来，显得很哈皮！作为年轻人，这完全可以理解，所以也有比较阴险的途径，可以变相突破这个禁忌，把这些函数，尤其是eval这个"邪恶"的函数引进来，调戏调戏MATLAB。
我们拿2.1节中曾经出现过的help实例：example1为例，它用单循环形成了1×15的两层结构数组，如果按照上述方法直接用eval，会出现“生成动态变量”的错误，那么，变通一下，可以把代码改为如下形式：

1. ii=1:15;
   
2. q=arrayfun(@(x)sprintf('s(%d).f1=10*rand(3+%d,4+%d)',x,x,x),...
   
3. ii,'UniformOutput',false);
   
4. cellfun(@eval,q)

复制代码

首先，利用arrayfun函数预先生成15个字符串，当然，这15个字符串就是定义的结构数组s，不过，我统战部卑鄙地将这15个生成字符串的数组分化、瓦解，并悄悄存在cell数组中，最后利用cellfun函数执行eval，由于此时的eval每次只进入一个cell中，并不知道其他cell的情况，于是各个击破，从而一举达到我们的邪恶目的。

图2-8 Cellfun函数中不怎么"邪恶"的eval


以同样的方法，我们可以写出前面出错例子的改正代码：

1. i=1:3;
   
2. str=arrayfun(@(x) sprintf('x%d=2*%d^2;',x,x),i,'Uniformoutput',false);
   
3. cellfun(@eval,str)
   
4. ResultCell=[x1,x2,x3]

复制代码

1. 
   
2. ResultCell =
   
3.      2     8    18

复制代码

从上述例程分析得出结论，这个动态赋值的做法让MATLAB感到很不爽，个人猜测(未必准确)原因在于：程序运行过程中动态乱加变量，是导致程序内存管理混乱、甚至程序崩溃的重要因素，打个比方，内存就是你的领地，作为君主，你当然希望你治下井井有条，人人和谐，这样才便于你大肆搜刮，可现在你的王国里经常莫名其妙以极快的速度冒出一些不明身份的人，而且每个人还都要独占一块儿地方，尽管不是每个人都会破坏安定团结，但不怕一万就怕万一，起码这种条例打开之后是隐患的开始，在这个领地里，能耍老百姓的只能是兄弟我一个帮派，多几个不认识的跟你讲你不喜欢的道理，显然不利于和谐社会大好局面的建设。所以，英明而**的MATLAB就替你做刽子手，干脆禁止这种在匿名函数和nested function里胡乱使用动态变量的定义以及赋值，于是，世界清静了。

bainhome

1.1 实例1：“参数化”绘图
这是很多人非常希望MATLAB可以轻松完成的实用功能，改变某函数的某具体参数，并绘制图形，容易令人看出相对某自变量在同取值范围下函数关系的变化走势，从而辅助完成对问题实质的把握，下面就几个简单图形的绘制，谈谈arrayfun函数的具体应用。


1.1.1 已知多组数据绘制普通2D图形
这里没有函数关系，只有多组数据，这些数据有可能是实验得到，有可能是其他仿真软件得到，有可能是MATLAB自己的数据，无论如何，在没有明确函数规律时而仅有数据的条件下，完全可以把它们分门别类地绘制出来。
例如，在模拟电子技术课程中，三极管是组成放大电路的最基础也是最重要元件之一，因此课本中作为重点中的重点，详细介绍了该元件的输入输出特性、失真方式和等效电路等各种理论内容，由于我们并不是随时随地可以用实物进行示波器测量，计算机仿真就成为该课程理论学习的必要补充手段之一，需要提起足够的重视。本例通过先期在multisim中建立三极管伏安特性仿真电路，如图 3&#8209;1所示，信号源内阻R2和集电极负载电阻R3设为可变电阻，外特性的变化就体现在这两个电阻阻值参数变化上，基极电源VBB=V1=3V，集电极电源VCC=V2=12V，基极电阻RB=R1=10kW，放置两个电流源和两个电压源测量静态工作点参数：基极电流IBQ、集电极电流ICQ、基极与发射极间电压UBEQ以及集电极与发射极间电压UCEQ，三极管放大倍数设定为80。
            
              图3-1 Multisim中的三极管放大电路
通过更改R2阻值，令IBQ依次得到[20 40 60 80 100 120]mA，打开仿真按钮进行仿真，当UCEQ为：[0 2 4 6 8]V时，测得电流源上ICQ的数值，汇总为如下数据：

1. 
   
2. VADataO =
   
3.        NaN         0    2.0000    4.0000    6.0000    8.0000
   
4.     0.0200    0.0037    1.1800    1.4650    1.7530    2.0400
   
5.     0.0400    0.0051    2.3960    2.9980    3.6020    4.2200
   
6.     0.0600    0.0054    3.6200    4.5480    5.4490    6.3710
   
7.     0.0800    0.0049    4.8470    6.0820    7.3100    8.5330
   
8.     0.1000    0.0041    5.7000    7.1590    8.6010   10.0000
   
9.     0.1200    0.0026    6.9720    8.6970   10.0000   12.0000  

复制代码

注意：这组数据为笔者在Multisim中通过仿真实验得到的测定数据，第一行为集电极与发射极间的电压值，第一列为基极电流IB，表内数据为ICQ。图形绘制代码如下：

1. 
   
2. clc;close all;clear all;
   
3. load VAData.mat
   
4. hold on;
   
5. arrayfun(@(ii) plot(VADataO(1,2:end),VADataO(ii+1,2:end),...
   
6. 'color',ii*.13*ones(1,3)),1:6)
   
7. legend('20\eta A','40\eta A','60\eta A','80\eta A','100\eta A','120\eta A')
   
8. hold off  

复制代码

                                            
                                 图3-2 三极管输出伏安特性曲线
可以看出UCE变化导致IC变化的特征，UCE=2V左右处，三极管自饱和区进入放大区，这是一个简化的伏安特性曲线。从执行代码上看出：arrayfun内部操作plot函数方法，在同一坐标轴上画图。
代码比较简单，因此没有加相应的注释，这个例子是典型的MATLAB在专业课学习过程中的应用，必须说明两点：首先，这种最基本的放大电路在MATLAB内部也可以用simulink实现，但不推荐大家用simulink仿真这种模拟电路，原因很简单：后续支持不够，它在模块数量、PCB电路板设计以及其他的专业库支持方面，还仍然远远不如Multisim、Proteus或者Altium Designer这样的软件，而且今后恐怕也肯定赶不上了，所以搞电路板设计的同仁们没必要在MATLAB上困扰太多时间；其次，在Multisim中有“伏安分析仪”这个工具，同样可以轻松实现伏安特性的仿真，这里只是为了说明当参数变化时arrayfun可以不用for循环实现图形的绘制，在绘图上增加一些代码的细节，仅此而已。故本例的先决限制范围比较窄，这两方面需要向大伙做个交待。
1.1.1 已知多组数据绘制普通3D图形
设计一个方式绘制一组平面，令所有平面交于一条直线，包括直线方程、平面方程等所有参数任意给定。



我们的示例代码如下：

1. clear;clc;close all
   
2. N=7;
   
3. colo=rand(N,3);
   
4. f=@(c)@(x,y)x+y/c;
   
5. hold on
   
6. h=arrayfun(@(x)ezmesh(f(x)),1:N);
   
7. arrayfun(@(x)set(h(x),{'FaceColor','FaceAlpha',...
   
8.     'EdgeAlpha','FaceColor'},{'interp',1,0.65,colo(x,:)}),1:N);
   
9. view(3)

复制代码

执行结果：

                                                   

                                         图3-3 经过同一直线的平面族
从代码和结果图 3&#8209;3可以看出：依然是利用arrayfun批量处理句柄的思路完成类似循环的多参数图形绘制，只是此处的函数方法变成“集束式”的多组平面，改变的参数是空间向量Z的y分量系数1/c，同时由于三维图形要先定网格，一个arrayfun搞起来不大愉快，所以弄两个，哥俩配合配合：一个骗ezmesh划网格，一个操纵前一个arrayfun生成的句柄向量设置轴属性。注意：第一个arrayfun很猛，比如用循环去完成同样的功效，需要使用如下办法：

1. 
   
2. 
   
3. N=7;
   
4. colo=rand(N,3);
   
5. f=@(c)@(x,y)x+y/c;
   
6. hold on
   
7. for x=1:N
   
8. h(x)=ezmesh(f(x));
   
9. end

复制代码

但arrayfun则根本不需要你操心到底生成了多少个h(x)，或者x到底是几，它只是把每次生成的内容，如果是单个数字，就保存成向量，如果是二维矩阵或者其他数据类型，则保存成cell数组，当然，这需要把在第2.2节提到的"Uniformoutput"属性变成"False"。
至于为什么经过同一直线，看看方程就知道，设两个任意非零常数c1和c2，构成两个平面：

由平面法向量联立得到相交直线方向向量：

得到公共直线的向量表达式(c1-c2,0,c1-c2)，此行列式可以用符号工具箱得到：

1. syms i j k c1 c2
   
2. Result=simple(det([i j k;c1 1 -c1;c2 1 -c2]))  

复制代码

1. 
   
2. Result =
   
3. (i + k)*(c1 - c2)

复制代码

y方向分量为零说明该直线落在XOZ平面上，一般方程为：

随便代个点到平面方程，既然直线落在XOZ面上，不妨一眼看出特殊值(方程很简单，所以可以看出来结果)：令y0=0，得到x0=0;z0=0。
有同学说一眼看不出怎么办？也简单，既然是直线，那么方程必然是连续函数且在所有实数范围内x和z有定义，随便令x=0，即可得到z=0，y=0.
说明直线经过原点，回到直线方程中，有x-z=0.

以上就是构造公共直线的方法，这与MATLAB无关，都是大学一年级的知识，写这些的原因是顺道儿提醒某些在论坛问出各种白痴问题却牛B哄哄如同大爷的同学，例如这个经典句型就是我长期以来的最爱：
“这么简单的问题，MATLAB也会算错”!!!
在这里小弟我谦卑、小心翼翼而且诚恳地提醒：“尊敬的老爷，计算机是工具，需要您具备些基本的背景知识，再好的武器，也需要人去瞄准某个方向，这个任务是您自己的，不要把什么事情都交给计算机去做，它只能干重复性工作，实际上计算机比您还傻(尽管您已经充分“傻B无极限、自恋无极限”了，可计算机的确还是不如你猛)”

qibbxxt

3.2.2 矩阵元素的复杂组合
题目：给出三个3*4矩阵：

1. 
   
2. A=[2 6 4 8;5 7 8 2;1 5 6 4]
   
3. B=[4 8 7 4;7 8 2 6;3 2 5 5]
   
4. C=[8 4 1 1;6 9 5 4;7 7 8 9]

复制代码

1. A =
   
2. 2 6 4 8
   
3. 5 7 8 2
   
4. 1 5 6 4
   
5. B =
   
6. 4 8 7 4
   
7. 7 8 2 6
   
8. 3 2 5 5
   
9. C =
   
10. 8 4 1 1
    
11. 6 9 5 4
    
12. 7 7 8 9
    

复制代码

希望得到下面的矩阵

1. D =
   
2. 2 0 0 6 0 0 4 0 0 8 0 0
   
3. 0 4 0 0 8 0 0 7 0 0 4 0
   
4. 0 0 8 0 0 4 0 0 1 0 0 1
   
5. 5 0 0 7 0 0 8 0 0 2 0 0
   
6. 0 7 0 0 8 0 0 2 0 0 6 0
   
7. 0 0 6 0 0 9 0 0 5 0 0 4
   
8. 1 0 0 5 0 0 6 0 0 4 0 0
   
9. 0 3 0 0 2 0 0 5 0 0 5 0
   
10. 0 0 7 0 0 7 0 0 8 0 0 9
    

复制代码

分析：首先说明排列规则，其实仔细看看很有规律，见下方D矩阵的表达形式：
1).A矩阵每个元素右面带2个0元素，变成3×12，这3行再依次填充于D矩阵的第1、4、7行；
2).B矩阵每个元素左右两侧各带一个0元素，也变成3×12，新形成的矩阵三行依次填充于D矩阵的2、5、8行；
3).C矩阵每个元素前面放2个0元素，当然也变成3×12，新形成的矩阵三行依次填充于D矩阵的3、6、9行。

向形成这样的矩阵做法很多，比如通过kron函数实现目的，先把矩阵分成三个部分：

1. 
   
2. D1=kron(A,[1 0 0;0 0 0;0 0 0])
   
3. D2=kron(B,[0 0 0;0 1 0;0 0 0])
   
4. D3=kron(C,[0 0 0;0 0 0;0 0 1])

复制代码

1. D1 =
   
2. 2 0 0 6 0 0 4 0 0 8 0 0
   
3. 0 0 0 0 0 0 0 0 0 0 0 0
   
4. 0 0 0 0 0 0 0 0 0 0 0 0
   
5. 5 0 0 7 0 0 8 0 0 2 0 0
   
6. 0 0 0 0 0 0 0 0 0 0 0 0
   
7. 0 0 0 0 0 0 0 0 0 0 0 0
   
8. 1 0 0 5 0 0 6 0 0 4 0 0
   
9. 0 0 0 0 0 0 0 0 0 0 0 0
   
10. 0 0 0 0 0 0 0 0 0 0 0 0
    
11. D2 =
    
12. 0 0 0 0 0 0 0 0 0 0 0 0
    
13. 0 4 0 0 8 0 0 7 0 0 4 0
    
14. 0 0 0 0 0 0 0 0 0 0 0 0
    
15. 0 0 0 0 0 0 0 0 0 0 0 0
    
16. 0 7 0 0 8 0 0 2 0 0 6 0
    
17. 0 0 0 0 0 0 0 0 0 0 0 0
    
18. 0 0 0 0 0 0 0 0 0 0 0 0
    
19. 0 3 0 0 2 0 0 5 0 0 5 0
    
20. 0 0 0 0 0 0 0 0 0 0 0 0
    
21. D3 =
    
22. 0 0 0 0 0 0 0 0 0 0 0 0
    
23. 0 0 0 0 0 0 0 0 0 0 0 0
    
24. 0 0 8 0 0 4 0 0 1 0 0 1
    
25. 0 0 0 0 0 0 0 0 0 0 0 0
    
26. 0 0 0 0 0 0 0 0 0 0 0 0
    
27. 0 0 6 0 0 9 0 0 5 0 0 4
    
28. 0 0 0 0 0 0 0 0 0 0 0 0
    
29. 0 0 0 0 0 0 0 0 0 0 0 0
    
30. 0 0 7 0 0 7 0 0 8 0 0 9
    

复制代码

最后实现三个矩阵的相加：

1. 
   
2. D=D1+D2+D3
   

复制代码

1. 
   
2. D =
   
3.      2     0     0     6     0     0     4     0     0     8     0     0
   
4.      0     4     0     0     8     0     0     7     0     0     4     0
   
5.      0     0     8     0     0     4     0     0     1     0     0     1
   
6.      5     0     0     7     0     0     8     0     0     2     0     0
   
7.      0     7     0     0     8     0     0     2     0     0     6     0
   
8.      0     0     6     0     0     9     0     0     5     0     0     4
   
9.      1     0     0     5     0     0     6     0     0     4     0     0
   
10.      0     3     0     0     2     0     0     5     0     0     5     0
    
11.      0     0     7     0     0     7     0     0     8     0     0     9
    

复制代码

可以看出，这些矩阵形成时的代码有类似之处，因此可以用如下方式实现：

1. 
   
2. d={A,B,C};
   
3. e=eye(length(d));
   
4. f=cell(1,1,length(d));
   
5. for i=1:length(d)
   
6.     f{i}=kron(d{i},diag(e(i,:)));
   
7. end
   
8. g=cell2mat(f);
   
9. D=sum(g,3)
   
10. 
    

复制代码

1. 
   
2. D =
   
3.      2     0     0     6     0     0     4     0     0     8     0     0
   
4.      0     4     0     0     8     0     0     7     0     0     4     0
   
5.      0     0     8     0     0     4     0     0     1     0     0     1
   
6.      5     0     0     7     0     0     8     0     0     2     0     0
   
7.      0     7     0     0     8     0     0     2     0     0     6     0
   
8.      0     0     6     0     0     9     0     0     5     0     0     4
   
9.      1     0     0     5     0     0     6     0     0     4     0     0
   
10.      0     3     0     0     2     0     0     5     0     0     5     0
    
11.      0     0     7     0     0     7     0     0     8     0     0     9
    
12. 
    

复制代码

同样，能够很方便地用arrayfun改写上述代码：

1. 
   
2. d=cat(3,A,B,C);
   
3. N=size(d,3);
   
4. e=eye(N);
   
5. f=arrayfun(@(i)kron(d(:,:,i),diag(e(i,:))),1:N,'UniformOutput',false);
   
6. D=sum(cell2mat(reshape(f,[1,1,N])),3)  
   
7. 
   

复制代码

1. 
   
2. D =
   
3.      2     0     0     6     0     0     4     0     0     8     0     0
   
4.      0     4     0     0     8     0     0     7     0     0     4     0
   
5.      0     0     8     0     0     4     0     0     1     0     0     1
   
6.      5     0     0     7     0     0     8     0     0     2     0     0
   
7.      0     7     0     0     8     0     0     2     0     0     6     0
   
8.      0     0     6     0     0     9     0     0     5     0     0     4
   
9.      1     0     0     5     0     0     6     0     0     4     0     0
   
10.      0     3     0     0     2     0     0     5     0     0     5     0
    
11.      0     0     7     0     0     7     0     0     8     0     0     9
    
12. 
    

复制代码

这样，利用arrayfun，替代了对kron函数的循环操作。当然，这种做法也不是唯一的，个人给出一个“矩阵味道”更浓的做法，即：行列变换

1. 
   
2. A(9,12)=0;                   %A矩阵扩维
   
3. B(9,12)=0;                   %B矩阵扩维
   
4. C(9,12)=0;                   %C矩阵扩维
   
5. D=A([1,4,5,2,6,7,3,8,9],[1,5,6,2,7,8,3,9,10,4,11,12])+...
   
6.   B([4,1,5,6,2,7,8,3,9],[5,1,6,7,2,8,9,3,10,11,4,12])+...
   
7.   C([4,5,1,6,7,2,8,9,3],[5,6,1,7,8,2,9,10,3,11,12,4])
   

复制代码

1. 
   
2. D =
   
3.      2     0     0     6     0     0     4     0     0     8     0     0
   
4.      0     4     0     0     8     0     0     7     0     0     4     0
   
5.      0     0     8     0     0     4     0     0     1     0     0     1
   
6.      5     0     0     7     0     0     8     0     0     2     0     0
   
7.      0     7     0     0     8     0     0     2     0     0     6     0
   
8.      0     0     6     0     0     9     0     0     5     0     0     4
   
9.      1     0     0     5     0     0     6     0     0     4     0     0
   
10.      0     3     0     0     2     0     0     5     0     0     5     0
    
11.      0     0     7     0     0     7     0     0     8     0     0     9
    
12. 
    

复制代码

总结：在矩阵的操作中arrayfun的实现方法往往不占优势，甚至还要牺牲很多效率，这从另一个角度说明这个函数的真正用意还是在函数句柄，而不是具体数字的批量操作，所以绝非循环的泛泛替代品。不过经过一系列函数的复杂套用，arrayfun基本也能够实现最终目的，说明这个函数与矩阵的基本操作之间同样有藕断丝连的瓜葛，往往在操作数字的时候，底层调用的依旧是那些基本的函数，只是在外面包了一个匿名函数的壳而已。通过最后给出的替代求解方法，有理由认为：当出现矩阵操作时，被很多人忽略了的最基本和最简单的命令往往就是解决问题的王道。
关于此题的更多做法参考：http://forum.simwe.com/thread-963252-1-2.html

rocwoods

上面分析了这么多，我想大家对*fun系列函数有了比较清楚的认识，下面再来一个arrayfun另类应用但是用循环很难实现的例子——求一般区域上的n重积分。事先说明一点的是虽然这里讨论的n理论上讲不限大小，但是鉴于n重积分的特点，如果积分重数高，同时积分区域范围较大，数值积分计算量十分巨大，如果各重积分原函数都能显式表示，对于多重积分，采用符号积分是第一选择。否则再只能采用数值计算情况下，多重积分（大于3重）一般是采用蒙特卡洛法近似计算。如果对精度有要求但对时间没有要求的情况下，我们可以借助arrayfun函数利用递归的形式把普通的求一、二重积分的数值函数调用起来求多重积分。
废话少说，下面介绍求解一般区域n重积分的函数——nIntegrate。
该函数的调用语法如下：

f = nIntegrate(fun,low,up)     

f为函数的返回值是n重积分积分结果。
fun是被积函数字符串形式，不同的变量依次以x1,x2,...xn表示，（需要注意的是，必须以x1,x2,...xn这种形式表示，这是由函数内部递归构造求解表达式决定的，其余像y1,y2,...yn或是其他表示方法都不行）。
low和up都是长度为n的cell数组，low存储从外到内各重积分的积分下限函数，up存储从外到内各重积分的积分上限函数（都是字符串形式）。low和up内的函数表示都要遵循一些原则，这些原则在程序注释里进行了说明，后续例子中将进一步详细说明。
下面给出nIntegrate函数的代码。

1. 
   
2. function f = nIntegrate(fun,low,up)
   
3. %f, 返回值，n重积分积分结果
   
4. %fun, 是被积函数字符串形式；
   
5. %low存储从外到内各重积分的积分下限函数；
   
6. %up存储从外到内各重积分的积分上限函数（都是字符串形式）
   
7. %为了保证函数正常运行，low和up内的函数遵循如下原则：设积分重数为m，最内层到最外层的
   
8. %积分变量依次以xm,...x2,x1来表示（只能以这样顺序，其他顺序或者别的字母表示变量都不可以）
   
9. %最内层的上下限函数不管是不是关于x(m-1)的函数都要求x(m-1)必须出现，不是其函数时
   
10. %可以写成“+0*x(m-1)”等形式,依次类推次内层要求x(m-2)必须要出现等等，一直到次外层
    
11. %要求变量x1出现，最外层才是常数。
    
12. N = length(low); %判断积分重数
    
13. fun = vectorize(fun); %将被积函数写成点乘形式，方便数值积分函数调用
    
14. if verLessThan('MATLAB','7.8') %判断当前运行该函数的MATLAB版本是否低于7.8即R2009a
    
15.     %低于7.8调用GenerateExpr_quadl函数递归构造求解表达式
    
16.     expr = GenerateExpr_quadl(N);
    
17. else %7.8或者以后的版本调用GenerateExpr_quad2d函数递归构造求解表达式
    
18.     if mod(N,2) == 0
    
19.         expr = GenerateExpr_quad2d(N);
    
20.     else
    
21.         expr = ['quadl(@(x1) arrayfun(@(x1)',...
    
22.             GenerateExpr_quad2d(N-1),',x1),',low{1},',',up{1},')'];
    
23.     end
    
24.    
    
25. end
    
26. %=======================================================
    
27. %主要利用quadl函数递归构造求解表达式，适用于R2009a以前的版本
    
28. %=======================================================
    
29.     function expr = GenerateExpr_quadl(n)
    
30.         if n == 1
    
31.             expr = ['quadl(@(x',num2str(N),')',fun,',',low{N},',',up{N},')' ];
    
32.         else
    
33.             expr = ['quadl(@(x',num2str(N-n+1),') arrayfun(@(x',...
    
34.               num2str(N-n+1),')',GenerateExpr_quadl(n-1),',x',num2str(N-n+1),...
    
35.                 '),',low{N-n+1},',',up{N-n+1},')'];
    
36.         end
    
37.     end
    
38. %============================================================
    
39. %主要利用quad2d函数递归构造求解表达式，适用于R2009a以及以后的版本
    
40. %============================================================
    
41.     function expr = GenerateExpr_quad2d(n)
    
42.         if n == 2
    
43.             expr = ['quad2d(@(x',num2str(N-1),',x',num2str(N),')',fun,',',...
    
44.                 low{N-1},',',up{N-1},',@(x',num2str(N-1),')',low{N},',@(x',...
    
45.                 num2str(N-1),')',up{N},')' ];
    
46.         else
    
47.             expr = ['quad2d(@(x',num2str(N-n+1),',x',num2str(N-n+2),')',...
    
48.                 'arrayfun(@(x',num2str(N-n+1),',x',num2str(N-n+2),')',...
    
49.                 GenerateExpr_quad2d(n-2),...
    
50.                 ',x',num2str(N-n+1),',x',num2str(N-n+2),'),',...
    
51.                 low{N-n+1},',',up{N-n+1},...
    
52.                 ',@(x',num2str(N-n+1),')',low{N-n+2},...
    
53.                 ',@(x',num2str(N-n+1),')',up{N-n+2},')' ];
    
54.         end
    
55.     end
    
56. f = eval(expr);
    
57. end
    
58. 
    

复制代码

分析nIntegrate函数的实现代码，大家会发现，这个函数就是利用arrayfun函数封装循环的特性，借助递归方法把MATLAB自带的quadl、quad2d等计算一重和二重积分的函数串联起来形成计算多重积分的字符串表达式（这是和arrayfun遍历循环高度简洁化、模板化分不开的），然后利用eval函数执行这个表达式来得到结果。
怎么样？这个用法很变态吧？如果采用普通循环来写出递归表达式恐怕不是很容易。

bainhome

既然前面第3.2节说arrayfun处理数据时有些不爽，加之使用又没有循环那么直观，正确的做法又该是怎样的呢？好问题！此时可以参看前面第2.1和第2.2节给出了help中关于arrayfun的“正宗嫡传”心法，这对我们准确把握“如何使用arrayfun系列函数”给出了更多的提示与启发，于是借用胡寿松《自动控制原理》第5版教材第3章：线性系统的时域分析法习题3-26的求解，用arrayfun做出一些调整，虽然里面多数内容可以通过for替代，但你会发现，用arrayfun等命令去做控制类的题目，也有很多便利的优点。

题目：喷气式战斗机自动驾驶仪中配有横滚控制系统，结构图如图 3&#8209;4所示，要求：
(1).确定闭环传递函数
(2).当K1分别等于0.7，3.0和6.0时，确定闭环系统特征根；
(3).在(2)所给条件下，应用主导极点概念，确定各二阶近似系统，估计原有系统的超调量和峰值时间；

(4).绘制原有系统的实际单位阶跃响应曲线，并与(3)中的近似结果进行比较

bainhome

题目本身是为练习主导极点的确定与应用，此处援引目的则是突出arrayfun函数的应用。因此前面部分与arrayfun无关内容快速通过，感兴趣者可参考胡寿松教材自行学习。
(1).闭环传递函数

1. 
   
2. den1=arrayfun(@(x) [1 10],1:3,'uniformoutput',false);
   
3. den2=arrayfun(@(x) [1 1.4 0],1:3,'uniformoutput',false);
   
4. num1={.7,3,6};num2={11.4,11.4,11.4};
   
5. sys1=tf(num1,den1);
   
6. sys2=tf(num2,den2);
   
7. for i=1:3
   
8. sys{i}=feedback(series(sys1(i),sys2(i)),1);
   
9. end
   
10. sys{:}  

复制代码

1. 
   
2. Transfer function:
   
3.             7.98
   
4. ----------------------------
   
5. s^3 + 11.4 s^2 + 14 s + 7.98
   
6.   
   
7. Transfer function:
   
8.             34.2
   
9. ----------------------------
   
10. s^3 + 11.4 s^2 + 14 s + 34.2
    
11. 
    
12. Transfer function:
    
13.             68.4
    
14. ----------------------------
    
15. s^3 + 11.4 s^2 + 14 s + 68.4
    

复制代码

点评:之所以在上述代码中求解系统开环增益不同时的闭环传递函数时，选择使用循环而没有用arrayfun，原因是arrayfun函数当前并不支持"tf"这种数据类型，所以从另一个角度看出：这个函数目前与MATLAB内部的其他函数之间相处也不像想象中那么融洽无间。
(2).闭环特征根

1. 
   
2. [Num,Den]=arrayfun(@(x) tfdata(sys{x}),1:3,'uniformoutput',false);
   
3. Root=arrayfun(@(x) roots(Den{x}{1}),1:3,'uniformoutput',false);
   
4. celldisp(Root)  
   

复制代码

1. 
   
2. Root{1} =
   
3. -10.0910         
   
4.   -0.6545 + 0.6020i
   
5.   -0.6545 - 0.6020i
   
6. Root{2} =
   
7. -10.3678         
   
8.   -0.5161 + 1.7414i
   
9.   -0.5161 - 1.7414i
   
10. Root{3} =
    
11. -10.6889         
    
12.   -0.3555 + 2.5045i
    
13.   -0.3555 - 2.5045i  
    

复制代码

(3).二阶近似系统及其动态性能
分析：稳定的高阶系统，其闭环极点和零点在左半s平面上虽有各种分布模式，但就与虚轴之间距离来说，只有远近之分。如果在所有闭环极点中，距离虚轴最近的极点周围无闭环零点，而其他闭环极点又远离虚轴，那么距离虚轴最近的闭环极点所对应相应分量随时间推移衰减缓慢，在系统的时间响应过程中起主导作用，这样的闭环极点称为闭环主导极点。在控制工程中，通常要求控制系统既具有较快的响应速度，又具有一定的阻尼程度，此外，还要求减少死区、间隙和库仑摩擦等非线性因素对系统性能的影响，因此高阶系统的增益常常调整为使系统具有一对闭环共轭主导极点。这时可以利用二阶系统的动态性能指标估算高阶系统的动态性能 。

以上述系统为例，当 K1=0.7时,绘制系统零极点图形,如下:

1. pzmap(sys{1})

复制代码

结合(2)中的闭环特征根求解，容易看出一个单实根s1=-10.091,实部模比主导极点大很多，因此是一对闭环共轭主导极点，故该三阶系统近似为如下系统：


同理，当K1=3.0和K1=6.0时的近似系统二阶闭环传递函数分别为：


和：

值得新学习自动控制原理的同学注意的是：近似二阶系统需要写成标准形式，在与原系统比较中容易发现其分子已经发生了变化，在分子上的项是近似系统的开环增益，这是容易理解的：
以为例，非主导极点剔除之后系统特征多项式：




结合二阶系统闭环传递函数标准形式：




得到：




另外两个系统原理同上。


动态性能中列出延迟时间td，峰值时间tp，超调量和调节时间ts，由于参量较多，不妨写成结构数组便于调用。
步骤1：根据近似系统传递函数，容易写出动态特性参数计算代码如下：

1. 
   
2. NewNum={0.7908,3.299,6.399};
   
3. NewDen={[11.309 .7908],[1 1.032 3.299],[1 .7111 6.399]};
   
4. wn=sqrt([NewDen{1}(3),NewDen{2}(3),NewDen{3}(3)]);
   
5. % 二阶系统自然频率
   
6. ksi=[NewDen{1}(2),NewDen{2}(2),NewDen{3}(2)]/2./wn;
   
7. % 阻尼比
   
8. sigma=100*exp(-pi*ksi./sqrt(1-ksi.^2));
   
9. % 超调量
   
10. tp=pi./wn./sqrt(1-ksi.^2);
    
11. % 峰值时间
    
12. td=(1+.6*ksi+.2*ksi.^2)./wn;
    
13. % 延迟时间
    
14. ts=3.5./ksi./wn;
    
15. % 调节时间
    
16. % 组合为结构数组
    
17. Struct=struct('K1',[.73 6],'wn',wn,'ksi',ksi,'sigma',sigma,'tp',tp,'td',td,'ts',ts)
    

复制代码

1. 
   
2. Struct=
   
3. K1: [0.7000 3 6]
   
4. wn: [0.88931.8163 2.5296]
   
5. ksi: [0.73600.2841 0.1406]
   
6. sigma: [3.2862 39.421764.0189]
   
7. tp: [5.21841.8040 1.2544]
   
8. td: [1.74290.6533 0.4302]
   
9. ts: [5.34766.7829 9.8439]
   

复制代码

(4).单位阶跃响应曲线
学过自动控制的同学都知道，我们用step函数来绘制LTI系统的单位阶跃响应曲线。当然还有其他方法，例如用反laplace变换把施加阶跃响应得到的传递函数写成时域形式，再用plot绘制，这种方式直观经典，但不新鲜：网上有很多人早就给过现成示例，并做成手写GUI的范例广泛传播，因此这里不再赘述。现在不做反laplace变换，仅使用step这个现成函数，结合arrayfun之后，也可以变化出几种具有启发性的代码。利用(1)中不同sys,可绘制不同K1条件下三阶系统阶跃响应曲线，为比较在(3)中的剔除非主导极点近似二阶系统，把它们的阶跃响应也一并绘制在同一张图上，代码如下：

1. ApproxNum={0.7908,3.299,6.399};
   
2. ApproxDen={[1 1.309 .7908],[1 1.032 3.299],[1 .7111 6.399]};
   
3. for x=1:3
   
4.        ApproxSys{x}=tf(ApproxNum{x},ApproxDen{x});
   
5. end
   
6. hold on;
   
7. step(ApproxSys{1},'b-.',ApproxSys{2},'g-.', ApproxSys{3},'r-.')
   
8. arrayfun(@(x) step(sys{x}),1:3)
   
9. legend('K_1=0.7','K_1=3.0','K_1=6.0','K_1=0.7','K_1=3.0','K_1=6.0')
   
10. set(gca,'fontsize',14)
    
11. h=findobj(gca,'type','line');
    
12. set(h,'linewidth',2)
    
13. hold off

复制代码

bainhome

可看出二阶近似系统与三阶系统的阶跃响应曲线虽然不重合，但也已经非常接近。
如果不想把不同取值条件下的曲线绘制在同一张图上，也有很多办法，如下给出利用arrayfun的方式，由于绘制多幅图形时需要每次新开一个figure，因此仅仅利用arrayfun去操纵step是不方便的，不过arrayfun对于句柄的胃口很好，不妨单另写一个function，然后让arrayfun操纵这个函数即可，步骤如下：
步骤1：新开一个M函数，保存如下代码：

1. 
   
2. function MultiSysStep(sys1,sys2)
   
3. figure;
   
4. step(sys1,'r',sys2,'b-.')
   
5. set(gca,'fontsize',14)
   
6. h=findobj(gca,'type','line');
   
7. set(h,'linewidth',2)  
   

复制代码

步骤2：在command windows中执行如下代码：

1. arrayfun(@(x) MultiSysStep(sys{x},ApproxSys{x}),1:3)  
   

复制代码

可以看出，保存一个function，然后用arrayfun来调用，这是一种比较灵活的形式，有点儿类似于Windows中的批处理命令。
现在撇开arrayfun的应用，谈点儿关于step的题外话：假如还想把不同 取值下的系统以子图形式进行保存，该如何做？其实这更加简单，因为step和tf命令本身就具备矢量操作的功能，只需把传递函数数据形式略做改变，代码如下：

1. ApproxNum={0.7908,3.299,6.399};
   
2. ApproxDen={[1 1.309 .7908],[1 1.032 3.299],[1 .7111 6.399]};
   
3. sysApprox=tf(ApproxNum,ApproxDen);
   
4. step(sysApprox)
   

复制代码

既然已经到这一步了，我们不妨再深入一些：假如想在图 3&#8209;7基础上再做点儿改进：把原系统和近似系统分别放在三幅子图中比较，又该怎么做呢？一般来说这样分组绘图的问题直接利用step绘图是有困难的，好在这个命令的调用格式灵活而强大，我们可以先用step存储图形相关数据，再利用subplot函数接力画图。
步骤1：用step存储数据但并不画图，代码如下：

1. 
   
2. 
   
3. [yOri,tOri]=arrayfun(@(x) step(sys{x}),1:3,'uniformoutput',false)
   
4. [yApp,tApp]=arrayfun(@(x) step(ApproxSys{x}),1:3,'uniformoutput',false)  
   
5. 
   

复制代码

1. 
   
2. 
   
3. yOri
   
4. =
   
5. 
   
6.     [127x1 double]    [139x1 double]    [391x1 double]
   
7. 
   
8. tOri =
   
9. 
   
10.     [127x1 double]    [139x1 double]    [391x1 double]
    
11. 
    
12. yApp =
    
13. 
    
14.     [127x1 double]    [104x1 double]    [127x1 double]
    
15. 
    
16. tApp =
    
17. 
    
18.     [127x1 double]    [104x1 double]    [127x1 double]  
    
19. 
    

复制代码

步骤2：所生成数据用subplot绘图：

1. for i=1:3
   
2. subplot(1,3,i);
   
3. plot(tApp{i},yApp{i},'b-.',tOri{i},yOri{i},'r')
   
4. end

复制代码

总结：从上述例子可以看出，arrayfun函数在控制系统分析的学习过程中很有作用，它既可以帮助直观而快速地表达调整某个环节参数对于系统本身的影响，加上输出参数则可以自动生成cell数据，与step、tf等命令能方便地结合，在很大程度上可以替代循环，使代码更简洁明了，同时从上述不同方式绘制阶跃响应图的过程也能发现：合理利用struct、cell、double等不同的数据类型，挖掘命令自身的潜力(比如step函数的不同调用格式可以生成阶跃图形，也可以用来提取绘图数据)，并把二者组合起来，对深入了解MATLAB的精髓很有好处。

bainhome

cellfun和structfun函数和arrayfun非常类似，了解了arrayfun，这两个函数就十分好理解，前者对cell内部每个单元，利用指定的函数方法挨个点名，后者则遍历结构数组中的每个域内的数据实现处理。
cellfun函数无论从其调用格式还是所给的实例来看，都是非常典型的"arrayfun系"成员。比如先看最常用的调用格式：

A = cellfun(fun, C, D, ...)

其中：
fun:函数句柄，可以是匿名函数、自定义函数、官方函数等等；
C,D:需要用fun处理的cell数组；
换通俗的话讲就是：cellfun是抗把子，令fun深入基层，与每个cell内部的群众亲密接触，而且cellfun是个很仗义的老大，不但命令fun和与群众接触，而且要有反馈结果，于是就有了A这个存储对不同群众反馈意见的结果箱。又或者邪恶地想：cellfun是“庸俗的金子”，fun是韦小宝，C,D这些cell好比教主夫人、曾柔、公主、双儿…，而A是韦虎头、韦铜锤…由于不可抗力因素cellfun的出现，当fun遇到C,D…想必是要“发生”一些事情，或者“发生”一些人。下面不妨从MATLAB，而不是生理卫生的角度简要分析一下，“发生”的原理：
不妨再次以第节中的控制系统阶跃响应图绘制为例，按照其(1)解答中那样构造sys的cell数据，然后绘图，代码如下：

1. 
   
2. 
   
3. den1=arrayfun(@(x) [1 10],1:3,'uniformoutput',false);
   
4. den2=arrayfun(@(x) [1 1.4 0],1:3,'uniformoutput',false);
   
5. num1={.7,3,6};num2={11.4,11.4,11.4};
   
6. sys1=tf(num1,den1);
   
7. sys2=tf(num2,den2);
   
8. for i=1:3
   
9. sys{i}=feedback(series(sys1(i),sys2(i)),1);
   
10. end
    
11. ApproxNum={0.7908,3.299,6.399};
    
12. ApproxDen={[1 1.309 .7908],[1 1.032 3.299],[1 .7111 6.399]};
    
13. for x=1:3
    
14.        ApproxSys{x}=tf(ApproxNum{x},ApproxDen{x});
    
15. end
    
16. hold on
    
17. cellfun(@step,sys,ApproxSys)
    
18. set(gca,'fontsize',14)
    
19. h=findobj(gca,'type','line');
    
20. set(h,'linewidth',2)

复制代码

图 4&#8209;1达到了和前面图 3&#8209;6一样的效果，这说明反馈的事情已经“发生”了，但是人还没有“发明”出来，于是再看如下代码：

1. 
   
2. 
   
3. [yOri,tOri]=cellfun(@step,sys,'uniformoutput',false)
   
4. 
   
5. [yApp,tApp]=cellfun(@step,ApproxSys,'uniformoutput',false)  

复制代码

1. yOri=
   
2.     [127x1 double]    [139x1 double]    [391x1 double]
   
3. tOri =
   
4.     [127x1 double]    [139x1 double]    [391x1 double]
   
5. yApp =
   
6.     [127x1 double]    [104x1 double]    [127x1 double]
   
7. tApp =
   
8.     [127x1 double]    [104x1 double]    [127x1 double]  

复制代码

有人会问，你干嘛不像前面一样，把sys和ApproxSys这两个cell写到同一个cellfun里得到全部数据？不妨试试：

1. [yApp,tApp]=
   
2. cellfun(@step,sys,ApproxSys,'uniformoutput',false)  
   

复制代码

1. ??? Error using ==> DynamicSystem.step at 60
   
2. The step" command operates on a single model when used with output arguments.  
   

复制代码

这段出错信息翻译过来大致意思是：当需要输出结果时，step命令一回只能操纵单一模型，饭一口一口吃，一般别把媳妇和小三名字混在一起乱喊，太危险。之所以说仅仅是“危险”而非禁忌，是因为某些函数也还允许被搞得比较荒唐：比如help中给出的如下求协方差例子：

1. C = {randn(10,1), randn(20,1), randn(30,1)};
   
2. D = {randn(10,1), randn(20,1), randn(30,1)};
   
3. CDcovs = cellfun(@cov, C, D, 'UniformOutput', false)  
   

复制代码

1. CDcovs =     [2x2 double]    [2x2 double]    [2x2 double]

复制代码

当然，熟悉cov的同仁很容易发现：这种情况需要两个前提条件才能产生正确结果：
1).cov本身根本就需要两个参数；
2).C和D必须同维
这样，每次cellfun命令cov同时从C中取C{i}、D中去D{i}，“发生”一个结果，放入CDcovs{i}。
另外，还可以看看cellfun的其他命令调用格式，因为cell的神秘性，MATLAB特意给cellfun开了个小灶，给出了如下命令调用方式：

1. A =cellfun('fname', C)

复制代码

这里的'fname'可以是如下表中的任意一个：
函数                          返回值
isempty        true for an empty cell element
islogical        true for a logical cell element
isreal            true for a real cell element
length           Length of the cell element
ndims           Number of dimensions of the cell element
prodofsize     Number of elements in the cell element

例如'length'可以得到所操纵cell的每个子cell内数据的列维数：

1. 
   
2. C={[1:10],rand(3,4),[],['abc','de']}
   
3. A=cellfun('length',C)
   

复制代码

1. C =
   
2.     [1x10 double]    [3x4 double]    []    'abcde'
   
3. A =
   
4.     10     4     0     5  

复制代码

4.2 structfun函数


有了前面两个函数的介绍，这部分内容基本已经不用多说，原理、数据存贮、函数调用方法…所有内容都相同，只不过cellfun在cell中遍历每个细胞的数据，而structfun则专找strcut结构数组的麻烦，以help中的例子做简要说明：

1. s.f1 = 'Sunday';
   
2. s.f2 = 'Monday';
   
3. s.f3 = 'Tuesday';
   
4. s.f4 = 'Wednesday';
   
5. s.f5 = 'Thursday';
   
6. s.f6 = 'Friday';
   
7. s.f7 = 'Saturday';
   
8. 
   
9. shortNames = structfun(@(x) (x(1:3)),s,'UniformOutput',false)  

复制代码

1. shortNames =
   
2.     f1: 'Sun'
   
3.     f2: 'Mon'
   
4.     f3: 'Tue'
   
5.     f4: 'Wed'
   
6.     f5: 'Thu'
   
7.     f6: 'Fri'
   
8.     f7: 'Sat'

复制代码

上述代码对于遍历s结构数组治下每个'field'内的字符串，取得前三个字符，即：三字母的缩写。假如带点儿想象力去看待这个函数，会先想到如下两个问题
1).遍历struct数组的函数生成struct类型结果数组，类型没变，但结果没有通用性；
2).把函数方法和struct内部的数据换一换，功能和cellfun和structfun相同。
不妨再拿前面的阶跃响应曲线绘制试试手，为便于理解，借鉴前面得到的传递函数结果，生成一个合适的结构数组供step调用，代码如下：

1. Num={7.98,34.2,68.4};
   
2. Den={[1 11.4 14 7.98],[1 11.4 14 34.2],[1 11.4 14 68.4]};
   
3. sys=tf(Num,Den);
   
4. ApproxNum={0.7908,3.299,6.399};
   
5. ApproxDen={[1 1.309 .7908],[1 1.032 3.299],[1 .7111 6.399]};
   
6. sysApprox=tf(ApproxNum,ApproxDen);
   
7. structSys=struct('sys',sys,'AppSys',sysApprox);
   
8. [SResultsY,SResultsT]=structfun(@step,structSys,'uniformoutput',false)
   

复制代码

1. SResultsY =
   
2.        sys: [391x1x3 double]
   
3.     AppSys: [196x1x3 double]
   
4. SResultsT =
   
5.        sys: [391x1 double]
   
6.     AppSys: [196x1 double]

复制代码

结果可以用plot来验证：

1. hold on;
   
2. arrayfun(@(x) plot(SResultsT.sys,SResultsY.sys(:,:,x)),1:3)
   
3. arrayfun(@(x) plot(SResultsT.AppSys,SResultsY.AppSys(:,:,x),'r-.'),1:3)
   
4. set(gca,'fontsize',14)
   
5. h=findobj(gca,'type','line');
   
6. set(h,'linewidth',2)
   
7. axis tight

复制代码

bainhome

经过上述内容的介绍，相比对于arrayfun、cellfun和structfun函数会有一个简单的了解，其实客观的说，这系列函数如果抛开函数句柄操作的内容不讲，的确有点儿类似for循环失散多年的孪生小兄弟，只不过集成化更好，例如可以不加循环次数，在cell等数据内部自行遍历；又或者输出数据自行存储等，都能够说明这个函数比他的"for"兄弟要聪明点儿。同时由于对句柄操纵不需要认为干涉太多，使得代码从表面看更加简单清爽，在一些特殊的地方，例如前述的多重积分、控制系统多参数绘图等，使用熟练后，这些函数都比循环更利于输出数据的管理和统一，当然用循环一样能完成这些工作，但在代码控制上就需要编程者给出更多的干预。
除了优点，这系列函数的缺点也是显而易见的，其首恶就是效率，在坛子里rocwoods已经给出了一个分析的帖子，它的效率比之循环在几乎所有函数操纵过程中都呈倍数降低，甚至更高。这是可以理解的，无论是数据的存储，句柄内容等等，都要比循环复杂得多。这点倒是有点儿类似winsows和lunix之间的区别，别人替你考虑的多、你自己手动少，必然带来效率的降低。这是一个也许很不恰当的联想，不过实在找不出更合适的比喻，将就凑合这个吧。
其次，代码没有for直观，因为循环被隐藏在函数内部，使得用户一眼看去，这些函数没有它兄弟“for”那么被待见，我也一直比较喜欢说：“如果一个MATLAB的命令你看help十分钟没弄明白它到底在干什么,就不需要弄明白:想法子用简单而容易理解的命令去代替或者部分代替它”。但这个"十分钟"随着水平的提高,会从MATLAB的help内看出日益丰富的内容,并非一成不变,因此当你对于MATLAB数据的结构相对有些理解之后,arrayfun有时不大直观的问题应该会迎刃而解。
最后，特殊的函数造成被特殊地理解，由于和句柄函数有着千丝万缕的瓜葛，某些数据类型如tf等和很多函数如eval、evalin、assignin等之间，现在配合起来还是不大默契。
================================================================

现在，目前我们所写内容已全部发出

ps：前面有个会员写了个参考的evalin做法，很不错，欢迎大家继续讨论。

taohe

这个帖子的楼上各位都很给力啊，这是个很好的新年礼物。谢谢。

关于匿名函数等和原先for循环代码的执行速度差异，或许可以这样理解：从matlab6.5起，matlab使用了JIT的优化、加速技术，使得for循环的执行得以大幅加速。而对于匿名函数等的加速效果可能太弱。

对于其他普通形式的运算，比如矩阵相乘等，使用匿名函数的方式和直接编写的代码的执行速度差别应该不大，请参见Loren的以下博文：
http://blogs.mathworks.com/loren ... of-function-inputs/

qibbxxt

17# taohe
恩，今天刚刚看了这篇文章
自从用arrayfun一来一直在思考速度的问题
上次运行8皇后，用for只用5s，但用arrayfun却用了23s，除了程序自身设计的问题外，看来匿名函数和arrayfun的速度还是有点慢
我觉得tao兄说的有理

rocwoods

纠正一个错误，我之前关于arrayfun效率比循环低很多的结论很草率。
就自身的效率而言，arrayfun应该是不输循环，甚至还要好些。http://forum.simwe.com/thread-965287-1-3.html 这个帖子对循环和arrayfun的比较不是同等条件的。因为arrayfun无可避免的要和函数句柄绑定在一起，arrayfun的开销中有相当部分是函数句柄调用开销。正确的和循环的PK应该如下：

1. 
   
2. >> f = @(x) x^2;
   
3. >> x = 1:1000000;
   
4. >> tic;z1 = arrayfun(f,x);toc
   
5. Elapsed time is 10.925279 seconds.
   
6. >> z2 = zeros(size(z1));
   
7. >> tic;for k =1:1000000;z2(k) = f(x(k));end;toc
   
8. Elapsed time is 13.356582 seconds.
   
9. >> isequal(z1,z2)
   
10. ans =
    
11.      1
    

复制代码

从上面可以看出，同样调用100万次同一个函数句柄，arrayfun要比for循环快20%左右。说实在的，匿名函数效率已经算是比较好的了，（虽然和直接调用MATLAB的加减乘除等built-in函数还是有巨大差距，见http://forum.simwe.com/thread-965287-1-3.html 中的比较）如果换成内联函数呢？为什么我们现在强烈倡导要摒弃内联函数，下面还是这个例子，大家看看内联函数的效率有多低。由于arrayfun不直接支持内联函数，所以需要改动一下：

1. 
   
2. >> clear
   
3. >> g = inline('x^2');
   
4. >> f = @(x) g(x);
   
5. >> x = 1:10000;
   
6. >> tic;z1 = arrayfun(f,x);toc
   
7. Elapsed time is 2.048061 seconds.
   
8. >> z2 = zeros(size(z1));
   
9. >> tic;for k =1:10000;z2(k) = f(x(k));end;toc
   
10. Elapsed time is 2.071208 seconds.
    

复制代码

仅仅调用1万次，时间已经超过2秒了。不过这个例子也同样说明了arrayfun本身的效率是不输循环的。只不过被匿名函数调用开销脱了后腿，导致误认为不如循环效率高。

bainhome

inline要在字符串、eval、句柄之间来回捣鼓，使用的过程基本都是效率毒药。
发现taohe兄现在也在看loren的blog，昨天我还给qibbxxt推荐来着。我准备再过段时间得空的话把她有些文章翻译+改写移过来几篇，比如那个讲cell array的文章，前几个例子写得就相当灵动飘逸。