ODE45: IF/loop function
Show older comments
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 Comments
Walter Roberson
on 9 Jan 2020
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
on 10 Jan 2020
Edited: Komal Rajana
on 10 Jan 2020
Walter Roberson
on 10 Jan 2020
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
on 10 Jan 2020
It might be easier to write your own RK method for this problem.
Komal Rajana
on 10 Jan 2020
Edited: Komal Rajana
on 10 Jan 2020
Meg Noah
on 10 Jan 2020
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
on 10 Jan 2020
There are a lot of examples online, here is on of them.
Accepted Answer
More Answers (0)
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!