Trouble with ODE45

4 Ansichten (letzte 30 Tage)
Christian Kosacki
Christian Kosacki am 13 Mär. 2017
Kommentiert: Christian Kosacki am 14 Mär. 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

Akzeptierte Antwort

Jan
Jan am 13 Mär. 2017
Bearbeitet: Jan am 13 Mär. 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 Kommentar
Christian Kosacki
Christian Kosacki am 14 Mär. 2017
Okay right, I will look into that. Thanks for the help!

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by