Trouble with ODE45

3 visualizaciones (últimos 30 días)
Christian Kosacki
Christian Kosacki el 13 de Mzo. de 2017
Comentada: Christian Kosacki el 14 de Mzo. de 2017
Hi,
I've got a problem solving my ODE45 function. I have the equation d^2x/dt^2 - sigma^2*x = F(t), meaning the right hand side is a vector. F(t) is an equation containing data from an experiment. To simplify lets say F = 1:1:81.
I have converted the ODE into a first order ODE with x1' = x2 and x2' = F - sigma^2*x1. sigma = 2*pi/t
t = linspace(3.6,7.2,81); %time vector 
F = 1:1:81; % Right hand side
initial_x = 0;
initial_dxdt = 0;
[t,x] = ode45(@funct,t,[initial_x,initial_dxdt]);
plot(t,x(:,1));
xlabel('t'); ylabel('x');
function dxdt = funct(t,x)
dxdt_1 = x(2);
dxdt_2 = F - (2*pi/t)^2*x(1);
dxdt = [dxdt_1;dxdt_2];
end
What I am trying to do is to change the value of F for every value of t, but it seems that I cannot insert any values of add variables into "funct" without getting massive error messages. Do I have to insert ODE45 into a for-loop? Or is there something that I am missing? Thanks in advance

Respuesta aceptada

Jan
Jan el 13 de Mzo. de 2017
Editada: Jan el 13 de Mzo. de 2017
It looks like F is not smooth. Then Matlab's integrators cannot handle this directly, see https://www.mathworks.com/matlabcentral/answers/59582-ode45-errer-input-argument-y-is-undefined#answer_72047.
The best way is a loop:
t = linspace(3.6,7.2,81); %time vector
F = 1:1:81; % Right hand side
x0 = [0, 0];
t0 = t(1);
for iStep = 2:numel(t)
Fs = F(iStep - 1); % For example
[ts, xs] = ode45(@(t, x) funct(t, x, Fs), [t0, t(iStep)], x0);
t0 = ts;
x0 = xs(end, :);
% Store ts and xs for the output...
end
function dxdt = funct(t,x, F)
dxdt = [x(2); ...
F - (2*pi/t)^2*x(1)];
end
This is untested and thought to show the concept only.
  1 comentario
Christian Kosacki
Christian Kosacki el 14 de Mzo. de 2017
Okay right, I will look into that. Thanks for the help!

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by