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!
