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);endfunction 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!
Best Answer