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

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)?

[Yi-xiao Liu]

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:

https://www.mathworks.com/help/matlab/math/solve-system-of-odes-with-multiple-initial-conditions.html

Answer 1

[Walter Roberson]

r = Y(1,:);
v = Y(2,:);
drdt = v;
dvdt = f./m * ones(size(v));

d2rdt2=[drdt;dvdt];

[Yi-xiao Liu]

Error using vertcat

Dimensions of arrays being concatenated are not consistent.

BTW what is it doing? I placed a breakpoint in the MotionODE and Y is passed in as a 6n-by-1 vector. so

r = Y(1,:);

v = Y(2,:);

only assigned them as a scalar.

[Walter Roberson]

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]

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!

Answer 2

[Jan]

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.

[Yi-xiao Liu]

I found the error to be on the magnitude of ~1e-10 for a n=500 using Walter's solution above, which was acceptable I guess.

Do you have any suggestion to improve performance?

[Jan]

The error can be small, if the trajectory is not very sensitive for a variation of the initial values. If you integrate a chaotic function like a double pendulum, the error will get much larger.

You can improve the performance by choosing another integrator or my modifying the tolerances. In your case the function to be uintegrated is tiny. Therefore saving function calls is not the standard way to improve the speed. If this function contains thousands of lines of code, the situation would be different.

[Yi-xiao Liu]

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]

Do you have the Parallel Processing Toolbox? Then try to replace the FOR by PARFOR.