help with ode solvers non-linear ode 1st law of thermodynamics

Hello all:

I am trying to solve the ODE's:

d2y/dt2=-g+(P-Pa)*Ap

dP/dt=k*P*((1/m) *(dm/dt) - (k/(L+y))*(dy/dt))

dm/dt=(m/V) [q-c*Ao*(Pa/P)^(1/k) * {2*k/(k-1) * P*V/m *(1-(Pa/P)^((k-1)/k))}^0.5

V(y)=Ap*(L+y)

where w,Ao,q,Ap,Pa,k,L,g are all constants. but P(t),m(t), y(t) are functions of t

I wrote the following M-file but it looks like it's not working for the m variable. Is there anyting wrong with code or state space variable?

##############################################

clear;close all
global Pa;global A;global Wp;global q;global Ao,global c;global L;global k;
Pa=101000;
d=.05;
Wp=7;
q=2;
Ao=0.01;
c=.1;
L=0.75;
k=1.4;
Te=300;
A=0.25*pi*d^2;
mass=(Pa*L*A)/(287*Te);
tic
W=200;
tf=10;
Fs = .1; %sampling rate
Ts = 1/Fs; %sampling time interval
tspan =0:Fs:tf; %sampling period
y0=[0,.0,101000,mass];
[t,y]=ode45(@Sample01,tspan,y0,[]); 

%y[YfreqDomain_d,frequencyRange_d]=positiveFFT(y(:,1),Fs);

subplot(2,2,1);plot(t,y(:,1),'k');
subplot(2,2,2);plot(t,y(:,2),'k');
subplot(2,2,3);plot(t,y(:,3)/Pa,'k');
subplot(2,2,4);plot(t,y(:,4),'k');

dlmwrite('tssst',[y]);
toc

HERE IS THE FUNCTION OF STATE SPACE

function y_dot=Sample01(t,y)
global Pa;global A;global Wp;global q;global Ao,global c;global L;global k;
y_dot=zeros(4,1);
y_dot(1)=y(2);
y_dot(2)=-Wp*9.81/Wp+(y(3)-Pa)*A/Wp;
y_dot(3)=y(3)*((k/y(4))*y_dot(4)-(k/(L+y(1)))*y(2));

y_dot(4)=y(4)/(L+y(1))*(q-c*Ao*(Pa/y(3))^(1/k)*((2*k/(k-1))*y(3)*(L+y(1))/y(4)*(1-(Pa/y(3))^((k-1)/k)))^0.5);

[Ahmed S. Sowayan]

when I run the code, someting wrong with the mass (m) results are not good..

ANY HELP?

[Star Strider]

I intend to let the Symbolic Math Toolbox re-derive your equations to check them.

This equation makes no sense:

dP/dt=k*P(1/m dm/dt -(k/(L+y)dy/dt)

Please clarify it, adding multiplication operators and parentheses where necessary.

[Ahmed S. Sowayan]

second equation has been corrected now, please check agian.

[Star Strider]

Please also correct .

I prefer to work from your original equations for this.

[Ahmed S. Sowayan]

I think this is exactly my problem dm/dt..

I dont know how to represent this derivative of m in the state variable equations. You may notice that I represented it by y_dot(4) ... is this correct?

[Star Strider]

That is what I am working to figure out. I have no idea what you are doing.

Right now, it is difficult for me to understand your equation because of the substitutions you have made to it in your code, and the lack of operators and parentheses in the way you wrote it symbolically.

I would have to consider your coding your ODEs as state variable equations, although that is certainly possible.

[Ahmed S. Sowayan]

[Ahmed S. Sowayan]

Did you see the original eqaution in the pic above?

[Star Strider]

I do not understand the superscript notation in and just outside the curly braces in the equation. What does that mean?

[Ahmed S. Sowayan]

the curly braces is raised to the power of 0.5

[Ahmed S. Sowayan]

if you run the code results are expected for y, P. But for m it is complex and huge numbers which does make sense.

I think something wrong with state variable representation..

[Ahmed S. Sowayan]

I mean does NOT make sense

[David Goodmanson]

Hi Ahmed,

there is certainly plenty of opportunity for one factor in the dm/dt expression to go funny. If it happens that P becomes less than Pa, then

(1-(Pa/P)^((k-1)/k))^0.5

becomes imaginary.

[Ahmed S. Sowayan]

Hello David

the ratio Pa/P is always of order 1 to 2 but it raised to power less than one therefor the result is less than or equal 1, so it is not must become always imaginary. Also the power 0.5 include more terms.

(2k/(k-1)*P*(V/m)*(1-(Pa/P)^((k-1)/k)))^0.5

Thanks

[David Goodmanson]

Hi Ahmed,

You appear to be saying that Pa/P is greater than one (if not, then this comment does not apply). A number greater than one taken to any positive power is still greater than one, so

PaovP = 1.5  % for example
k = 1.4;
a = (k-1)/k
  a =  0.2857
PaovP^a
  ans =  1.1228
(1-PaovP^a)^(1/2)
  ans = 0.0000 + 0.3505i

>>

Answer 1

[Star Strider]

I let the Symbolic Math Toolbox have its way with your equations:

syms w Ao q Ap Pa k L g t y(t) P(t) m(t) V c Y 
Dy = diff(y);
D2y = diff(Dy);
DP = diff(P);
Dm = diff(m); 
V = L+Dy;
Eq1 = D2y == -g+(P-Pa)*Ap
Eq2 = DP == k*P*((1/m) *(Dm) - (k/V)*(Dy))
Eq3 = Dm == m/V*(q-c*Ao*(Pa/P)^(1/k)*(2*k/(k-1) * (P*V/m)*(1-(Pa/P)^((k-1)/k)))^0.5)
[VF,Sbs] = odeToVectorField(Eq1, Eq2, Eq3)
Sample01 = matlabFunction(VF, 'Vars',{t,Y,k,q,L,Ao,c,Pa,Ap,g})

Sample01 = @(t,Y,k,q,L,Ao,c,Pa,Ap,g)[-(k.^2.*Y(1).*Y(3)-k.*q.*Y(1)+Ao.*c.*k.*Y(1).*(Pa./Y(1)).^(1.0./k).*sqrt((k.*(L+Y(3)).*(Pa-Y(1).*(Pa./Y(1)).^(1.0./k)).*(Pa./Y(1)).^(-1.0./k).*-2.0)./(Y(4).*(k-1.0))))./(L+Y(3));Y(3);-g+Ap.*Y(1)-Ap.*Pa;((q-Ao.*c.*(Pa./Y(1)).^(1.0./k).*sqrt((k.*(L+Y(3)).*(Pa-Y(1).*(Pa./Y(1)).^(1.0./k)).*(Pa./Y(1)).^(-1.0./k).*-2.0)./(Y(4).*(k-1.0)))).*Y(4))./(L+Y(3))];

Pa=101000;
d=.05;
Wp=7;
q=2;
Ao=0.01;
c=.1;
L=0.75;
k=1.4;
Te=300;
A=0.25*pi*d^2;
Ap = 0.42;                                                  % <— Supply Correct Value
g = 9.81;                                                   % <— Supply Correct Value
mass=(Pa*L*A)/(287*Te);
tic
W=200;
tf=10;
Fs = .1; %sampling rate
Ts = 1/Fs; %sampling time interval
tspan =0:Fs:tf; %sampling period
y0=[0,.0,101000,mass];
[t,y]=ode45(@(t,Y)Sample01(t,Y,k,q,L,Ao,c,Pa,Ap,g),tspan,y0); 

subplot(2,2,1);plot(t,y(:,1),'k');
subplot(2,2,2);plot(t,y(:,2),'k');

The problem is that with this code, all the results (except for the initial conditions) are NaN.

Before you run your ode45 code, you will need to comment-out the Symbolic Math Toolbox code to prevent it from seeing the arguments to ‘Sample01’ as symbolic variables, and throwing an error.

Look over the equations I entered here, make any necessary changes, and if you have access to the Symbolic Math Toolbox, run this code yourself with those changes. I am confident that the Symbolic Math Toolbox is correct! However, I am not confident that I coded your equations correctly, and now all except have disappeared,.so I cannot check them.

The ‘Sbs’ output from odeToVectorField are the substitutions it made to create the ‘Y’ vector, and are equivalent to its rows.

[Ahmed S. Sowayan]

can you just post the code without the symbolic stuff. I commented the line you suggested but the answers is NaN..!!!!!!

[Star Strider]

I did. It begins with the anonymous function version of ‘Sample01’, after the matlabFunction call. You can comment-out everything from the syms call to the matlabFunction call, unless you want to re-derive them using corrected functions. I copied the matlabFunction output as the separate ‘Sample01’ anonymous function and pasted it to the code I posted.

I included the Symbolic Math Toolbox code so: (1) you could check the original equations to be sure they are correct, and (2) you would not have to re-type all of it yourself.

The NaN values indicate a ‘0/0’ result, so you need to determine what that problem is. Note that the first 2 of your initial conditions are 0 (apparently corresponding to ‘P(t)’ and ‘y(t)’ in my code), so you might want to be certain those are correct. You may need to re-order your initial conditions vector to reflect the ordering in my code.

Note that:

Sbs =
  P
  y
 Dy
  m

corresponding to ‘Y(1:4)’.

[Ahmed S. Sowayan]

Hi Star Strider;

I run your code but I got an empty figures. The values of y, P, and m are all NaN...

I think the equations that have been generated from the symbolic manipulator are NOT correct.

I made your suggestion by eliminating dy/dm from second equation i.e. I substituted the third equation into the second equation and instead of having sinesoidal graph for y, ydot, and P I got huge numbers which obviously not correct. Here is the code after the update.

[Ahmed S. Sowayan]

function y_dot=Sample01(t,y)

global Pa;global A;global Wp;global q;global Ao,global c;global L;global k;

y_dot=zeros(5,1);

y_dot(1)=y(2);

y_dot(2)=-Wp*9.81/Wp+(y(3)-Pa)*A/Wp;

y_dot(3)=y(3)*(k/(A*(L+y(1)))*(q-c*Ao*(Pa/y(3))^(1/k)*((2*k/(k-1))*y(3)*A*(L+y(1))/y(4)*(1-(Pa/y(3))^((k-1)/k)))^0.5)-(k/(L+y(1)))*y(2));

%y_dot(4)=y(4)/(A*(L+y(1)))*(q-c*Ao*(Pa/y(3))^(1/k)*((2*k/(k-1))*y(3)*A*(L+y(1))/y(4)*(1-(Pa/y(3))^((k-1)/k)))^0.5);

y_dot(4)=y(5);

% incompressible case

%y_dot(4)=y(4)/(L+y(1))*q-c*Ao*((2*y(4)/(L+y(1)))*(y(3)-Pa))^0.5;

[Star Strider]

Check my code to be certain the initial equations (‘Eq1‘, ‘Eq2’, and ‘Eq3’) are correct. If so, run my code with the correct initial conditions. It should work.

And avoid global variables at all costs! Pass the additional parameters as arguments as I did in my code.

‘I run your code but I got an empty figures. The values of y, P, and m are all NaN...’

I know. I mentioned that earlier.

‘I think the equations that have been generated from the symbolic manipulator are NOT correct.’

I believe they are correct, if ‘Eq1‘, ‘Eq2’, and ‘Eq3’ are correct. That’s the reason I asked you to check them to be sure that I entered them correctly, change them if necessary, then run the code again.

Meanwhile, try this (with very small but non-zero initial conditions):

Sample01 = @(t,Y,k,q,L,Ao,c,Pa,Ap,g)[-(k.^2.*Y(1).*Y(3)-k.*q.*Y(1)+Ao.*c.*k.*Y(1).*(Pa./Y(1)).^(1.0./k).*sqrt((k.*(L+Y(3)).*(Pa-Y(1).*(Pa./Y(1)).^(1.0./k)).*(Pa./Y(1)).^(-1.0./k).*-2.0)./(Y(4).*(k-1.0))))./(L+Y(3));Y(3);-g+Ap.*Y(1)-Ap.*Pa;((q-Ao.*c.*(Pa./Y(1)).^(1.0./k).*sqrt((k.*(L+Y(3)).*(Pa-Y(1).*(Pa./Y(1)).^(1.0./k)).*(Pa./Y(1)).^(-1.0./k).*-2.0)./(Y(4).*(k-1.0)))).*Y(4))./(L+Y(3))];

Pa=101000;
d=.05;
Wp=7;
q=2;
Ao=0.01;
c=.1;
L=0.75;
k=1.4;
Te=300;
A=0.25*pi*d^2;
Ap = 0.42;                                                  % <— Supply Correct Value
g = 9.81;                                                   % <— Supply Correct Value
mass=(Pa*L*A)/(287*Te);
tic
W=200;
tf=10;
Fs = .1; %sampling rate
Ts = 1/Fs; %sampling time interval
tspan =0:Fs:tf; %sampling period
y0=[0,.0,101000,mass]+1E-8;
[t,y]=ode45(@(t,Y)Sample01(t,Y,k,q,L,Ao,c,Pa,Ap,g),tspan,y0); 
% [t,y]=ode15s(@(t,Y)Sample01(t,Y,k,q,L,Ao,c,Pa,Ap,g),tspan,y0); 
y = real(y);

figure
subplot(2,2,1);plot(t,y(:,1),'k');
title('P(t)')
subplot(2,2,2);plot(t,y(:,2),'k');
title('y(t)')
subplot(2,2,3);plot(t,y(:,3)/Pa,'k');
title('dy/dt')
subplot(2,2,4);plot(t,y(:,4),'k');
title('m')

[Ahmed S. Sowayan]

Dear Star;

Well I corrected your code according to the correct equations but still I am not confident in the results, at least the huge numbers are gone now, Oscillation is expected in all parameters P, y, ydot, and m. Results for m is increasing linearly with time??? I doubt this result. Any comments?

This is the code after corrections.

[Ahmed S. Sowayan]

clear;

Sample08 = @(t,Y,k,q,L,Ao,c,Pa,Ap,g,Wp)[

Y(1)*(k/(Ap*(L+Y(2)))*(q-c*Ao*(Pa/Y(1))^(1/k)*((2*k/(k-1))*Y(1)*Ap*(L+Y(2))/Y(4)*(1-(Pa/Y(1))^((k-1)/k)))^0.5)-(k/(L+Y(2)))*Y(3));

Y(3);

-g+(Y(1)-Pa)*Ap./Wp;

Y(4)/(Ap*(L+Y(2)))*(q-c*Ao*(Pa/Y(1))^(1/k)*((2*k/(k-1))*Y(1)*Ap*(L+Y(2))/Y(4)*(1-(Pa/Y(1))^((k-1)/k)))^0.5)];

Pa=101000;

d=.05;

Wp=7;

q=20;

Ao=0.01;

c=.9;

L=0.75;

k=1.4;

Te=300;

A=0.25*pi*d^2;

Ap = 0.42; %<— Supply Correct Value

g = 9.81; % <— Supply Correct Value %%% This is correct value of gravitational acceleration

mass=(Pa*L*A)/(287*Te);

tic

W=200;

tf=10;

Fs = .1; %sampling rate

Ts = 1/Fs; %sampling time interval

tspan =0:Fs:tf; %sampling period

y0=[101000,0,0,1]+1E-8;

[t,y]=ode45(@(t,Y)Sample08(t,Y,k,q,L,Ao,c,Pa,Ap,g,Wp),tspan,y0);

% [t,y]=ode15s(@(t,Y)Sample01(t,Y,k,q,L,Ao,c,Pa,Ap,g),tspan,y0);

y = real(y);

figure

subplot(2,2,1);plot(t,y(:,1),'k');

title('P(t)')

subplot(2,2,2);plot(t,y(:,2),'k');

title('y(t)')

subplot(2,2,3);plot(t,y(:,3),'k');

title('dy/dt')

subplot(2,2,4);plot(t,y(:,4),'k');

title('m')

[Ahmed S. Sowayan]

please run the corrected code and look to m ...it is linear, I think, it should not.....wait is Y(4) m or mdot?

[Star Strider]

The code runs without error and appears to give an acceptable result. That’s good enough for me.

I have absolutely no idea what you’re doing, so I have no idea what ‘m(t)’ should be doing.

If you are using my code (that the Symbolic Math Toolbox derived), the ‘Sbs’ output would be interpreted as:

Sbs =
  P
  y
 Dy
  m

so:

P(t)     = Y(1)
y(t)     = Y(2)
dy(t)/dt = Y(3)
m(t)     = Y(4)

[Ahmed S. Sowayan]

Thanks Star;

Your code though gives complex numbers results. If you don't see any errors in modeling the state space variable in my code, then the proper sellection of the parameters q,A,c etc. will overcome this issue.

Thanks anyways

[Star Strider]

My pleasure.

The code are your differential equations that the Symbolic Math Toolbox has composed into first-order differential equations to use with the ODE solvers. You are raising (possibly negative) values to non-integer powers, so that likely generates the complex values. Also, in MATLAB, one complex value in an array results in the entire array being complex, even if all the other imaginary values are zero.

The imaginary parts of the complex numbers are vanishingly small, so when I plotted them, I just use real on the ‘y’ matrix, and that primarily to avoid the warning messages plot throws with complex numbers.

It is relatively easy to optimize your code with respect to the individual parameters (see Optimizing a Simulation or Ordinary Differential Equation (link)), or fit them to observed data if you want to estimate the parameters (see: Fit an Ordinary Differential Equation (link) and Monod kinetics and curve fitting (link)).

[Ahmed S. Sowayan]

The other problem is that y is incresing with time, and it should fluctuate around an equilibrium value same as pressure and ydot. and I think m also should behave the same.

I don't know what to do... :((

[Star Strider]

I’ve no idea what you’re doing (you’ve not described it), much less being certain that I have the background to understand the physics behind it. I don’t recognize the equations, so I’m just helping with the code.

If ‘y’ is increasing and shouldn’t be, consider that you could be integrating a constant that shouldn’t be part of ‘y’. (Maybe one of the initial conditions should be much closer to zero than it is?) That’s the only explanation that comes to mind. It’s also possible that if you find and correct that problem, the other problems will solve themselves.

[Ahmed S. Sowayan]

This is an air spring model in a piston cylinder arrangement. The piston is moving by the effect of gravity and oter forces. The cylinder has two holes on each side of it. One hole is for air inlet (q) by means of fan, and the other is small hole (orifice Ao) for air out of the cylinder.

m_dot=m_dot(in) - m_dot(out)

the balance of forces on the piston is derived using Newton's 2nd law. Equation 1 in the Pic.

so Pressure inside the cylinder is related to the volume change of cylinder and other independent variables (mass of air m, y (distance of piston traveled), y_dot.

This a quik explaination of the problem.

[Star Strider]

Thank you.

I can’t contribute any more than I already have. If your equations are coded correctly in the symbolic code (you need to check that to be certain), the differential equations should do what you expect them to do.

Try this:

syms w Ao q Ap Pa k L g t y(t) P(t) m(t) V c Y 
Dy = diff(y);
D2y = diff(Dy);
DP = diff(P);
Dm = diff(m); 
V = L+y;
Eq1 = D2y == -g+(P-Pa)*Ap/w
Eq2 = DP == k*P*((1/m)*Dm - (1/V)*Dy)
Eq3 = Dm == m/V*(q-c*Ao*(Pa/P)^(1/k)*(2*k/(k-1) * (P*V/m)*(1-(Pa/P)^((k-1)/k)))^0.5)
[VF,Sbs] = odeToVectorField(Eq1, Eq2, Eq3)
Sample01 = matlabFunction(VF, 'Vars',{t,Y,k,q,L,Ao,c,Pa,Ap,g,w})

Pa=101000;
d=.05;
Wp=7;
q=2;
Ao=0.01;
c=.1;
L=0.75;
k=1.4;
Te=300;
A=0.25*pi*d^2;
Ap = 0.42;                                                  % <— Supply Correct Value
g = 9.81;                                                   % <— Supply Correct Value
mass=(Pa*L*A)/(287*Te);
tic
w=200;
tf=10;
Fs = .1; %sampling rate
Ts = 1/Fs; %sampling time interval
tspan =0:Fs:tf; %sampling period
y0=[0,.0,101000,mass]+1E-8;
[t,y]=ode45(@(t,Y)Sample01(t,Y,k,q,L,Ao,c,Pa,Ap,g,w),tspan,y0); 
% [t,y]=ode15s(@(t,Y)Sample01(t,Y,k,q,L,Ao,c,Pa,Ap,g),tspan,y0); 
y = real(y);

figure
subplot(2,2,1);plot(t,y(:,1),'k');
title('P(t)')
subplot(2,2,2);plot(t,y(:,2),'k');
title('y(t)')
subplot(2,2,3);plot(t,y(:,3)/Pa,'k');
title('dy/dt')
subplot(2,2,4);plot(t,y(:,4),'k');
title('m')

[Ahmed S. Sowayan]

Yes thanks, this program is givining the same output as my program.. at least from the symbolic code I know now that my state space equations are correct.

Thanks for your time,,

[Star Strider]

My pleasure.