Matlab - ode 45 y componer un vector con valores de una matriz

   
Vista:

ode 45 y componer un vector con valores de una matriz

Publicado por Ana (1 intervención) el 17/03/2014 11:28:57
Buenas a todos,

Estoy usando ode45 para resolver una ecuación diferencial ordinaria. En este caso lo que quiero es obtener el valor de la Presión (Pr) en función del ángulo theta (theta3). Como ya sabe, ode45 da como resultados dos columnas, que en este caso serán [theta, Pr]

Lo que yo quiero es que los valores numéricos obtenidos en la columna Pr, se me almacenen como componentes de un vector "P" definido anteriormente en el código, y cuyas componentes i=1:280 son conocidas. Quiero que estos valores se almacenen en i=281:340

Además, no sé cómo definir 'theta3' para que se considere como variable en el problema

phi=0.8 ;
rc=15;
lambda=0.25;
b=0.086;
s=0.086;
r=0.043;
R=287;
cv=717;
x2=280;
x3=340;
%-------------------------------------------------------------------

k1= 2 ;
k2= 5000 ;
k3= 14.2./phi.^0.644 ;
k4= 0.79.*k3.^0.25 ;
beta= 0.25 ;

%------------------------------------------------------------------
syms theta2;

fpre=1-(1-((theta2-x2)./b).^k1).^k2 ;
fdif=1-exp(-k3.*((theta2-x2)./b).^k4) ;

Xb= beta.*fpre+(1-beta).*fdif ;
dXb=diff(Xb);
%---------------------------------------------------------------------
DPt=@(theta3,Pr)((dQ-Pr.*dVt).*R./cv-Pr.*dVt)./V ;

V=(pi.*b^2.*s./4).*(1./(rc-1)+1./2.*(1./lambda+1-cosd(theta3)-(1./lambda.^2-sind(theta3.^2).^0.5))) ;
dVt=pi.*b.^2.*r./4.*sind(theta3).*(1+lambda.*cosd(theta3)./(sqrt(1-lambda.^2.*sind(theta3).^2))) ;

dXbsust=subs(dXb, theta2, theta3);
dQ=Qt.*dXbsust ;
tspan=[x2+1 x3];
%P(i)=dsolve('DP(i)=((dQ(i)-P(i).*dVt(i)).*R./Cv-P(i).*dVt(i))./V(i)', 'P(0)=100' , 'theta')
[theta3,Pr]=ode45(DPt,tspan, 100);

for i=281:340

P(i)=[theta3,Pr]

V=(pi.*b^2.*s./4).*(1./(rc-1)+1./2.*(1./lambda+1-cosd(theta(i))-(1./lambda.^2-sind(theta(i).^2).^0.5))) ;
%Derivada de V respecto de theta
dVt=pi.*b.^2.*r./4.*sind(theta(i)).*(1+lambda.*cosd(theta(i))./(sqrt(1-lambda.^2.*sind(theta(i)).^2))) ;
T(i)=P(i).*V(i)./(m.*R) ;
end

figure(1)
title ('Presión')
hold all
grid on
plot (theta,P(i))



Muchas gracias de antemano.
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder