solving a system of equations
1 view (last 30 days)
Show older comments
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 Comments
Answers (1)
Torsten
on 1 Nov 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 Comments
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!