MATLAB: Help with ode solvers non-linear ode 1st law of thermodynamics

MATLABode45state-space

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

Best Answer

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.
Related Question