solving a system of equations

1 view (last 30 days)
joshua payne
joshua payne on 1 Nov 2022
Answered: Torsten on 1 Nov 2022
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

Answers (1)

Torsten
Torsten on 1 Nov 2022
M1=3
M1 = 3
P1=101325 %pa or 1 atm
P1 = 101325
theta2 = .349 %rad or 20 degrees
theta2 = 0.3490
theta3 =.262 %rad or 15 degrees
theta3 = 0.2620
%first get P2 and P3
%area 2
gamma=1.4
gamma = 1.4000
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)))
beta2 = 0.6590
M1n2=M1*sin(beta2)
M1n2 = 1.8370
M2n=(((M1n2^2) +(2/(gamma-1)))/((((2*gamma)/(gamma-1))*(M1n2^2)-1)))^(1/2)
M2n = 0.6084
M2=M2n/sin(beta2-theta2)
M2 = 1.9943
P2=P1*(1+((2*gamma)/(gamma+1))*(M1n2^2 -1))
P2 = 3.8204e+05
% 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)))
beta3 = 0.5629
M1n3=M1*sin(beta3)
M1n3 = 1.6009
M3n=(((M1n3^2) +(2/(gamma-1)))/((((2*gamma)/(gamma-1))*(M1n3^2)-1)))^(1/2)
M3n = 0.6682
M3=M3n/sin(beta3-theta3)
M3 = 2.2543
P3=P1*(1+((2*gamma)/(gamma+1))*(M1n3^2 -1))
P3 = 2.8609e+05
%solve for 5 variables 5 unknowns
syms P4 P4p theta4 theta4p beta4 beta4p
eqn1 = P4==P4p
eqn1 = 
eqn2 = P4==P3*(1+((2*gamma)/(gamma+1))*((M3^2)*((sin(beta4))^2)-1))
eqn2 = 
eqn3 = P4p==P2*(1+((2*gamma)/(gamma+1))*((M2^2)*((sin(beta4p))^2)-1))
eqn3 = 
eqn4 = theta4==atan(2*cot(beta4)*((((M3^2)*((sin(beta4))^2)-1))/((M3^2)*(gamma+cos(2*beta4))+2)))
eqn4 = 
eqn5 = theta4p==atan(2*cot(beta4p)*((((M2^2)*((sin(beta4p))^2)-1))/((M2^2)*(gamma+cos(2*beta4p))+2)))
eqn5 = 
eqn6 = theta2+theta4p==theta3+theta4
eqn6 = 
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])
S = struct with fields:
P4: 286008.97317717601854893168488697 P4p: 286008.97317717601854893168488697 theta4: 0.00008460289275915508457608111229808 theta4p: -0.086915397107240844915423918887702 beta4: -163.82235899238345755662732994989 beta4p: 9.8849808344785758940601584211874

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!