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

Asked by Ahmed S. Sowayan

### Ahmed S. Sowayan (view profile)

on 18 Feb 2019 at 9:08
Latest activity Commented on by Star Strider

### Star Strider (view profile)

on 20 Feb 2019 at 12:19
Accepted Answer by Star Strider

### Star Strider (view profile)

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);

David Goodmanson

### David Goodmanson (view profile)

on 18 Feb 2019 at 20:26
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

### Ahmed S. Sowayan (view profile)

on 18 Feb 2019 at 20:55
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

### David Goodmanson (view profile)

on 18 Feb 2019 at 22:19
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
>>

### Release

R2017b ### Star Strider (view profile)

Answer by Star Strider

### Star Strider (view profile)

on 18 Feb 2019 at 20:43

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.

Star Strider

### Star Strider (view profile)

on 19 Feb 2019 at 21:33
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

### Ahmed S. Sowayan (view profile)

on 20 Feb 2019 at 9:05
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

### Star Strider (view profile)

on 20 Feb 2019 at 12:19
My pleasure.