Solving vectorized ODE (Solving same ODE with many initial conditions at once).

9 views (last 30 days)
I am tring to apply 2nd newton law to many points, where m,r,v,f is a n-by-1, n-by-3, n-by-3, n-by-3, n-by-3, matrixes for mass, position,velocity, and force for n points. The 3 columns are the x,y,z components for those values.
options = odeset('Vectorized','on');
y0=[r,v];
[t,YSol]=ode45(@(t,y) MotionODE(t,y,m,f),[0,dt],y0,options);
r=YSol(t==dt,1:3)';
v=YSol(t==dt,4:5)';
function d2rdt2= MotionODE(t,Y,m,f)
%ODE function for motion
r=Y(:,1:3);
v=Y(:,4:5);
drdt=v;
dvdt=f./m;
d2rdt2=[drdt,dvdt];
end
However I got an error: "Index in position 2 exceeds array bounds (must not exceed 1).". Is the option "vectorized" designed for what I think it is for? How do I get this run correctly (other than a for loop)?
  1 Comment
Yi-xiao Liu
Yi-xiao Liu on 21 Mar 2021
Walter's answer was accepted because it answers my original question. But Jan's answer is also informative regarding accuracy of this approach and I encourge everyone to give a read.
I also found this link to be revelant:

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 17 Mar 2021
r = Y(1,:);
v = Y(2,:);
drdt = v;
dvdt = f./m * ones(size(v));
d2rdt2=[drdt;dvdt];
  3 Comments
Walter Roberson
Walter Roberson on 17 Mar 2021
Okay, I misread earlier. 'Vectorized' is not useful for your situation.
%m is n x 1
%r is n x 3
%v is n x 3
%f is n x 3
%the boundary conditions are arranged in memory as
%all of the position x coordinates, then all of the position y, then all
%of the position z, then all of the velocity x, then all of the velocity y,
%then all of the velocity z
y0 = [r,v]; %[n x 3, n x 3] --> n x 6
f_over_m = f ./ m; %n x 3 ./ n x 1 -> n x 3
[t,YSol] = ode45(@(t,y) MotionODE(t, y, f_over_m), [0,dt], y0);
R = reshape(YSol(end,1:end/2), [], 3); %px, py, pz
V = reshape(YSol(end,end/2+1:end), [], 3); %vx, vy, vz
function d2rdt2= MotionODE(t, Y, f_over_m)
%ODE function for motion
Y = reshape(Y, [], 6); %px, py, pz, vx, vy, vz
r = Y(:,1:3); %n x 3
v = Y(:,4:5); %n x 3
drdt = v; %n x 3
dvdt = f_over_m; %n x 3
d2rdt2 = reshape([drdt, dvdt], [], 1); %MUST return column vector
end
Yi-xiao Liu
Yi-xiao Liu on 17 Mar 2021
it should be v=Y(:,4:6) in the ODE function, which is the mistake I originally made in the question.
Other than that it works great, Thanks!

Sign in to comment.

More Answers (1)

Jan
Jan on 17 Mar 2021
same ODE with many initial conditions at once
This is a bad idea. Remember that the step size control is triggered by the most susceptible component of the trajectory. This reduces the stepsize for all components and in consequence increases the accumulated rounding errors. The total number of calculations can be larger than running the integration for each initial value separately.
  4 Comments
Yi-xiao Liu
Yi-xiao Liu on 18 Mar 2021
I guess my question is more like how to solve this ODE with different initial conditions independently, as you suggested for better accuracy, while also concurrently?
Using a for loop takes forever and most of the CPU is idel at the time
Jan
Jan on 18 Mar 2021
Do you have the Parallel Processing Toolbox? Then try to replace the FOR by PARFOR.

Sign in to comment.

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!