How can I vectorize this code?

3 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
zspan=[0 l];
startval=[b; a];
[z1,y1]=ode23(@(z,y) famplifire(sigma_pa,N0,sigma_pe,sigma_se,k,y,z),zspan,startval);
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);
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 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))];
function dy = famplifier(z, y)
y = reshape(y, [], 2);
y1 = y(:,1);
y2 = y(:,2);
dy = [-y2, sin(y1) + cos(y2)];
dy = dy(:);
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
George Bashkatov on 5 May 2021
There is no ideas anymore?(

Sign in to comment.




Community Treasure Hunt

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

Start Hunting!