ODE45: IF/loop function
古いコメントを表示
Hello,
I would like to know if its correct to write my differential equations as follow in the 'if' command in the function file. If not, please advice. I am trying to model a simple damper.
Thanks.
function [dy] = SDOF2(t, u);
dy = zeros(2,1);
m=1000; %Mass (lb.sec^2/in.)
k =100000; %Stiffness (lb/in.)
omega = sqrt(k/m); %Natural Frequency
c=2000;%Damping coefficient (lb.sec/in.)
g =386;%Acceleration of gravity (in./sec^2)
c_cr=2*m*omega; %Critical damping coefficient
xi = c/c_cr;%Damping ratio
%%Define the forcing function
% if t<=0.5
% F = sin(4*pi*t);
% else
% F =0;
F = 0;
% b=1*m*(((-omega*omega*u(1)-2*xi*omega*u(2)+F))/u(2)>0);
% end %%%%-ESTIMATION-%%%%
if ((((-omega*omega*u(1))-(2*xi*omega*u(2))+F/m)/u(2))>0)
dy(1) = u(2);
dy(2)= -omega*omega*u(1)-2*xi*omega*u(2)+F/m-0.5*m*((-omega*omega*u(1))-(2*xi*omega*u(2))+F/m);
else
dy(1) = u(2);
dy(2)= -omega*omega*u(1)-2*xi*omega*u(2)+F/m;
end
% dy(2)= -omega*omega*u(1)-2*xi*omega*u(2)+F/m;
end
7 件のコメント
Walter Roberson
2020 年 1 月 9 日
You need to switch to using event functions to terminate the ode45 call and then restart it.
ode45 requires that the function be differentiable. When you have an if in the function it is rarely differentiable. You need to stop at the discontinuity and restart.
Komal Rajana
2020 年 1 月 10 日
編集済み: Komal Rajana
2020 年 1 月 10 日
Walter Roberson
2020 年 1 月 10 日
Strictly speaking, no. However you can subtract eps(realmin) from the value as the result will be negative if the value was exactly 0.
Aquatris
2020 年 1 月 10 日
It might be easier to write your own RK method for this problem.
Komal Rajana
2020 年 1 月 10 日
編集済み: Komal Rajana
2020 年 1 月 10 日
Meg Noah
2020 年 1 月 10 日
There's a solution here on the web:
Also, I've implemented a Runge-Kutta to solve baseball motion under atmosphere drag and lift:
It's a 3-D version.
Aquatris
2020 年 1 月 10 日
There are a lot of examples online, here is on of them.
採用された回答
その他の回答 (0 件)
カテゴリ
ヘルプ センター および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!