How can I vectorize this code?

Question

George Bashkatov 2021-3-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/770097-how-can-i-vectorize-this-code

评论： George Bashkatov 2021-5-5

i had a loop:

for j=1:length(Pa)
for i=1:(length(r)-1)
    Ip=2*Pa(j)*exp(-2*r(i)^2/wa^2)/(pi*wa^2); % Интенсивность накачки, Вт/м^2
Is=2*Pin*exp(-2*r(i)^2/wa^2)/(pi*wl^2); % Интенсивность сигнала, Вт/м^2
a=Is/Issat;
b=Ip/Ipsat;
zspan=[0 l];
startval=[b; a];
[z1,y1]=ode23(@(z,y) famplifire(sigma_pa,N0,sigma_pe,sigma_se,k,y,z),zspan,startval);
xrx=y1;
itog=itog+y1(end,2)*(r(i+1)^2-r(i)^2)/2;
%tt(i)=y1(end,2)
%rer=sum(tt,'all')
end
qwwer(1,j)=itog*2*pi*Issat;
itog=0;
end

I tryed to transform the first loop on j into something like j=1:1:length(Pa). But I had a problem. Variable startval is a variable that contain initial values for the solver of differential equation. When I change my loop to the vector, no matter either i do this for i or j, I obtain vector b(1,26) or even two vectors a(1,26) and b(1,26). As far as I know, startval should contain only (1,1) vectors or just numbers (also matlab writes that there are wrong dimension of startval if I try to make vectorization). So I need to write another loop specially for a and b. Something like:

for j=1:length(Pa)
startval=[b[j] a];
[z1,y1]=ode23(@(z,y) famplifire(sigma_pa,N0,sigma_pe,sigma_se,k,y,z),zspan,startval);
end

I don't want to do this, because I think that shouldn't make my code work faster.

So, can you told me, can I use vectorization here and how I can do that?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Rik 2021-3-11

在 MATLAB Online 中打开

startval=[b[j] a];

This is not valid Matlab syntax. It is also unclear to me what it attempt to do.

If you want to run ode23 over a grid of values, there might not actually be a shortcut.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2021-3-11

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/770097-how-can-i-vectorize-this-code#answer_645538

编辑：Walter Roberson 2021-4-22

在 MATLAB Online 中打开

Modify famplifire so that it accepts a vector of y values, and reshapes the values to 2 columns, and then computes values for each row. For example, convert

function dy = famplifier(z, y)
  dy = [-y(2); sin(y(1)) + cos(y(2))];
end

into

function dy = famplifier(z, y)
   y = reshape(y, [], 2);
   y1 = y(:,1);
   y2 = y(:,2);
   dy = [-y2, sin(y1) + cos(y2)];
   dy = dy(:);
end

Now create startval as an array with two columns of the boundary conditions

Ip = 2*Pa(:).*exp(-2*r.^2./wa.^2)./(pi.*wa.^2); % Интенсивность накачки, Вт/м^2
Is = 2*Pin.*exp(-2*r.^2./wa.^2./(pi.*wl.^2); % Интенсивность сигнала, Вт/м^2
a = Is/Issat;
b = Ip/Ipsat;
startval = [b(:), a(:)];
[z1,yout]=ode23(@(z,y) famplifire(sigma_pa,N0,sigma_pe,sigma_se,k,y,z),zspan,startval);

yout will then be a 2D array that you can

yout2 = reshape(yout, size(yout,1), [], 2);

First pane, yout2(:,:,1) corresponds to the old y1(:,1), second pain yout2(:,:,2) corresponds to the old y1(:,2), with the rows corresponding to time, the columns corresponding to the different b and a value pairs.

9 个评论
显示 7更早的评论隐藏 7更早的评论

George Bashkatov 2021-4-23

在 MATLAB Online 中打开

1)The old function looks like (before vectorization):

function dydz = famplifire(sigma_pa,N0,sigma_pe,sigma_se,k,y,z)
k=N0*sigma_pa/(sigma_pa+sigma_pe);
dy1 = y(1).*(-sigma_pa.*N0+(sigma_pa+sigma_pe).*(k.*y(1)./(y(1)+y(2)+1))); 
dy2 = y(2).*sigma_se.*(k.*y(1)./(y(1)+y(2)+1)); 
 dydz = [dy1; dy2];

after:

 function dydz = famplifire11(sigma_pa,N0,sigma_pe,sigma_se,k,y,z)
y = reshape(y, [], 2);
k=N0*sigma_pa/(sigma_pa+sigma_pe); 
dy1 = y(:,1).*(-sigma_pa.*N0+(sigma_pa+sigma_pe).*(k.*y(:,1)./(y(:,1)+y(:,2)+1))); 
dy2 = y(:,2).*sigma_se.*(k.*y(:,1)./(y(:,1)+y(:,2)+1)); 
 dydz = reshape([dy1, dy2],[],1);

the function expects 7 parameters.

2) So, the revised code.

clear, close all;
sigma_pa=0.6*10^(-18)*10^(-4);
sigma_pe=0.3*10^(-18)*10^(-4);
sigma_se=0.9*10^(-18)*10^(-4); 
lambda_s=2.5*10^(-6);
lambda_p=1.9*10^(-6);  
c=2.99792458*10^8;
Vp=c/lambda_p; 
Vs=c/lambda_s;
h=6.626*10^(-34);
tau=4*10^(-6); 
wa=80*10^(-6); 
wl=80*10^(-6);
N0=1*10^19*10^6; 
l=7.8*10^-3; 
k=N0*sigma_pa/(sigma_pa+sigma_pe);
Pin=10;
Pa=(0:1:25);
Issat=(h*Vs)/(sigma_se*tau); 
Ipsat=(h*Vp)/((sigma_pa+sigma_pe)*tau);
step=2*wa/10000;
r=0:step:2*wa;
itog=0;
j=1:length(Pa);
for i=1:(length(r)-1)
 Ip = 2*Pa(:).*exp(-2*r.^2./wa.^2)./(pi.*wa.^2); 
Is = repmat(2*Pin.*exp(-2*r.^2./wa.^2./(pi.*wl.^2)), length(Pa), 1); 
a = Is/Issat;
b = Ip/Ipsat;
zspan=[0 1];
startval = [b(:), a(:)];
[z1,yout]=ode23(@(z,y) famplifire11(sigma_pa,N0,sigma_pe,sigma_se,k,y,z),zspan,startval);
yout2 = reshape(yout, size(yout,1), [], 2);
xrx=yout2;
itog=itog+yout2(end,:,2)*(r(i+1)^2-r(i)^2)/2;
end
qwwer(1,j)=itog*2*pi*Issat;
hold on
plot(Pa,qwwer);
xlabel('Мощность накачки, Вт')
ylabel('Выходная мощность генерации, Вт')
grid on