Скорость и оптимизация кода 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