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

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?

Answer 1

[sam0037]

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

[Jan]

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

Answer 2

[Jan]

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.