solving a system of equations
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
clc
clear all
M1=3
P1=101325 %pa or 1 atm
theta2 = .349 %rad or 20 degrees
theta3 =.262 %rad or 15 degrees
%first get P2 and P3
%area 2
gamma=1.4
lambda2=sqrt((((M1.^2)-1).^2)-3.*(1+((gamma-1)/2).*(M1.^2)).*(1+((gamma+1)/2).*(M1.^2)).*(tan(theta2)).^2);
chi2=(1./(lambda2.^3)).*((((M1.^2)-1).^3)-9.*(1+((gamma-1)./2).*(M1.^2)).*(1+((gamma-1)/2).*(M1.^2)+((gamma+1)./4).*(M1.^4)).*((tan(theta2)).^2));
beta2= atan(((M1.^2)-1+2.*lambda2.*cos((4.*pi+acos(chi2))./3))./(3.*(1+((gamma-1)/2).*(M1.^2)).*tan(theta2)))
M1n2=M1*sin(beta2)
M2n=(((M1n2^2) +(2/(gamma-1)))/((((2*gamma)/(gamma-1))*(M1n2^2)-1)))^(1/2)
M2=M2n/sin(beta2-theta2)
P2=P1*(1+((2*gamma)/(gamma+1))*(M1n2^2 -1))
% repeat for P3
lambda3=sqrt((((M1.^2)-1).^2)-3.*(1+((gamma-1)/2).*(M1.^2)).*(1+((gamma+1)/2).*(M1.^2)).*(tan(theta3)).^2);
chi3=(1./(lambda3.^3)).*((((M1.^2)-1).^3)-9.*(1+((gamma-1)./2).*(M1.^2)).*(1+((gamma-1)/2).*(M1.^2)+((gamma+1)./4).*(M1.^4)).*((tan(theta3)).^2));
beta3= atan(((M1.^2)-1+2.*lambda3.*cos((4.*pi+acos(chi3))./3))./(3.*(1+((gamma-1)/2).*(M1.^2)).*tan(theta3)))
M1n3=M1*sin(beta3)
M3n=(((M1n3^2) +(2/(gamma-1)))/((((2*gamma)/(gamma-1))*(M1n3^2)-1)))^(1/2)
M3=M3n/sin(beta3-theta3)
P3=P1*(1+((2*gamma)/(gamma+1))*(M1n3^2 -1))
%solve for 5 variables 5 unknowns
syms P4 P4p theta4 theta4p beta4 beta4p
P4=P4p
P4=P3*(1+((2*gamma)/(gamma+1))*((M3^2)*((sin(beta4))^2)-1))
P4p=P2*(1+((2*gamma)/(gamma+1))*((M2^2)*((sin(beta4p))^2)-1))
theta4=atan(2*cot(beta4)*((((M3^2)*((sin(beta4))^2)-1))/((M3^2)*(gamma+cos(2*beta4))+2)))
theta4p=atan(2*cot(beta4p)*((((M2^2)*((sin(beta4p))^2)-1))/((M2^2)*(gamma+cos(2*beta4p))+2)))
theta2+theta4p=theta3+theta4
above is my code, i need to solve for theta4, theta4p, beta4, beta4p, and either P4 or P4p, since these two are equal.
not sure which solver would work for this
0 comentarios
Respuestas (1)
Torsten
el 1 de Nov. de 2022
M1=3
P1=101325 %pa or 1 atm
theta2 = .349 %rad or 20 degrees
theta3 =.262 %rad or 15 degrees
%first get P2 and P3
%area 2
gamma=1.4
lambda2=sqrt((((M1.^2)-1).^2)-3.*(1+((gamma-1)/2).*(M1.^2)).*(1+((gamma+1)/2).*(M1.^2)).*(tan(theta2)).^2);
chi2=(1./(lambda2.^3)).*((((M1.^2)-1).^3)-9.*(1+((gamma-1)./2).*(M1.^2)).*(1+((gamma-1)/2).*(M1.^2)+((gamma+1)./4).*(M1.^4)).*((tan(theta2)).^2));
beta2= atan(((M1.^2)-1+2.*lambda2.*cos((4.*pi+acos(chi2))./3))./(3.*(1+((gamma-1)/2).*(M1.^2)).*tan(theta2)))
M1n2=M1*sin(beta2)
M2n=(((M1n2^2) +(2/(gamma-1)))/((((2*gamma)/(gamma-1))*(M1n2^2)-1)))^(1/2)
M2=M2n/sin(beta2-theta2)
P2=P1*(1+((2*gamma)/(gamma+1))*(M1n2^2 -1))
% repeat for P3
lambda3=sqrt((((M1.^2)-1).^2)-3.*(1+((gamma-1)/2).*(M1.^2)).*(1+((gamma+1)/2).*(M1.^2)).*(tan(theta3)).^2);
chi3=(1./(lambda3.^3)).*((((M1.^2)-1).^3)-9.*(1+((gamma-1)./2).*(M1.^2)).*(1+((gamma-1)/2).*(M1.^2)+((gamma+1)./4).*(M1.^4)).*((tan(theta3)).^2));
beta3= atan(((M1.^2)-1+2.*lambda3.*cos((4.*pi+acos(chi3))./3))./(3.*(1+((gamma-1)/2).*(M1.^2)).*tan(theta3)))
M1n3=M1*sin(beta3)
M3n=(((M1n3^2) +(2/(gamma-1)))/((((2*gamma)/(gamma-1))*(M1n3^2)-1)))^(1/2)
M3=M3n/sin(beta3-theta3)
P3=P1*(1+((2*gamma)/(gamma+1))*(M1n3^2 -1))
%solve for 5 variables 5 unknowns
syms P4 P4p theta4 theta4p beta4 beta4p
eqn1 = P4==P4p
eqn2 = P4==P3*(1+((2*gamma)/(gamma+1))*((M3^2)*((sin(beta4))^2)-1))
eqn3 = P4p==P2*(1+((2*gamma)/(gamma+1))*((M2^2)*((sin(beta4p))^2)-1))
eqn4 = theta4==atan(2*cot(beta4)*((((M3^2)*((sin(beta4))^2)-1))/((M3^2)*(gamma+cos(2*beta4))+2)))
eqn5 = theta4p==atan(2*cot(beta4p)*((((M2^2)*((sin(beta4p))^2)-1))/((M2^2)*(gamma+cos(2*beta4p))+2)))
eqn6 = theta2+theta4p==theta3+theta4
S = vpasolve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6],[P4 P4p theta4 theta4p beta4 beta4p])%,[0 inf;0 inf;0 2*pi;0 2*pi;0 2*pi;0 2*pi])
0 comentarios
Ver también
Categorías
Más información sobre Linear Algebra en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!