MATLAB Answers

Translated by

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

0

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

Ahmed S. Sowayan さんによって質問されました 2019 年 2 月 18 日 9:08
最新アクティビティ Star Strider
さんによって コメントされました 2019 年 2 月 20 日 12:19
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);

  15 件のコメント

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

  17 件のコメント

Star Strider
2019 年 2 月 19 日 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')
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
2019 年 2 月 20 日 12:19
My pleasure.

サインイン to comment.



Translated by