need to solve for a variable numerically

need to solve for M1 based off data and equations.

gamma=1.4;

M=((((Po./P).^((gamma-1)./gamma)-1)*2)/(gamma-1)).^(1/2);
Msonic=find(abs(M-1)<.005)

values=M(305)
values=M(696)
   
Po1=P.*(1+((gamma-1)/2).*M.^2);


index1=1:304; %subsonic
index2=697:1000;%subsonic
index3=305:696; %supersonic

figure
plot(t(index1),M(index1), 'b')
hold on
plot(t(index2),M(index2), 'b')
plot(t(index3),M(index3), 'r')

M2=((((M1).^(2))+(2/(gamma-1)))./((((2*gamma)/(gamma-1)).*(M1).^2)-1))
pr=Po(index3)./P(index3)
preq=((1+((gamma-1)/2).*M2.^2)^(gamma/gamma-1)).*(((gamma+1).*M1.^2)/(2+(gamma-1).*M1.^2))

syms M1 pr gamma

eqn=preq==pr

pretty(eqn)



superM=vpasolve(eqn,M1, 'real')

for index3 i need to solve for M1. pr is calculated using the attached excel document data. Po is stagnation pressure, P is static pressure, t are the indexes.

Answer 1

[David Hill]

m=readmatrix('HW_3_problem_2_data.xlsx');
index3=305:696;
pr=m(index3,2)./m(index3,3);
gamma=1.4;
syms M1
M2=((((M1).^(2))+(2/(gamma-1)))./((((2*gamma)/(gamma-1)).*(M1).^2)-1));
preq=((1+((gamma-1)/2).*M2.^2)^(gamma/gamma-1)).*(((gamma+1).*M1.^2)/(2+(gamma-1).*M1.^2));
for k=1:length(pr)
 superM(k,:)=vpasolve(preq==pr(k),M1);
end