How can I vectorize this code?

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?

[Rik]

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]

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.

[George Bashkatov]

Matlab writes that error in line

startval = [b(:), a(:)];

Dimensions of arrays being concatenated are not consistent.

In workspace I see, that a is a 1x10001 double and b is a 26x10001 double

[Walter Roberson]

Ah, I see why that would happen, but it will take me some thinking on how best to fix it.

[Walter Roberson]

Ip = 2*Pa(:).*exp(-2*r.^2./wa.^2)./(pi.*wa.^2); % Интенсивность накачки, Вт/м^2
Is = repmat(2*Pin.*exp(-2*r.^2./wa.^2./(pi.*wl.^2)), length(Pa), 1); % Интенсивность сигнала, Вт/м^2
a = Is/Issat;
b = Ip/Ipsat;
startval = [b(:), a(:)];

[George Bashkatov]

Thank you for your help, but the code doesn't work((

1) You forget to write

zspan=[0 l];

before

[z1,yout]=ode23(@(z,y) famplifire(sigma_pa,N0,sigma_pe,sigma_se,k,y,z),zspan,startval);

2) We vectorize the startval.My first code gives one value during each iterration. So, after the code to be executed it gives me two values: [0.0029;0.0015]. After vectorization my code gives me an array 260026x2. And the last values in each column should coincide with [0.0029;0.0015]. I see, that the value in the first column coincide with 0.0029, but the second is not. And I don't know how to fix that.

Because of this problem, my programs works infinitely without any result. Maybe the code has another problems too, but I can't analyse them by reason of startval.

[Walter Roberson]

I cannot test as you did not post your input values.

[George Bashkatov]

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;

After all calculations I need to obtain qwwer (in the first code variant you can see how it can be calculated) and plot it by the command:

plot(Pa,qwwer);

[Walter Roberson]

You have

[z1,yout]=ode23(@(z,y) famplifire(sigma_pa,N0,sigma_pe,sigma_se,k,y,z),zspan,startval);

which passes 7 parameters to famplifire but that function only expects two parameters ?

Also please post your revised code to calculate qweer now that you are no longer looping the ode23 call.

[George Bashkatov]

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]

There is no ideas anymore?(