# How to change the ODE to be solved for different ranges in ODE45

17 views (last 30 days)
Marius on 9 Feb 2016
Answered: Jan on 15 Feb 2016
Hi
I have to solve a dynamics problem with MATLAB. I have two differential equations which are only valid for positive resp. negative values of the angle theta. I already transformed them in to the state space, where
y=[theta(t); dtheta(t)] and so dy/dt = dy = [dtheta(t); ddtheta(t)]
with theta(t) the angle, dtheta(t) the velocity and ddtheta(t) the acceleration.
function dy = cylinder_DGL_(t,y,Sys)
% CYLINDER_DGL evaluates the equation of motion for the rocking
% cylinder for given time instants t, excitation ug and vg and system
% properties R and H
% Save geometry properties into new variables
m = Sys.m;
r = Sys.r;
h = Sys.h;
R = Sys.R;
value = Sys.value; %negative or positive for angle alpha, represents the two differential equations about Point O and O'
alpha = Sys.alpha;
g=9.81;
%ODE 1: valid for theta(t) > 0
% Evaluate the equation of motion with help of the state space
% formulation
% State space formulation: row 1
dy(1,1)= y(2);
% State space formulation: row 2
dy(2,1)= -3/4*g/R*sin(alpha-y(1));
%ODE 2: valid for theta(t) < 0
% Evaluate the equation of motion with help of the state space
% formulation
% State space formulation: row 1
dy(1,1)= y(2);
% State space formulation: row 2
dy(2,1)= -3/4*g/R*sin(-alpha-y(1));
end
And then the integration in the script with:
[t,y] = ode45(@(t,y)cylinder_DGL(t,y,Sys),time_r,IC);
What do I have to do now, so that only the valid is used during ODE 45 integration?

sam0037 on 15 Feb 2016
Hi Marius,
In this case both ODE1 and ODE2 can be represented in a single function by a minor change. Replace alpha with aplha*(sign of theta(t)) in the definition of dy(2,1).
Following is the updated code:
function dy = cylinder_DGL_(t,y,Sys)
% CYLINDER_DGL evaluates the equation of motion for the rocking
% cylinder for given time instants t, excitation ug and vg and system
% properties R and H
% Save geometry properties into new variables
m = Sys.m;
r = Sys.r;
h = Sys.h;
R = Sys.R;
value = Sys.value; %negative or positive for angle alpha, represents the two differential equations about Point O and O'
alpha = Sys.alpha;
g=9.81;
% updated code below
k=[];
if(y(1)>=0)
k = alpha;
else
k = -alpha;
end
% State space formulation: row 1
dy(1,1)= y(2);
% State space formulation: row 2
dy(2,1)= -3/4*g/R*sin(k-y(1));
end
(NOTE: I have avoided using a signum function as sign(0)=0)
Thanks,
Shamim

#### 1 Comment

Jan on 15 Feb 2016
Don't do this! Matlab's ODE integrators cannot handle discontinuities reliably. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047

Jan on 15 Feb 2016
You can add an event function, which stops the integration, when the criterion is met. The enclose the integration in a while loop, use the final values of an piecewise integration as initial value for the next integration with a changed paramter.