Solving ODEs with different sets of initial conditions
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Amir Kleiner
el 13 de Nov. de 2018
Comentada: Torsten
el 15 de Nov. de 2018
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!
0 comentarios
Respuesta aceptada
Torsten
el 14 de Nov. de 2018
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
2 comentarios
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!