### Translated by

このページのコンテンツは英語から自動翻訳されています。自動翻訳をオフにする場合は「<a class="turn_off_mt" href="#">ここ</a>」をクリックしてください。

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

Ahmed S. Sowayan

### Ahmed S. Sowayan (view profile)

さんによって質問されました 2019 年 2 月 18 日

### Star Strider (view profile)

さんによって コメントされました 2019 年 2 月 20 日
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)

2019 年 2 月 18 日
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)

2019 年 2 月 18 日
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)

2019 年 2 月 18 日
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
>>

サインイン to comment.

R2017b

## 1 件の回答

### Star Strider (view profile)

2019 年 2 月 18 日
採用された回答

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)

2019 年 2 月 19 日
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)

2019 年 2 月 20 日
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.