# How do you properly vectorize an ODE for integration?

10 views (last 30 days)
John Mahoney on 21 Dec 2013
Commented: Jan on 5 Dec 2017
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
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 on 21 Dec 2013
"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.

Jan on 21 Dec 2013
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.
Jan on 5 Dec 2017
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.