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

3 次查看（过去 30 天）
Marius
Marius 2016-2-9
回答： Jan 2016-2-15
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?

请先登录，再回答此问题。

回答（2 个）

sam0037
sam0037 2016-2-15
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
Jan 2016-2-15
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.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by