How can I vectorize this code?

4 次查看（过去 30 天）
George Bashkatov
George Bashkatov 2021-3-11
评论： 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 个评论
Rik
Rik 2021-3-11
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.

请先登录，再进行评论。

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

回答（1 个）

Walter Roberson
Walter Roberson 2021-3-11
编辑：Walter Roberson 2021-4-22
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更早的评论
George Bashkatov
George Bashkatov 2021-4-23
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spline Postprocessing 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by