# Solving Eqn with Varying Variable (Ms)

1 view (last 30 days)
Benneth Perez on 27 Sep 2021
Commented: Alan Stevens on 28 Sep 2021
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.

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
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