Some MATLAB tips (二)

 

tips4.关于MATLAB中内联函数,匿名函数,嵌套函数,M-文件函数等的效率问题.

 

从根本上来说,，MATLAB函数分为两类匿名函数（Anonymous Function）和M-文件函数（M-File Function）。下面， 我们从一个例子来说明各类函数地应用与效率比较。

 

问题描述：现在我要画一个数值算法的稳定区域，即要求出在复平面中矩形区域[-2,2;-30;30]内所有能使方程a(u)*x²+b(u)*x+c(u)=0的跟的模小于1的u的值并画图。

 

以匿名函数写的程序如下：

clear all
tic
% 系数；去掉分母后可以认为与原方程同解
a = @(u) exp(u)*(-1+exp(u)-2*u+u*exp(u));
b = @(u) -exp(u)*(-1-2*u+exp(2*u));
c = @(u) exp(2*u)*(-1+exp(u)-u);
% 求根公式
root1 = @(u) -b(u)+sqrt(b(u)^2-4*a(u)*c(u));
root2 = @(u) -b(u)-sqrt(b(u)^2-4*a(u)*c(u));
% 变量初始化
t = 1;
Result=zeros(400*600,1);
% 遍历区域【－2：2；－30：30】；若点满足条件，则存储该值
for ii = -2:0.01:2;
    for jj = -30:0.1:30;       
        if (   abs(root1(ii + i*jj)) < abs(2*a(ii + i*jj))...
            && abs(root2(ii + i*jj)) < abs(2*a(ii + i*jj))  )
            Result(t) = ii + i*jj;
            t = t+1;
        end
    end
end
plot(Result,'.')
toc

 

运行结果：Elapsed time is 86.485000 seconds.

 

 

M－文件函数和嵌套函数程序：

function R1 = root12(u)
% 用求根公式计算方程的两个根
R1 = [ (-b(u)+sqrt(b(u)^2-4*a(u)*c(u)))/(2*a(u)+eps),
       (-b(u)-sqrt(b(u)^2-4*a(u)*c(u)))/(2*a(u)+eps) ]; 

function funa = a(u)   % Subfuntion
% 计算系数a(u)
    funa = exp(u)*(-1+exp(u)-2*u+u*exp(u));
end 

function funb = b(u)   % Subfunction
% 计算系数b(u)
    funb = -exp(u)*(-1-2*u+exp(2*u));
end  

function func = c(u)   % Subfunction
% 计算系数c(u)
    func = exp(2*u)*(-1+exp(u)-u);
end

end

 

clear all
tic
% 变量初始化
t = 1;
Result=zeros(400*600,1);
% 遍历区域【－2：2；－30：30】；若点满足条件，则存储该值
for ii = -2:0.01:2;
    for jj = -30:0.1:30;       
        if (   min(abs(root12(ii + i*jj))<1) ~=0  )
            Result(t) = ii + i*jj;
            t = t+1;
        end
    end
end
plot(Result,'.')
toc

 

运行结果：Elapsed time is 36.328000 seconds.

 

 

不用函数直接把算式写在M脚本里：

clear all;

tic

% 变量初始化
t = 1;
Result=zeros(400*600,1);

% 遍历区域【－2：2；－30：30】；若点满足条件，则存储该值
for ii = -2:0.01:2;
    for jj = -30:0.1:30;       
        %if (   abs(root1(ii + i*jj)) < abs(2*a(ii + i*jj))...
        %    && abs(root2(ii + i*jj)) < abs(2*a(ii + i*jj))  )
        if (    abs( exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj)))...
                + sqrt(((-exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj))))^2)...
                - 4*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj)))*(exp(2*(ii + i*jj))*(-1+exp(ii + i*jj)-(ii + i*jj)))))...
                < abs( 2*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj))))...
            &&  abs( exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj)))...
                - sqrt(((-exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj))))^2)...
                - 4*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj)))*(exp(2*(ii + i*jj))*(-1+exp(ii + i*jj)-(ii + i*jj)))))...
                < abs( 2*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj)))))
            Result(t) = ii + i*jj;
            t = t+1;
        end
    end
end

plot(Result,'.')

toc

 

运行结果：Elapsed time is 26.578000 seconds.

 

结果是很明显，直接将公式写在M文件里比使用函数要快的多。这主要是因为调用函数时系统要做额外的进栈出栈和传值工作。这里没有考虑内联函数（Inline Function），因为他是基于符号计算的，效率比匿名函数要低一些。但是将公式直接写在M文件了，程序的可读性，重用性，可维护性都大大降低。折中的办法也只能是对于调用次数很多的复杂的公式还是要做一个M文件函数的，当然如果掉用次数小公式本身又比较简单，用匿名函数也是自然的事情。

 

 

tips5.关于上述代码的另一点小tip.

 

这一点是无意中发现的。将上述程序做如下小改动（细致点看哦）；  

clear all;

tic

% 变量初始化
t = 1;
%Result=zeros(400*600,1);

% 遍历区域【－2：2；－30：30】；若点满足条件，则存储该值
for ii = -2:0.01:2;
    for jj = -30:0.1:30;       
        %if (   abs(root1(ii + i*jj)) < abs(2*a(ii + i*jj))...
        %    & abs(root2(ii + i*jj)) < abs(2*a(ii + i*jj))  )
        if (    abs( exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj)))...
                + sqrt(((-exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj))))^2)...
                - 4*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj)))*(exp(2*(ii + i*jj))*(-1+exp(ii + i*jj)-(ii + i*jj)))))...
                < abs( 2*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj))))...
            &&  abs( exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj)))...
                - sqrt(((-exp(ii + i*jj)*(-1-2*(ii + i*jj)+exp(2*(ii + i*jj))))^2)...
                - 4*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj)))*(exp(2*(ii + i*jj))*(-1+exp(ii + i*jj)-(ii + i*jj)))))...
                < abs( 2*(exp(ii + i*jj)*(-1+exp(ii + i*jj)-2*(ii + i*jj)+(ii + i*jj)*exp(ii + i*jj)))))
            Result(t) = ii + i*jj;
            t = t+1;
        end
    end
end

plot(Result,'.')

toc

 

看出来没有？其实就是把 %Result=zeros(400*600,1); 关掉了；把 && 改为 &.

 

接下来，看你的 Current Directory 窗口 上方有一个键上面有个小对号√ ：

 

点一下进去回出现下面的界面，ZMY_StableArea_2.m 就是我们上面的脚本的文件名：

 

 

看到了吧！其实这是一个Code Check工具，一些不是太好的做法它会提示修改。例如上面的程序，如果所提示的第31行不做修改的话，运行结果为： Elapsed time is 63.984000 seconds. 这远远比做修改后的26.578000 seconds要慢！这其实容易理解，一个是为数组元素动态分配内存，一种是预分配，显然预分配要快一些，这是典型的空间换时间的做法。