Solving Eqn with Varying Variable (Ms)

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.

Answer 1

[Alan Stevens]

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

[Benneth Perez]

Hello Alan,

This is great! I'ts exactly what I was looking for. For my own knowledge and if its not too much to ask. Could you explain to me a little bit about the part below. I really want to get better at loops, I usually result in copying and pasting what I want to do multiple times to get the same response.

I under understand that J is just a counter thats basically going to run 4 times.

I'm guessing m=1.01 is the first assumption so that the equation has an m to start with. The part that I'm having the most trouble with is understanding is the commented section, for reference:

Thank you!

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

[Alan Stevens]

The i-loop runs through the values of p4p1 one at a time (so 1, 2, 3,..., 10000).

For each of these, in order to get fzero to find a value of M that makes the function fn equal to zero, we need to give it an initial guess. Sometimes you can get away with a single initial guess for every step in the loop. However, for your equation, it's better to use a guess that comes from the previously converged value. i.e. the initial guess for step i is set to be the converged value of M from step i-1. The very first guess of 1.01 was arbitrary, but, luckily, it worked!

fzero calls the function fn with two parameters in addition to M0. It passes the j'th value of a1a4, and the i'th value of p4p1. fzero then does it's calculations to determine, as best it can, the value of M which gets the value of fn as close to zero as possible.

The returned values of M are then stored in a 4x10000 array (the rows representing the different values of a1a4 and the columns the different values of p4p1).

I hope this helps.

[Benneth Perez]

Thank you I somewhat understand I need more practice and I've looked over your code several times.

If I wanted to find the p4p1 value assuming M=3. In otherwords sort of work backwards. I tried calling the function fn to get an answer but I'm running into some trouble since we were initially solving for all the M values when p4p1 was going from 1-10000.

I went back to my inital simple solve code but its giving me a crazy result

sympref('FloatingPointOutput',true)

gamma1 = 1.667;

gamma4 = 1.667;

M=3;

syms p4p1

a1a4 = 1./[1, 2, 4, 10];

eqn = (((2.*gamma1.*M.^2-gamma1+1)./(gamma1+1)).*(1-(a1a4).*((gamma4-1)./(gamma1+1)).*((M.^21)./M)).^(-2.*gamma4./(gamma4-1)))-p4p1 == 0;

p4p1 = solve(eqn,p4p1)

From the graphs I can see p4p1 should be around 2600 but I want an exact solution for p4p1 when M=3.

Could you stear me in the right direction?

[Alan Stevens]

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