22 views (last 30 days)

Show older comments

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

Jan
on 21 Dec 2013

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.

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

Start Hunting!