Trouble with ODE45

4 vues (au cours des 30 derniers jours)
Christian Kosacki
Christian Kosacki le 13 Mar 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

Réponse acceptée

Jan
Jan le 13 Mar 2017
Modifié(e) : Jan le 13 Mar 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 commentaire
Christian Kosacki
Christian Kosacki le 14 Mar 2017
Okay right, I will look into that. Thanks for the help!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by