Export changing variable with ode45 to solve system of equations in separate function file

I have rather lengthy system of equations saved in a function file 'system3'. I need one of the parameters in this system to change dependant on time. I have created a second function 'calculate_a1' that produces a vector of the parameter a1 at each of my 401 time points.
tResult = [];
xResult = [];
tStep = linspace(0,400,401);
y0 = [IC];
alpha = calculate_a1();
for index = 2:numel(tStep)
% Integrate:
a1 = alpha(1,index);
t = tStep(index-1:index);
sol = ode45(@system3,t,y0,a1)
% Collect the results:
tResult = cat(1, tResult, t);
xResult = cat(1, xResult, x);
% Final value of x is initial value for next step:
y0 = x(end, :);
end
Up until the 'sol' line, this works fine, but I'm struggling to export a1 with ode45 so that it can be used to solve system3.
Any help would be greatly appreciated!

 Respuesta aceptada

Torsten
Torsten el 5 de En. de 2016
Editada: Torsten el 5 de En. de 2016
Forget the time-loop and use interp1 in system3 instead:
tStep = linspace(0,400,401);
alpha = calculate_a1();
y0 = [IC];
tspan = [tStep(1) tStep(end)];
[T,Y]=ode45(@(t,y)system3(t,y,tStep,alpha),tspan,y0);
...
function dy = system3(t,y,tStep,alpha)
alpha_at_t =interp1(tStep,alpha,t);
...
Best wishes
Torsten.

2 comentarios

That seems to run now, but I get the 'Unable to meet integration tolerances without reducing the step size below the smallest value allowed' error. Is this anything to do with the above code, or is that more likely to be an issue with my system3 file?
Torsten
Torsten el 6 de En. de 2016
Editada: Torsten el 6 de En. de 2016
If the (tStep,alpha)-function is smooth, the interpolation should not be a reason for the trouble. You should check the time derivatives "dy" if they contain large values, "NaN" values or "Inf" values.
Best wishes
Torsten.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Mathematics en Centro de ayuda y File Exchange.

Etiquetas

Preguntada:

avg
el 5 de En. de 2016

Editada:

el 6 de En. de 2016

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by