MatLab代码速度和优化。如何提高？

Question

有人可以为我跑这个，告诉我需要多长时间吗？我的笔记本电脑采用了60年代。我无法分辨这是我的笔记本电脑是蹩脚的还是我的代码。可能都是。 我刚开始学习MatLab，所以我还不熟悉哪些功能比其他功能更适合特定任务。如果您对如何改进此代码有任何建议，将不胜感激。

function gbp 
clear; clc; 

zi = 0;         % initial position 
zf = 100;        % final position 

Ei = 1;         % initial electric field 
c = 3*10^8;        % speed of light 
epsilon = 8.86*10^-12;     % permittivity of free space 
lambda = 1064*10^-9;      % wavelength 
k = 2*pi/lambda;       % wave number 
wi = 1.78*10^-3;       % initial waist width (minimum spot size) 
zr = (pi*wi^2)/lambda;     % Rayleigh range 
Ri = zi + zr^2/zi;      % initial radius of curvature 
qi = 1/(1/Ri-1i*lambda/(pi*wi^2));  % initial complex beam parameter 
Psii = atan(real(qi)/imag(qi));   % Gouy phase 
mat = [1 zf; 0 1];      % transformation matrix 

A = mat(1,1); B = mat(1,2); C = mat(2,1); D = mat(2,2); 
qf = (A*qi + B)/(C*qi + D);    % final complex beam parameter 
wf = sqrt(-lambda/pi*(1/imag(1/qf))); % final spot size 
Rf = 1/real(1/qf);      % final radius of curvature 
Psif = atan(real(qf)/imag(qf));   % final Gouy phase 

% Hermite - Gaussian modes function 
u = @(z, x, n, w, R, Psi) (2/pi)^(1/4)*sqrt(exp(1i*(2*n+1)*Psi)/(2^n*factorial(n)*w))*... 
      hermiteH(n,sqrt(2)*x/w).*exp(-x.^2*(1/w^2+1i*k/(2*R))-1i*k*z); 

% Complex amplitude coefficients function 
a = @(n) exp(1i*k*zi)*integral(@(x) Ei.*conj(u(zi, x, n, wi, Ri, Psii)),-2*wi,2*wi); 

%---------------------------------------------------------------------------- 
xlisti = -0.1:1/10000:0.1;    % initial x-axis range 
xlistf = -0.1:1/10000:0.1;    % final x-axis range 
nlist = 0:2:20;       % modes range 

    function Eiplot 
     Efieldi = zeros(size(xlisti)); 
     for nr = nlist 
      Efieldi = Efieldi + a(nr).*u(zi, xlisti, nr, wi, Ri, Psii)*exp(-1i*k*zi); 
     end 
     Ii = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldi); 
    end 

    function Efplot 
     Efieldf = zeros(size(xlistf)); 
     for nr = nlist 
      Efieldf = Efieldf + a(nr).*u(zf, xlistf, nr, wf, Rf, Psif)*exp(-1i*k*zf); 
     end 
     If = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldf); 
    end 

Eiplot 
Efplot 

plot(xlisti,real(Ii),xlistf,real(If)) 

xlabel('x(m)')      % x-axis label 
ylabel('I(W/m^2)')     % y-axis label 
end 

Answer 1

的费用从电话来hermiteH - 为每一个电话，这个符号变量的使用创建了一个新的功能，然后评估在你输入的功能。加快速度的关键是预先计算厄米多项式函数，然后评估这些函数，而不是每次从头开始创建它们（在我的计算机上从〜26秒加速到大约0.75秒）。

有了变化：

function gbp 

x = sym('x'); 

zi = 0;         % initial position 
zf = 100;        % final position 

Ei = 1;         % initial electric field 
c = 3*10^8;        % speed of light 
epsilon = 8.86*10^-12;     % permittivity of free space 
lambda = 1064*10^-9;      % wavelength 
k = 2*pi/lambda;       % wave number 
wi = 1.78*10^-3;       % initial waist width (minimum spot size) 
zr = (pi*wi^2)/lambda;     % Rayleigh range 
Ri = zi + zr^2/zi;      % initial radius of curvature 
qi = 1/(1/Ri-1i*lambda/(pi*wi^2));  % initial complex beam parameter 
Psii = atan(real(qi)/imag(qi));   % Gouy phase 
mat = [1 zf; 0 1];      % transformation matrix 

A = mat(1,1); B = mat(1,2); C = mat(2,1); D = mat(2,2); 
qf = (A*qi + B)/(C*qi + D);    % final complex beam parameter 
wf = sqrt(-lambda/pi*(1/imag(1/qf))); % final spot size 
Rf = 1/real(1/qf);      % final radius of curvature 
Psif = atan(real(qf)/imag(qf));   % final Gouy phase 


% Hermite - Gaussian modes function 
nlist = 0:2:20;  % modes range 

% precompute hermite polynomials for nlist 
hermites = {}; 
for n = nlist 
    if n == 0 
     hermites{n + 1} = @(x)1.0; 
    else 
     hermites{n + 1} = matlabFunction(hermiteH(n, x)); 
    end 
end 

u = @(z, x, n, w, R, Psi) (2/pi)^(1/4)*sqrt(exp(1i*(2*n+1)*Psi)/(2^n*factorial(n)*w))*... 
      hermites{n + 1}(sqrt(2)*x/w).*exp(-x.^2*(1/w^2+1i*k/(2*R))-1i*k*z); 

% Complex amplitude coefficients function 
a = @(n) exp(1i*k*zi)*integral(@(x) Ei.*conj(u(zi, x, n, wi, Ri, Psii)),-2*wi,2*wi); 

%---------------------------------------------------------------------------- 
xlisti = -0.1:1/10000:0.1;    % initial x-axis range 
xlistf = -0.1:1/10000:0.1;    % final x-axis range 

    function Eiplot 
     Efieldi = zeros(size(xlisti)); 
     for nr = nlist 
      Efieldi = Efieldi + a(nr).*u(zi, xlisti, nr, wi, Ri, Psii)*exp(-1i*k*zi); 
     end 
     Ii = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldi); 
    end 


    function Efplot 
     Efieldf = zeros(size(xlistf)); 
     for nr = nlist 
      Efieldf = Efieldf + a(nr).*u(zf, xlistf, nr, wf, Rf, Psif)*exp(-1i*k*zf); 
     end 
     If = 1/2*c*epsilon*arrayfun(@(x)x.*conj(x),Efieldf); 
    end 

Eiplot 
Efplot 

plot(xlisti,real(Ii),xlistf,real(If)) 

xlabel('x(m)')      % x-axis label 
ylabel('I(W/m^2)')     % y-axis label 
end