Trouble with ODE45

3 ビュー (過去 30 日間)
Christian Kosacki
Christian Kosacki 2017 年 3 月 13 日
コメント済み: Christian Kosacki 2017 年 3 月 14 日
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

採用された回答

Jan
Jan 2017 年 3 月 13 日
編集済み: Jan 2017 年 3 月 13 日
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 件のコメント
Christian Kosacki
Christian Kosacki 2017 年 3 月 14 日
Okay right, I will look into that. Thanks for the help!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

製品

Community Treasure Hunt

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

Start Hunting!

Translated by