How do you properly vectorize an ODE for integration?

I am interested in ODE integration for a large (10^3 ish) number of initial conditions. As is known, but not terribly well documented, one can do this by "vectorizing" the ODE. I have been doing this with my own hack for awhile, and decided I should finally figure out how to do it properly.

Based on the MATLAB page

http://www.mathworks.com/help/matlab/ref/odeset.html

I thought this would be straightforward.

The following seems to demonstrate that it is not. I would like to know if I am doing something wrong, if this feature is broken in MATLAB, or what.

I should add that I am using MATLAB version 8.0.0.783 (R2012b).

Here are two ODE files, which I have attempted to vectorize to accommodate multiple initial conditions, and a documented script demonstrating that one works and one does not.

Any comments are most appreciated. - John

Here is the one that does not work.

function dydt = vdp1000_problem(t,y)
% This is taken straight from MATLAB documentation page
% http://www.mathworks.com/help/matlab/ref/odeset.html
% It appears to not work in the way described.
dydt = [y(2,:); 1000*(1-y(1,:).^2).*y(2,:)-y(1,:)];

Here is the ODE that does work (although I am not sure it works in the way MATLAB wants it to).

function dydt = vdp1000_works(t,y)
% This function recognizes that even though we pass in an array of column
% vectors, MATLAB reshapes this into a 1D array.

% Dimension of the ODE
dim = 2;
% Number of points we are evolving forward
num_pts = numel(y)/2;
% Reshape in a way that is more convenient to code the ODE
y = reshape(y, [dim, num_pts]);
dydt = [y(2,:); 1000*(1-y(1,:).^2).*y(2,:)-y(1,:)];
% Reshape to match the input (after it has been columnized by ode**).
dydt = reshape(dydt, [dim*num_pts, 1]);

Here is the script.

% This script shows the use of a vectorized ODE.
% Two functions are tested: one, taken from MATLAB doc, does not appear to
% work; the other, by modifying, does appear to work correctly.

% Setup
clear all

vf_problem = @(t, y) vdp1000_problem(t, y);
vf_works = @(t, y) vdp1000_works(t, y);

tspan = linspace(0,0.2,1000);

%%For completeness, test ODE integration with one initial condition.

% yinit = [0.5; 0.3];
% 
% % The "problem" version *does* work when given only one input.
% %[~, xs] = ode45(vf_problem, tspan, yinit);
%
% % This version works too.
% [~, xs] = ode45(vf_works, tspan, yinit);
% 
% figure
% hold on
% plot(xs(:,1), xs(:,2), '.b-')

%%Here is the interesting part. Test ODE integration with many initial conditions
Npts = 10;
yinit = rand([2, Npts]);

% This version does not seem to work at all.
% Tried different integrators (45, 23, ...) with and without
% 'vectorized'='on'.
% opts = odeset('vectorized', 'on');
% [~, xs] = ode15s(vf_problem, tspan, yinit, opts);

% This works for all ode solvers 23, 45, ...
[~, xs] = ode45(vf_works, tspan, yinit);

% % This works for ode 23, 45, 113
% % but not 15s, 23s, 23t, 23tb
% opts = odeset('vectorized', 'on');
% [~, xs] = ode45(vf_works, tspan, yinit, opts);

% The output from the examples that work (with vf_works) are not in the
% same shape as the input.
size(xs)
% The second index goes over coordinates (xA, yA, xB, yB, xC, yC, ...)
% Reshape for easier use.
xs = reshape(xs, [numel(tspan), 2, Npts]);

figure
hold on
cmap = jet(size(xs,3));
for ind = 1:size(xs, 3)
    plot(xs(:,1,ind), xs(:,2,ind), '-', 'color', cmap(ind,:))
end

[Jan]

"Does not work", "seems to demonstrate that it is not", "I am not sure it works in the way MATLAB wants it to": Such terms does not allow to recognize, what you think the problem is. Do you get error messages or do the results differ from your expectations? Why do you think, that your method differs from what Matlab wants? We cannot guess these important details.

Answer 1

[Jan]

While a vectorization should reduce the runtime remarkably, it is a bad idea from a numerical point of view: The initial conditions change the trajectory and in consequence the stepsize control of the integrator. This is required to control the accuracy of the integration. When you run the integration over a set of different initial conditions, the step size is limited by the most sensitive coordinate and therefore the accumulated discretization error is much higher than possible.

[robtheslob]

I am aware that this is an old thread but, as the OP mentioned, information on vectorized ODEs in MATLAB is scarce. I make use of this nearly undocumented feature frequently, so your answer to this question is quite concerning to me.

Could you expand upon your point here? Aren't all of the solutions with differing initial conditions bound in terms of error by AbsTol and RelTol? Isn't it the case that the smallest step size is taken at each iteration for all solutions in order maintain those error bounds for every solution with differing ICs? Is the accumulated discretization error that may result from this procedure not controlled by the aforementioned tolerances?

I understand that you've been sceptical of how error is controlled vis-à-vis ODE integration (I've visited this post of yours more than once) but that was also several years ago now, and you alluded to potential improvements to the code in a 2012 release. Are your points in that post relevant to this discussion?

Thank you for the information.

[Jan]

@robtheslob: If you want to integrate the function

function dy = fcn(t, y)
dy = f(y(:));
end

for different initial conditions, there are two ways:

1. Start the integrator with one start point after the other in a loop.
2. Start the integrator with all initial conditions defined as a vector. Then f(y(:)) must be written such, that it handles vectors correctly in addition.

In both cases, the integrator chooses the size of each step such, that the specified absolute and relative error is not exceeded, but near to this limit to reduce the number of function evaluations. This reduces the processing time and the accumulated rounding errors.

For case 2 the y contains the components for all start points. If then one of the components has large values in higher derivatives, the stepsize is reduced for all components. In consequence, the smoother parts of the trajectory is evaluated too often, which increases the run time and decreases the accuracy or the result.

Think of the trajectory of an asteroid: The initial conditions are a position far away from the sun, and some hundred meters above the earth. Both parts follow the same physical laws, when the density of the air is considered. If you integrate the equations of motion together, the step size will consider the near-earth part with maybe steps of 10e-8 seconds. But for the point far away from the earth, step sizes of hours or weeks might be much better because the trajectory is almost a straight line.

Therefore a vectorization of the integration can counter-productive for the run time and for the accuracy of the output. Different initial conditions or parameters mean different mathematical problems. Then using the same step sizes for the integration is a bad choice.

[robtheslob]

Thank you for the information.

I understand that sets of different initial conditions or different parameter values solved in the "vectorized" manner can lead to the solving of essentially different mathematical problems with the same step size, and this can result in some issues. One issue that you mentioned is that the solver may use unnecessarily small step sizes for problems that do not require it (as in your example of the asteroid) and, as a result, the solver run time increases.

What is less clear to me is the other issue, which is captured in your statement that when "the smoother parts of the trajectory [are] evaluated too often... [this] decreases the accuracy [of] the result." Could you elaborate on what, exactly, the sources of error are when a trajectory is evaluated too often (i.e. the step size is small)? I am aware of round-off error. Also, any information on what can be done to reduce this error (besides taking the largest step size possible that still meets the error tolerances) would be appreciated.

Finally, I have started using book "Solving ODEs with MATLAB" by Shampine et al. as a reference. I would also appreciate any literature recommendations you may have to better understand the theory of numerical ODE solvers. Thanks!

[Jan]

In the evaluation of each step, you get a certain rounding error. During the integration, these errors accumulate: a small deviation in the trajectory tends to increase in the next step.

So smaller step sizes increase the accumulated rounding error, and decrease the local discretization error (the error from the approximation of the derivatives by a quotient of differences). The stepsize controller tries to run the integrator at the optimal stepsize, a compromise between the local discretization error and accumulated rounding error. For this job, several assumptions and approximations are required, and here creating an "optimal controller" becomes to be more art than mathematics. Example: If a step is rejected, the stepsize is usually reduced by the factor 1/2. Many experiments have shown, that this is a good choice. But it will be very hard to find enough examples, which demonstrate clearly, that 1/exp(1) or 1/pi is worse. You will find a lot of papers, which discuss only the value of such parameters.