1 view (last 30 days)

Show older comments

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

Start Hunting!