Solving ODEs with different sets of initial conditions

I have a system of ODEs which I solve using ode45.

In my case, I am interested in solving the same equation hundreds and even thousands of times, each time with slightly different randomly generated initial conditons, as follows (the sovled function is just an example):

n_times = 1000;
t_i = 1; t_f = 10; tspan = [t_i,t_f];
for i = 1:n_times
    y0 = randn(2, 1);
    [t, y] = ode45(@(t,y) harmosc(t,y), tspan, y0);
end

function dy = harmosc(t,y)
    dy = zeros(size(y));
    dy(1)=y(2);
    dy(2)=-y(1);
end

In my current code I run a loop for all the initial conditions, but this is rather slow for the large amount of time I have to solve the equations. I have tried solving this entering the initial conditions as a matrix instead of looping over a vector, something like:

n_times = 1000;
t_i = 1; t_f = 10; tspan = [t_i,t_f];
y0 = randn(2, n_times);
[t, y] = ode45(@(t,y) harmosc(t,y), tspan, y0);


function dy = harmosc(t,y)
    dy = zeros(size(y));
    dy(1,:)=y(2,:);
    dy(2,:)=-y(1,:);
end

But this doesn't work - I manage to get a y vector of the wrong size ([length(t), 2*n_times]), and with only one solution (y(1,:), y(2,:)) while all other columns of the matrix are just some constant.

Does someone know how this can be done correctly?

I am aware of the fact that the loop results in better accuracy of the solution, but that's a price I'm willing to pay for the reduced time in the "matrix" way.

Thanks!

Answer 1

[Torsten]

function main
  n_times = 1000;
  t_i = 1; t_f = 10; tspan = [t_i,t_f];
  y0 = randn(2, n_times);
  y0 = reshape(y0, 2*n_times, 1);
  [t, y] = ode45(@(t,y) harmosc(t,y,n_times), tspan, y0);
  y = reshape(y,size(y,1),2,n_times)
  plot(t,y(:,1,15),t,y(:,2,15))
end

function dy = harmosc(t,y,n_times)
    y = reshape(y,2,n_times);
    dy = zeros(size(y));
    dy(1,:)=y(2,:);
    dy(2,:)=-y(1,:);  
    dy = reshape(dy,2*n_times,1);
end

[Amir Kleiner]

Thanks for the answer!

Why is it like this? why do I need to reshape the y0 vector?

[Torsten]

Because the ODE integrators work with vectors, not with matrices.