Having a problem troubleshooting code (relatively simple code)

Running a relatively simple code to plot 2 variables against each other.

clear;clc;

syms A1 B2 SWR

R=0.5;
Cp=0.5;
Cv=0.3;

R = Cv^2*((1/Cv-tan(A1))^2 - (tan(B2))^2)/2/(1-Cv*tan(A1) - Cv*tan(B2)); %%Eq 1

Cp = 1 - (tan(B2)/(1/Cv - tan(A1)))^2; %%Eq2

i = 1;
while Cv<101
      i = i + 1;
      Cv= Cv + 0.01;

      [A1,B2] = vpasolve([R, Cp], [A1,B2]); %%Solving Eq1 and Eq2 simultaneously for A1 and B2

      SWR = 1 - Cv*(tan(A1) + tan(B2)); %%using A1 and B2 to solve for SWR

      X(i,2) = Cv;
      Y(i,2) = SWR;

end

plot(X(:,2),Y(:,2),'g-');
hold on

My end goal here is to plot Cv on the X axis and SWR on the Y axis. But for some reason, it seems to be entering into some sort of infinite loop somewhere since its taking way too long to give an output or even an error. I have tried debugging it but I still am not able to find where the problem is occurring.

Answer 1

[Star Strider]

The Symbolic Math Toolbox is not the best way to solve iterative problems.

I coded your equations as anonymous functions and gave them to the Optimization Toolbox fsolve function that works best for problems like these. It works, but the ‘SWR’ value is negative (and since I believe that stands for ‘standing wave ratio’ on a transmission line, absolutely does not make sense, at least in that context).

Look over your equations to be sure they are set up correctly:

Cp=0.5;
% Cv=0.3;
Cvv = [0.3:0.1:101];

R = @(A1,B2,Cv) Cv.^2.*((1./Cv-tan(A1)).^2 - (tan(B2)).^2)./2./(1-Cv.*tan(A1) - Cv.*tan(B2)); %%Eq 1
Cp = @(A1,B2,Cv) 1 - (tan(B2)./(1./Cv - tan(A1))).^2; %%Eq2

% MAPPING: A1 = p(1),    B2 = p(2)
Eqn = @(p,Cv) [R(p(1),p(2),Cv), Cp(p(1),p(2),Cv)];

for k1 = 1:length(Cvv)
    Cv = Cvv(k1);
    [P(:,k1),Fval] = fsolve(@(p)Eqn(p,Cv), [-145; 23]);
    SWR(k1) = 1 - Cv*(tan(P(1,k1)) + tan(P(2,k1)));
end

figure(1)
plot(Cvv, SWR)
grid
xlabel('C_v')
ylabel('SWR')

I did the symbolic solution first with the initial value of ‘Cv’ (0.3), and used that solution ‘[-145,23]’ for the initial estimates in my loop. If you believe they should be other values, change them, run the loop, and see what works best. Since your equations are periodic, the ‘correct’ initial parameter estimates are important.

[Vidhan Malik]

yeah I'm troubleshooting those equations, something is wrong with the expressions. Also, "SWR" is Stage Work Ratio for compressor blade staging and "Cv" is the simplified flow coefficient! :)

thanks for telling me another way to do these sorts of problems!

[Star Strider]

My pleasure!

There are several numerical solvers that make these sorts of iterative solutions efficient. The Symbolic Math Toolbox is best for initial derivations and simplifications. Its specific strength is in not making stupid algebra errors (that used to occupy a too-large fraction of my time tracking down and correcting). You can use the matlabFunction function to create anonymous functions out of the equations you let the Symbolic Math Toolbox derive for you.

Jet engines aren’t an area of my expertise, so I can’t help with that.

[Vidhan Malik]

I honestly wish I could use the Symbolic Math Toolbox in this case! I'm working with 2 equations (R and Cp) with 4 unknowns, which I have to simplify in terms of just 3, and one of those 3 is Cv and the other two are A1 and B2. So relative to my problem, I have to go through multiple values of Cv and find A1 and B2 for each value of Cv. I don't think the symbolic match toolbox would be able to help me in this case? But I might be wrong! Still not too experienced with matlab just yet!

[Star Strider]

The Symbolic Math Toolbox would likely not be able to help you, except to derive your equations and simplify them, making the final equations more efficient.

I wish I’d thought of this before:

syms A1 B2 Cv

R = Cv^2*((1/Cv-tan(A1))^2 - (tan(B2))^2)/2/(1-Cv*tan(A1) - Cv*tan(B2)); %%Eq 1
Cp = 1 - (tan(B2)/(1/Cv - tan(A1)))^2; %%Eq2

R = simplify(R, 'steps', 10)
Cp = simplify(Cp, 'steps', 10)

R =
(Cv*tan(B2))/2 - (Cv*tan(A1))/2 + 1/2

Cp =
1 - tan(B2)^2/(tan(A1) - 1/Cv)^2

That should make things considerably easier!

[Vidhan Malik]

wow, that is such a neat feature!!! No more toiling over simplifications in the future, THANK YOU!!!

[Star Strider]

As always, my pleasure!

[Vidhan Malik]

Took a break from this problem and came back to it today. I was told by someone who looked over what I was doing and said that I had to initialize the system. Scrolling through what you have said before, that is what you have said/done also at some point but I honestly don't understand why I need to initialize the equations and even then, how it would effect the rest of the iterations. Also, yes the equations were wrong and now am using new ones!

 R= Cv/2*(tan(B1) + tan(B2));
 Cp = 1- (1+tan(B2)^2)/(1+tan(B1)^2);

With the changes to the equations and then modifying my code to be similar to what you had posted before using mapping and fsolve function, I am getting that it never stops displaying an error of :

 Equation solved.

 fsolve completed because the vector of function values is near zero
 as measured by the default value of the function tolerance, and
 the problem appears regular as measured by the gradient.

 <stopping criteria details>

Where the error seems to be happening at the line of

 SWR(k1) = Cv.*(tan(P(1,k1)) - tan(P(2,k1)));

So why exactly is it giving that error? I also updated the initial conditions to what they should be (B1>45 and B2<45 where I solved for a random value of B1 and used it to solve for a value of B2 and make that set the initial condition).

[Vidhan Malik]

Also its really weird but when using the equation

 R= Cv/2*(tan(B1) + tan(B2));

With Cv = 0.4 and B1=60, you are supposed to find that B2 should be around 37.5 but matlab displays a value of -0.3097 so why is matlab getting a different answer?

I used this:

 [B2] = vpasolve([R], [B2]);

[Star Strider]

1. It’s not an error. It just tells you that fsolve did what you wanted it to,
2. the Symbolic Math Toolbox trigonometric functions require their arguments be in radians, not degrees, so

R= Cv/2*(tan(B1) + tan(B2));

should be:

R= Cv/2*(tan(B1*pi/180) + tan(B2*pi/180));

if you want ‘B1’ and ‘B2’ in degrees.

I apologise for the delay on my side — I was away for a while.

[Vidhan Malik]

I tried using

 R= Cv/2*(tan(B1*pi/180) + tan(B2*pi/180));

But matlab is still not reporting the correct value. It just reports B2 to be the negative of B1 which isn't the correct answer.

So for R=0.5, Cv=0.4 and B1=65, B2 should come out to be 19.57 but matlab just displays -65

[Vidhan Malik]

Ah I figured it out! The way I was doing it, I wasn't defining the equation to R at all!

 R= 0.5 == Cv/2*(tan(B1*pi/180) + tan(B2*pi/180));

Using that and solving it gives the right values! Sorry for the trouble! And thank you once again for your amazing help! Now I just had to figure out how to figure the rest out! :)

[Star Strider]

I was just going to suggest that when I read your previous Comment. The fsolve function solves for the roots of its argument functions, such as F(x) = 0. So you have to cast your equations to create that ‘reality’ for fsolve.

No trouble at all! I’m happy got got it sorted.

[Vidhan Malik]

The last problem that I am having is that in the while loop that is setup to solve equations R and Cv, it doesn't seem seem to be increasing Cv for equation R for each iteration, it is keeping it constant for what I had initially defined it.

 R= Cv/2*(tan(B1) + tan(B2));
 Cp = 1- (1+tan(B2)^2)/(1+tan(B1)^2);

So how do I have the solving function inside the while loop increase the value of Cv which appears in the equation that is out of the while loop? I know that you used anonymous functions before I think to address this issue but what if I choose to not do it that way?

The code so far:

Cv=0.39;

R = 0.5 == Cv/2*(tan(B1*pi/180) + tan(B2*pi/180));
Cp = 0.5 == 1- (1+tan(B2*pi/180)^2)/(1+tan(B1*pi/180)^2); 

i = 0;
while i<61
    i = i + 1;
    Cv= Cv + 0.01;

      [B1,B2] = vpasolve([R, Cp], [B1,B2]); 

      SWR = Cv*(tan(B1*pi/180) - tan(B2*pi/180)); 

      X(i,2) = Cv;
      Y(i,2) = SWR;

end

[Star Strider]

With the anonymous functions and with ‘Cv’ as an additional argument, you would not have to include the equations inside the loop. You could just pass the new value of ‘Cv’ to it each time.

Using the Symbolic Math Toolbox, you can either declare ‘R’ and ‘Cp’ as symbolic functions (see the documentation for the symfun function), or you have to put the equations inside the loop so that they will pick up the different values of ‘Cv’ from your workspace.

I would go with the symbolic functions, but either way would be correct.

(However, I would actually go with anonymous functions and fsolve and numerical solutions rather than doing this in the Symbolic Math Toolbox, but the Symbolic Math Toolbox approach will work, if significantly more slowly.)

[Vidhan Malik]

I wish I had the time to learn anonymous functions! But for me currently, who is still learning matlab as a whole, my first instinct is to do it the symbolic math toolbox way since currently it makes more sense to me.

Also when I tried using the symfun function method I get the error:

 Symbolic parameters are not allowed in
 nonpolynomial equations.

So I guess I can't do it that way. Also I tried putting the equations inside the loop, but since we are solving for B1 and B2, we also end up defining it so in the next iteration of the loop, there are no longer any variables in the equation.

[Vidhan Malik]

I think I have to bring the equations inside the loop and then make arrays for Cv, B1 and B2 just like I have with X and Y. I definitely don't want to do it this way but I guess since I am not using anonymous functions, this is my only option.

[Vidhan Malik]

Thank you once again!

[Vidhan Malik]

Ah doing it that way is leading to a whole host of errors haha. So now I am attempting to do it the anonymous function method via the following:

syms B1 B2 SWR Cv n

Cv=0.39;

R = @(Cv) 0.5 == Cv/2*(tan(B1*pi/180) + tan(B2*pi/180)); 

Cp = @(Cv) 0.5 == 1- (1+tan(B2*pi/180)^2)/(1+tan(B1*pi/180)^2); 

i = 0;
while i<60
    i = i + 1;

      Cv= Cv+ 0.01;

      [B1,B2] = fsolve([R(Cv), Cp(Cv)], [B1,B2]); 

      SWR = Cv*(tan(B1*pi/180) - tan(B2*pi/180));

      X(i,2) = Cv;
      Y(i,2) = SWR;

end

This unfortunately gives me the error of:

Error using fsolve (line
148)
FSOLVE requires the
following inputs to be of
data type double: 'X0'.

Error in s (line 16)
    [B1,B2] =
    fsolve([R(Cv), Cp(Cv)],
    [B1,B2]);

Any ideas why?

[Star Strider]

Yes. I use the new equations with my previous fsolve code to get this:

Cvv = [0.3:0.1:101];

R = @(B1,B2,Cv) -0.5 + Cv/2*(tand(B1) + tand(B2)); %%Eq 1
Cp = @(B1,B2,Cv) -0.5 + (1+tand(B2).^2)/(1+tand(B1).^2); %%Eq2

% MAPPING: A1 = p(1),    B2 = p(2)
Eqn = @(p,Cv) [R(p(1),p(2),Cv), Cp(p(1),p(2),Cv)];

for k1 = 1:length(Cvv)
    Cv = Cvv(k1);
    [P(:,k1),Fval] = fsolve(@(p)Eqn(p,Cv), [1; 1]);
    SWR(k1) = Cv*(tand(P(1,k1)) - tand(P(2,k1)));
end

figure(1)
plot(Cvv, SWR)
grid
xlabel('C_v')
ylabel('SWR')

Note that I substituted tand (tangent with degree arguments) for the radian-argument conversion in the code required for the Symbolic Math Toolbox, since you’re using degree arguments. The degree-argument functions are more accurate, especially for small angles.

I use [1; 1] as the initial (degree) estimates for the angles. Change them if that is not correct, since I’ve no idea what they should be.

It takes about 23 seconds to run on my machine.