How can I vectorize this code?

4 views (last 30 days)
George Bashkatov
George Bashkatov on 11 Mar 2021
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 Comment
Rik
Rik on 11 Mar 2021
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.

Sign in to comment.

Answers (1)

Walter Roberson
Walter Roberson on 11 Mar 2021
Edited: Walter Roberson on 22 Apr 2021
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 Comments
George Bashkatov
George Bashkatov on 23 Apr 2021
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
George Bashkatov
George Bashkatov on 5 May 2021
There is no ideas anymore?(

Sign in to comment.

Categories

Find more on Spline Postprocessing in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!