Asked by Ahmed S. Sowayan
on 18 Feb 2019

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

Answer by Star Strider
on 18 Feb 2019

Accepted 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.

Star Strider
on 19 Feb 2019

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
on 20 Feb 2019

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
on 20 Feb 2019

My pleasure.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 15 Comments

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672496

## Star Strider (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672506

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672516

## Star Strider (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672542

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672564

## Star Strider (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672568

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672570

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672578

## Star Strider (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672581

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672586

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672589

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672590

## David Goodmanson (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672608

## Ahmed S. Sowayan (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672613

## David Goodmanson (view profile)

Direct link to this comment:https://ch.mathworks.com/matlabcentral/answers/445531-help-with-ode-solvers-non-linear-ode-1st-law-of-thermodynamics#comment_672639

Sign in to comment.