Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
ODE stiffness problems

Subject: ODE stiffness problems

From: Matthew

Date: 13 Mar, 2013 17:20:09

Message: 1 of 6

Hi, I am working on a project simulating a moving hydraulic cylinder with fluid being pushed out and drawn in through different ends. The fluid moves through valves which are opening and shutting.

I didn't know anything about functions and ODEs until recently, when someone tried to help me solve the problem. However they do not know what to do now.

I'm getting this error 'Warning: Failure at t=9.649647e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t. '

There are 4 integrals which need to be calculated by the ODE, position of the piston head, velocity of the piston head, pressure inside side A of the chamber and pressure inside chamber B of the cylinder.

I am following an example that I found, and I get the correct answers up until 1 second. When I simulate more than 1 second I get the error. I have tried messing with the tolerances and using all the different ODE solvers. Can someone please help me. Thanks.

The function is listed below followed by the main script.

function [dz, Av, Q] = cylinderoderhs(t,z,design)

% Calculating the size of the valve open area as a function of time
if t<design.t1
    Av = design.w * design.xvmax * (t/design.t1);
elseif t<design.t2
    Av = design.w * design.xvmax;
elseif t<=design.t3
    Av = design.w * design.xvmax * ((1-t)/design.t1);
else
    Av = 0;
end

% Working out the flow rate through either side of the valve.
Q(1,1) = sign(design.po-z(1)) * design.cd * Av * ((2*abs(design.po-z(1))/design.rho)^(1/2));

Q(1,2) = sign(design.ps-z(2)) * design.cd * Av * ((2*abs(design.ps-z(2))/design.rho)^(1/2));

dz = zeros(4,1);

% Calculating the change in pressure in chamber side 1, dp1/dt
dz(1) = (design.b / (design.Vt - design.Acap * z(3))) * (Q(1) + design.Acap * z(4));

% Calculating the change in pressure in chamber side 2, dp2/dt
dz(2) = (design.b /(design.Vt + design.Acap * z(3))) * (Q(2) - design.Acap * z(4));

% Calculating the position of the cylinder, xa. Which is the integrated
% va.
dz(3) = z(4);

% calculating the velocity of the cylinder, va.
dz(4) = (z(2) * design.Acap - z(1)*design.Arod - design.c*z(4) - design.ma*design.g) / design.ma;

end




% The times for the simulation start and end.
span = [0 1];

% System parameters
design.ma = 750; % Mass
design.da = 0.05; % Diameter of the piston
design.drod = 0.025; % Diameter of the rod
design.V = 15.7e-6; % Volume from the valve to the cylinder (in pipework)
design.la = 0.5; % Stroke length of the cylinder
design.lsl = 0.01;
design.lgap = 20e-6; % Anulus gap
design.w = 0.01; % Valve width opening
design.xvmax = 0.5e-3; % Maximum valve opening (Vertical)
design.cd = 0.62; % Orafice discharge coefficient
design.g = 9.81; % Gravity
design.rho = 855; % Density
design.b = 1.2e9; % Bulk modulus
design.mu = 15e-3; % Absolute viscosity
design.ps = 20e6; %System pressure
design.po = 0; %Drain pressure

% Annulus areas are the same.
design.Acap = (design.da^2 * pi/4) - (design.drod^2*pi/4);
design.Arod = design.Acap;

% The total volume in each side of the cylinder at x=0 (Center)
design.Vt = design.V + (design.Acap * design.la)/2;

% The damping coefficent.
design.c = pi * design.da * design.lsl * design.mu / design.lgap;

% Starting postion, xa, velocity, va and pressures either side of the
% chamber.
sxa = 0;
sva = 0;
sp2 = 750*9.81/(0.05^2*pi/4-(0.025^2*pi/4));
sp1 = 0;

% Times that the valve opens. From 0-t1 the valve is opening, it is open
% fully until t2 is reached then takes 0.05 seconds to close.
design.t1 = 0.05;
design.t2 = 0.95;
design.t3 = 1;

% Altering the relative tolerance of the ODE from default 1e-6.
%odeoptions = odeset('RelTol', 1e-3);
odeoptions = [];

% solve the ode function.
[t, z] = ode15s(@cylinderoderhs, span, [sp1 sp2 sxa sva], odeoptions, design);

% Extracting the internally calculated values
for i =1:numel(t)
    [dz(i,:), Av(i,1), Q(i,1:2)] = cylinderoderhs(t(i), z(i,:), design);
end

Subject: ODE stiffness problems

From: Torsten

Date: 14 Mar, 2013 08:52:06

Message: 2 of 6

"Matthew" wrote in message <khqcg9$anl$1@newscl01ah.mathworks.com>...
> Hi, I am working on a project simulating a moving hydraulic cylinder with fluid being pushed out and drawn in through different ends. The fluid moves through valves which are opening and shutting.
>
> I didn't know anything about functions and ODEs until recently, when someone tried to help me solve the problem. However they do not know what to do now.
>
> I'm getting this error 'Warning: Failure at t=9.649647e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t. '
>
> There are 4 integrals which need to be calculated by the ODE, position of the piston head, velocity of the piston head, pressure inside side A of the chamber and pressure inside chamber B of the cylinder.
>
> I am following an example that I found, and I get the correct answers up until 1 second. When I simulate more than 1 second I get the error. I have tried messing with the tolerances and using all the different ODE solvers. Can someone please help me. Thanks.
>
> The function is listed below followed by the main script.
>
> function [dz, Av, Q] = cylinderoderhs(t,z,design)
>
> % Calculating the size of the valve open area as a function of time
> if t<design.t1
> Av = design.w * design.xvmax * (t/design.t1);
> elseif t<design.t2
> Av = design.w * design.xvmax;
> elseif t<=design.t3
> Av = design.w * design.xvmax * ((1-t)/design.t1);
> else
> Av = 0;
> end
>
> % Working out the flow rate through either side of the valve.
> Q(1,1) = sign(design.po-z(1)) * design.cd * Av * ((2*abs(design.po-z(1))/design.rho)^(1/2));
>
> Q(1,2) = sign(design.ps-z(2)) * design.cd * Av * ((2*abs(design.ps-z(2))/design.rho)^(1/2));
>
> dz = zeros(4,1);
>
> % Calculating the change in pressure in chamber side 1, dp1/dt
> dz(1) = (design.b / (design.Vt - design.Acap * z(3))) * (Q(1) + design.Acap * z(4));
>
> % Calculating the change in pressure in chamber side 2, dp2/dt
> dz(2) = (design.b /(design.Vt + design.Acap * z(3))) * (Q(2) - design.Acap * z(4));
>
> % Calculating the position of the cylinder, xa. Which is the integrated
> % va.
> dz(3) = z(4);
>
> % calculating the velocity of the cylinder, va.
> dz(4) = (z(2) * design.Acap - z(1)*design.Arod - design.c*z(4) - design.ma*design.g) / design.ma;
>
> end
>
>
>
>
> % The times for the simulation start and end.
> span = [0 1];
>
> % System parameters
> design.ma = 750; % Mass
> design.da = 0.05; % Diameter of the piston
> design.drod = 0.025; % Diameter of the rod
> design.V = 15.7e-6; % Volume from the valve to the cylinder (in pipework)
> design.la = 0.5; % Stroke length of the cylinder
> design.lsl = 0.01;
> design.lgap = 20e-6; % Anulus gap
> design.w = 0.01; % Valve width opening
> design.xvmax = 0.5e-3; % Maximum valve opening (Vertical)
> design.cd = 0.62; % Orafice discharge coefficient
> design.g = 9.81; % Gravity
> design.rho = 855; % Density
> design.b = 1.2e9; % Bulk modulus
> design.mu = 15e-3; % Absolute viscosity
> design.ps = 20e6; %System pressure
> design.po = 0; %Drain pressure
>
> % Annulus areas are the same.
> design.Acap = (design.da^2 * pi/4) - (design.drod^2*pi/4);
> design.Arod = design.Acap;
>
> % The total volume in each side of the cylinder at x=0 (Center)
> design.Vt = design.V + (design.Acap * design.la)/2;
>
> % The damping coefficent.
> design.c = pi * design.da * design.lsl * design.mu / design.lgap;
>
> % Starting postion, xa, velocity, va and pressures either side of the
> % chamber.
> sxa = 0;
> sva = 0;
> sp2 = 750*9.81/(0.05^2*pi/4-(0.025^2*pi/4));
> sp1 = 0;
>
> % Times that the valve opens. From 0-t1 the valve is opening, it is open
> % fully until t2 is reached then takes 0.05 seconds to close.
> design.t1 = 0.05;
> design.t2 = 0.95;
> design.t3 = 1;
>
> % Altering the relative tolerance of the ODE from default 1e-6.
> %odeoptions = odeset('RelTol', 1e-3);
> odeoptions = [];
>
> % solve the ode function.
> [t, z] = ode15s(@cylinderoderhs, span, [sp1 sp2 sxa sva], odeoptions, design);
>
> % Extracting the internally calculated values
> for i =1:numel(t)
> [dz(i,:), Av(i,1), Q(i,1:2)] = cylinderoderhs(t(i), z(i,:), design);
> end

1. Call the integrator as
[t, z] = ode15s(@(t,z)cylinderoderhs(t,z,design), span, [sp1 sp2 sxa sva], odeoptions);
2. Define your function as
 function dz = cylinderoderhs(t,z,design)
3.Insert Q(1,1) and Q(1,2) in the calculation of dz(1) and dz(2) instead of Q(1) and
Q(2).
4. For clarity, change your if-statement in the calculation of Av to
 if t<design.t1
     Av = design.w * design.xvmax * (t/design.t1);
 elseif t>= design.t1 & t<design.t2
     Av = design.w * design.xvmax;
 elseif t>= t.design.t2 & t<design.t3
     Av = design.w * design.xvmax * ((1-t)/design.t1);
 else
     Av = 0;
 end
5.The sign and abs-functions may introduce problems because of their non-differentiability.

Good luck!

Best wishes
Torsten.

Subject: ODE stiffness problems

From: Matthew

Date: 14 Mar, 2013 17:49:17

Message: 3 of 6

> 1. Call the integrator as
> [t, z] = ode15s(@(t,z)cylinderoderhs(t,z,design), span, [sp1 sp2 sxa sva], odeoptions);
> 2. Define your function as
> function dz = cylinderoderhs(t,z,design)
> 3.Insert Q(1,1) and Q(1,2) in the calculation of dz(1) and dz(2) instead of Q(1) and
> Q(2).
> 4. For clarity, change your if-statement in the calculation of Av to
> if t<design.t1
> Av = design.w * design.xvmax * (t/design.t1);
> elseif t>= design.t1 & t<design.t2
> Av = design.w * design.xvmax;
> elseif t>= t.design.t2 & t<design.t3
> Av = design.w * design.xvmax * ((1-t)/design.t1);
> else
> Av = 0;
> end
> 5.The sign and abs-functions may introduce problems because of their non-differentiability.
>
> Good luck!
>
> Best wishes
> Torsten.

I made the changes you mentioned however, it still does not work. Is there anyway to get past the sign and abs? It is required to get the equation working. I might be able to get past it by introducing if statements.

The goal of this is create a hydraulic system for a wave energy converter and I am using this example as a basis as it is very similar. However, the wave energy converter will require quite a few more integrals in the ODE.

Subject: ODE stiffness problems

From: Matthew

Date: 14 Mar, 2013 23:11:06

Message: 4 of 6

I removed the abs and signs by putting a constant low flow rate just to see if it would work. However it doesn't make any change it still stops at the same point.

Subject: ODE stiffness problems

From: Torsten

Date: 15 Mar, 2013 08:54:07

Message: 5 of 6

"Matthew" wrote in message <khtlea$ba3$1@newscl01ah.mathworks.com>...
> I removed the abs and signs by putting a constant low flow rate just to see if it would work. However it doesn't make any change it still stops at the same point.


I solved your problem with another software and it worked without problems.

The problem with ODE15s seems to be the abrupt change of the variables
in the time span from t=0.95 to t=1. This is caused by the parameter Av.

A remedy (I also used within my ODE integrator) is to call ODE15s seperately for
the time span from t=0 to t=0.05, then from t=0.05 to t=0.95 and then from t=0.95 to t=1.
Do it in a loop and set the inital conditions of the variables for a new time span to the final values of the variables from the preceeding time span.

Best wishes
Torsten.

Subject: ODE stiffness problems

From: Marc

Date: 23 Apr, 2013 01:49:13

Message: 6 of 6

"Torsten" wrote in message <khunje$be2$1@newscl01ah.mathworks.com>...
> "Matthew" wrote in message <khtlea$ba3$1@newscl01ah.mathworks.com>...
> > I removed the abs and signs by putting a constant low flow rate just to see if it would work. However it doesn't make any change it still stops at the same point.
>
>
> I solved your problem with another software and it worked without problems.
>
> The problem with ODE15s seems to be the abrupt change of the variables
> in the time span from t=0.95 to t=1. This is caused by the parameter Av.
>
> A remedy (I also used within my ODE integrator) is to call ODE15s seperately for
> the time span from t=0 to t=0.05, then from t=0.05 to t=0.95 and then from t=0.95 to t=1.
> Do it in a loop and set the inital conditions of the variables for a new time span to the final values of the variables from the preceeding time span.
>
> Best wishes
> Torsten.

Torsten,

Is this "other" integrator you speak of Limex? If so, I found setting ODE15s MaxOrder to 2 sometimes help get similar/almost identical results. For my system, ODE23t also got me similar results.

I also used a similar "trick" of breaking up the time span manually.

Another set of integrators I have had nice luck with is the sundials ode integrators, code and ida (for DAEs) which are free and have a nice mex setup in Matlab.

Not 100% sure why setting the MaxOrder to 2 worked but it did.

My 2 cents....

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us