help with ode solvers non-linear ode 1st law of thermodynamics
8 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Ahmed S. Sowayan
le 18 Fév 2019
Commenté : Star Strider
le 20 Fév 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);
15 commentaires
David Goodmanson
le 18 Fév 2019
Modifié(e) : David Goodmanson
le 18 Fév 2019
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
>>
Réponse acceptée
Star Strider
le 18 Fév 2019
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 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!