MATLAB Answers

Solving Eqn with Varying Variable (Ms)

1 view (last 30 days)
Hello,
I'm tasked to solve for the equation, see attached image, I coded up the equation. I need to find the solution for Ms. In my case I'm solving for M which is Ms. I get a solution, however I need to find multiple solution as i change P4/p1 {1-10000) and a4/a1=[1,2,4,10].
I'm not too good at the loops stuff so I'm very confused as to what loops to use or where to place them. I used a for loop and it just ended up finding the same solution without updating the variables.
I deleted the stuff i felt was wrong and left what give me a solution for one instance
sympref('FloatingPointOutput',true)
gamma1 = 1.667;
gamma4 = 1.4;
syms M
a4a1=1;
a1a4 = 1/a4a1;
eqn = ((2*gamma1*M^2-gamma1+1)/(gamma1+1))*(1-(a1a4)*((gamma4-1)/(gamma1+1))*((M^21)/M))^(-2*gamma4/(gamma4-1)) == 0;
Ms = solve(eqn,M)
this returns
Ms =
-0.4473
0.4473
I need a solution for p4/p1=[10^0-10^4] and a/4/a1=[1, 2, 4, 10] notice the equation is actually asking for a1/a4 which is why i flipped it.

Accepted Answer

Alan Stevens
Alan Stevens on 27 Sep 2021
Like this:
a1a4 = 1./[1, 2, 4, 10];
n = 10000;
p4p1 = 1:n;
M = zeros(numel(a1a4),n);
for j = 1:numel(a1a4)
m = 1.01;
for i = p4p1
M0 = m; % use previous converged value for next initial guess
M(j,i) = fzero(@(M) fn(M,a1a4(j),p4p1(i)),M0);
m = M(j,i);
end
end
plot(p4p1,M),grid
xlabel('p4p1'), ylabel('M')
legend('a4/a1 = 1', 'a4/a1 = 2', 'a4/a1 = 4', 'a4/a1 = 10');
function Z = fn(M,a1a4,p4p1)
gamma1 = 1.667;
gamma4 = 1.4;
t1 = (2*gamma1*M^2-(gamma1-1))/(gamma1 + 1);
t2 = a1a4*(gamma4-1)/(gamma1-1)*(M-1/M);
t3 = -2*gamma4/(gamma4 - 1);
Z = t1*(1 - t2)^t3 - p4p1;
end
  4 Comments
Alan Stevens
Alan Stevens on 28 Sep 2021
Note that p4p1 will depend on the value you choose for a4a1 as well as M. You can rearrange the program as follows:
a1a4 = 1./[1, 2, 4, 10];
n = 10000;
p4p1 = 1:10:n;
M = zeros(numel(a1a4),numel(p4p1));
for j = 1:numel(a1a4)
m = 1.01;
for i = 1:numel(p4p1)
M0 = m; % use previous converged value for next initial guess
M(j,i) = fzero(@(M) fn(M,a1a4(j),p4p1(i)),M0);
m = M(j,i);
end
end
plot(p4p1,M),grid
xlabel('p4p1'), ylabel('M')
legend('a4/a1 = 1', 'a4/a1 = 2', 'a4/a1 = 4', 'a4/a1 = 10');
% Call function fn2 with the desired value of M and a1a4
M_desired = 3; a4a1_desired = 3; a1a4_desired = 1/a4a1_desired;
disp(['M = 3, a4a1 = ' num2str(a4a1_desired) ', p4p1 = ' num2str(fn2(3,a1a4_desired))])
M = 3, a4a1 = 3, p4p1 = 2273.2071
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Z = fn(M,a1a4,p4p1)
Z = fn2(M,a1a4) - p4p1;
end
function pratio = fn2(M,a1a4)
gamma1 = 1.667;
gamma4 = 1.4;
t1 = (2*gamma1*M^2-(gamma1-1))/(gamma1 + 1);
t2 = a1a4*(gamma4-1)/(gamma1-1)*(M-1/M);
t3 = -2*gamma4/(gamma4 - 1);
pratio = t1*(1 - t2)^t3;
end

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!