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)
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)
dim = 2;
num_pts = numel(y)/2;
y = reshape(y, [dim, num_pts]);
dydt = [y(2,:); 1000*(1-y(1,:).^2).*y(2,:)-y(1,:)];
dydt = reshape(dydt, [dim*num_pts, 1]);
Here is the script.
clear all
vf_problem = @(t, y) vdp1000_problem(t, y);
vf_works = @(t, y) vdp1000_works(t, y);
tspan = linspace(0,0.2,1000);
Npts = 10;
yinit = rand([2, Npts]);
[~, xs] = ode45(vf_works, tspan, yinit);
size(xs)
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
1 Comment
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/110356-how-do-you-properly-vectorize-an-ode-for-integration#comment_186494
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/110356-how-do-you-properly-vectorize-an-ode-for-integration#comment_186494
Sign in to comment.