vpasolve with an array

Hi all,

I was wondering if I could solve for p3 with a given array.

For example,

x1 = 1:0.01:14; 
% or x1 = linspace(1,14,0.01)

syms p3
S = vpasolve(p3*(1 -((gamma-1)*(p3-1))./(2*gamma*(2*gamma+(gamma+1)*(p3-1))).^.5).^(-2.*gamma./(gamma-1))-1 == x1, p3)

% and with given p3 array, I should get Ms = sqrt((gamma+1)/(2.*gamma))*(p3-1)+1, then I need to plot(x1,Ms) so that 
% I could get a continuous curve. 
but, it gives me this error (More equations than variables is only supported for polynomial systems.)

What can I possibly fx for this?  

Answer 1

[Sun Kyoo Choi]

Thank you.

I have successfully solved Ms, but the thing is I have to obtain Mr as well.

I have used the exact same formula that you have provided to me,

x1 = linspace(1,14,1000);

P3G = 3;  %initial guess 
gamma = 1.4; 

S = arrayfun(@(X) fzero(@(p3) p3*(1 -((gamma-1)*(p3-1))./sqrt(2*gamma*(2*gamma+(gamma+1)*(p3-1)))).^(-2.*gamma./(gamma-1)) - X, P3G), x1);

Ms = (((gamma+1)/(2.*gamma))*(S-1)+1).^0.5;

% So far, it is good. I gain S and Ms Matrix. 
% But, I need to plug in Ms array to gain Mr.

S1 = arrayfun(@(X) fzero(@(Mr) ((Ms.*(Mr.^2-1))/(Ms.^2-1))*((1+((2.*(gamma-1)*(Ms.^2-1)))./(gamma+1)^2).*(gamma+1./Ms.^2)).^.5 - Mr - X, P3G), Ms);
...
    it gives me "Operands to the || and && operators must be convertible to logical scalar values".

%I am trying to find Mr arrays with respect to Ms. 

[Walter Roberson]

Ms is a vector so you are trying to fzero with an expression that is a vector.

[Sun Kyoo Choi]

So, what exactly should I fix?

[Walter Roberson]

Notice that the right hand side does not involve M_R. You are defining a quadratic equation and you do not need to approximate solutions for that.

right = Ms./(Ms.^2 - 1) .* sqrt(1 + 2.*(gamma-1)./(gamma+1)^2.*(Ms.^2-1).*(gamma+1./Ms.^2));

You now have M_R/(M_R^2-1) = right which is the same as M_R = right * (M_R^2-1) which is 0 = right*M_r^2 - M_R - right which is clearly a quadratic and can be solved with the usual closed form for quadartics, yielding two answers for each Ms value.

[Sun Kyoo Choi]

0 = right*M_R^2 - M_R - right

% Since M_R is unknown matrix or array, I was wondering if I had to use
% Mr_theo = [right, -1, -right] 
% Mr = roots(Mr_theo)
it gives me NaN or Inf... 

[Walter Roberson]

I think your right will be a vector.

Check, though:

disp(length(right))
disp(nnz(isfinite(right)))

[Sun Kyoo Choi]

        1000

   999
   
   This was my results for
   disp(length(right))
disp(nnz(isfinite(right)))

What does this represent? 

If right is a vector, how am I supposed to write code for Mr? 

[Walter Roberson]

One of your right values is nan or infinite . That could happen if Ms is infinite or 1, and possibly under other circumstances.

When you accidentally put all of the right values into a single vector to find the roots of, you were putting that non-finite value into the same vector, and roots() of that extended vector would have been nan or infinite.

arrayfun(@(V) roots([V, -1, V]), right, 'uniform', 0)

will give you a cell array of the roots, one cell for each Ms value.

Or you could use the standard (-b±sqrt(b^2-4ac))/2a formula to calculate the roots in vectorized form. b=-1, a=right, c=-right

[Sun Kyoo Choi]

Thank you!

I really appriciate it! :)

Answer 2

[Walter Roberson]

Replace

syms p3

With

p3 = sym('p3_', [1, length(x1)]) ;

This will make p3 into an array of different symbols, p3_1, p3_2 and so on, and will solve for each of those.

This will not necessarily be efficient, but it does satisfy your requirements that vpasolve solve for an array of equation.

Remember that vpasolve attempts to find values for the variables that solve all of the equations simultaneously, so your code was trying to find a single p3 that worked for the entire x1 vector.

My suspicion is that it would be more efficient to arrayfun a vpasolve or fzero or fsolve call.

Answer 3

[Sun Kyoo Choi]

Then, what would be a possible way to fix this problem, because when I replace p3 to what you have provided, it gives me this error. How could I use arrayfun for this?

x1 = linspace(1,14,1000);

gamma = 1.4;

p3 = sym('p3_', [1, length(x1)]) ;

S = vpasolve(p3*(1 -((gamma-1)*(p3-1))./(2*gamma*(2*gamma+(gamma+1)*(p3-1))).^.5).^(-2.*gamma./(gamma-1))-1 == x1, p3);

Error using symengine

Dimensions do not match.

[Walter Roberson]

Remember to change the * to .* and the / to ./

[Sun Kyoo Choi]

Seems like it is working, but it is taking forever to run. It has been taking over 10mins to run this.

It is because it has 1x1000 matrix? (huge)

I doubt it though.

[Walter Roberson]

My suspicion is that it would be more efficient to arrayfun a vpasolve or fzero or fsolve call.

[Sun Kyoo Choi]

I was trying somewhat like this, but it still gave me an error.

Since I have never used arrayfun, I do not really know the format.

In my case what would I have to modify to use arrayfun?

Thank you for keep answering my quesiton.

S = arrayfun(@vpasolve,eqn,'uniform',0);

[Walter Roberson]

P3G = 2;  %initial guess 
S = arrayfun(@(X) fzero(@(p3) p3*(1 -((gamma-1)*(p3-1))./(2*gamma*(2*gamma+(gamma+1)*(p3-1))).^.5).^(-2.*gamma./(gamma-1))-1 - X, P3G), x1);

[Sun Kyoo Choi]

why did yoiu put P3G all of sudden?

Is it not possble to create an array of p3 without P3G?

[Walter Roberson]

What is your goal here? Is it to find P3 to more than 16 decimal places accuracy such that vpasolve() is needed for extended precision? Is it to find P3 to double precision but you must have the last couple of bits as accurate as possible so that is why you are using vpasolve to get the extended digits that you can then truncate with double()? Is it to write the simplest program possible and you find vpasolve to be easier to read than the numeric alternatives? Is it to solve to reasonable double precision accuracy in which case double() is fine without using vpasolve? Is performance an issue?

My arrayfun solution was created around the most common need: reasonable double precision accuracy and reasonable performance combined. It uses the numeric solver fzero to find P3. fzero happens to need an initial guess.

If you feel strongly about the matter I could recode it in terms of the higher precision but slower vpasolve, which does not require an initial guess (but usually is faster if you supply a guess.) If you do not pass an initial guess to vpasolve it generates its own initial guess, probably 0 or 1.

[Sun Kyoo Choi]

Okay, that was correct. Sorry my mistake.

My final question is that when I plot x1 vs. Ms, where Ms is defined

%Ms = (((gamma+1)/(2.*gamma))*(S-1)+1).^0.5;
%S is p3 an array that I have created. 
P3G = 3;  %initial guess 
S = arrayfun(@(X) fzero(@(p3) p3*(1 -((gamma-1)*(p3-1))./(2*gamma*(2*gamma+(gamma+1)*(p3-1))).^.5).^(-2.*gamma./(gamma-1))-1 - X, P3G), x1);
Ms = (((gamma+1)/(2.*gamma))*(S-1)+1).^0.5;
plot(x1,Ms);

%It give me this shape of graph, which is correct, but is there any way that I can
%start from (1,1)? How could I adjust it? Currently, it starts from (1,1.1595).
%I would like to fix so that it could start from (1,1) instead of (1,1.1595). 

[Walter Roberson]

In order for MS to be 1, S would have to be 1. In order for S to be 1, p3 would have to be 1. You want MS to be 1 when x1 is 1, so you can substitute X = 1 and p3 = 1 into

p3*(1 -((gamma-1)*(p3-1))./(2*gamma*(2*gamma+(gamma+1)*(p3-1))).^.5).^(-2.*gamma./(gamma-1))-1 - X

The answer falls out to -1 no matter what the gamma value is. -1 is never equal to 0, not even for very small values of -1.

Therefore, No, it would be inconsistent for MS to be 1 at x1 = 1, no matter what the gamma value is.