Calculate issue with my code

Can someone help me to fix this issue?

I am a building design student, and I wish use matlab to create a formula then plot the grahpic.

But I can't find the issue with my code

I want to create a formula to calculate the answer like the red box in the table.

But it always happen the error answer after twelve variable.

Can someone help me to fix this issue thank you.

[Walter Roberson]

Which release are you using? In R2019b I do not get that result.

Your output does not appear to match the program. In the output you display a symbolic expression that you do not display in the program.

I notice you have

     D2SBPu(i+1)=min(abs(D2SBPu0(:)));

min(abs()) is a bit risky: it potentially makes complex-valued results equivalent to real-valued results in importance. Also the equation has two roots, the second of which falls faster than the first one; if the second were to fall lower than the first one, then you could end up switching to the second, and then if the second were to fall negative then you could end up switching back to the first one again because the absolute value of the negative could become large.

[Steve Chou]

Dear Roberson

Thank you for your response.

I use R2017a to do this.

Ok, I understand what you mean.

So, I need to fix the min(abs()) problem, but if the formula was correct, how I can code it on matlab to prevent this issue.

The other parts is two root, if I fix the issue with min(abs()), can it possible to fixed the problem with two roots by itself?

Thank you very much.

[darova]

It is a long code. WHat part should be checked?

[Walter Roberson]

I just installed R2017a to check, but the result looks good to me.

Is it possible that you are using the older Symbolic Toolbox, the one that is Maple based?

[Steve Chou]

Dear Roberson

Can you show your result and plot graphic to me?

I want to compare the result different between you and me.

thank you!

[Walter Roberson]

Ah, I did not read the output carefully enough. I will do more debugging.

[Walter Roberson]

I was looking at the wrong equation before. The equation giving problems gets a solve() result as if it were a polynomial of degree 3. The subs() after that has the effect of simplify() . min(abs(BPu0)) of the result appears to have a bug in R2017a. The work-around is min(abs(vpa(BPu0)))

However when I go over to Maple to try to investigate, I see that Maple finds 6 solutions rather than 3 solutions, and back-substitution seems to show all of them as being valid. This is important because the min(abs()) is not going to select the correct result because the solution ranges overlap for the full 6 solutions.

I am investigating to see if there is a better approach.

[Walter Roberson]

In going through the equations, I notice that the three vectors your produce are independent of each other.

The first variable is much easier for the solver to deal with because it does not have abs() around the part that is equated to 1. It appears there are probably three solutions to it, roots of a cubic.

The second variable seems to be a bit harder for solvers to deal with. With the abs() in place, it appears that there are six solutions, only three of which MATLAB can find, corresponding to ignoring either the positive or negative branch that abs() to the same thing. It takes Maple a bit of thought to find the general solution, but not so bad.

The third variable has abs() as well, and for whatever reason is taking Maple rather a long time to deal with. I have reason to suspect that it might have 6 solutions as well.

A question I would have for you is whether all of the named variables involved in the formulae are guaranted to be real-valued? That can eliminate some possibilites. For example it can save me from having to worry about the possibility that the abs() of the expression is having to deal with complex variables, thus permitting me to resolve abs(f(x))=y to f(x)=y union -f(x)=y

[Walter Roberson]

I wrote the following Maple code to try to find the general solution for the variables. The outputs of this code end up being expressed in terms of min() of abs() of two RootOf() structures of cubic polynomials. In Maple you can then extend to exact solutions by using allvalues()

The idea here is to express the forms symbolically, finding the form of the solution once, after which you can substitute in all of your input values simultaneously to get the numeric solutions simultaneously.

But aside from that, you should probably follow the same basic computation strategy in MATLAB: break up the solve(abs(expression)==1) calls into the union of solve(expression==1) and solve(expression==-1) -- assuming the inputs are all real-valued, that is.

fiPn0 := phi_Pn[i + 1];
fi_Mnw0 := phi_Mnw[i + 1];
Pe1w0 := Pe1w[i + 1];
fi_Mnz0 := phi_Mnz[i + 1];
Pe1z0 := Pe1z[i + 1];
Lni[i + 1] := i;

Eq1 := SBPu0/fiPn0 + 8/9*(-SBPu0*Dew/fi_Mnw0 + SBPu0*Dez/(fi_Mnz0*(1 - SBPu0/Pe1z0))) = 1;
sol1 := [solve](Eq1, SBPu0, AllSolutions, explicit = false);
SBPu0_sol := map(simplify, sol1);
SBPu[i + 1] := min(map(abs, SBPu0_sol));

Eq2lhs := BPu0/fiPn0 - BPu0*Dew/(fi_Mnw0*(1 - BPu0/Pe1w0)) + BPu0*Dez/(fi_Mnz0*(1 - BPu0/Pe1z0));
sol2p := [solve](Eq2lhs = 1, BPu0, AllSolutions, explicit = false);
BPu0p := map(simplify, sol2p);
lprint("BPu0 positive done");
sol2n := [solve](Eq2lhs = -1, BPu0, AllSolutions, explicit = false);
BPu0n := map(simplify, sol2n);
lprint("BPu0 negative done");
BPu[i + 1] := min(map(abs, [op(BPu0p), op(BPuon)]));

Eq3lhs := D2SBPu0/fiPn0 + 8/9*(-D2SBPu0*D2ew/fi_Mnw0 + D2SBPu0*D2ez/(fi_Mnz0*(1 - D2SBPu0/Pe1z0)));
sol3p := [solve](Eq3lhs = 1, D2SBPu0, AllSolutions, explicit = false);
D2SBPu0p := map(simplify, sol3p);
lprint("D2SBPu0 positive done");
sol3n := [solve](Eq3lhs = -1, D2SBPu0, AllSolutions, explicit = false);
D2SBPu0n := map(simplify, sol3n);
lprint("D2SBPu0 negative done");
D2SBPu[i + 1] := min(map(abs, [op(D2SBPu0p), op(D2SBPu0n)]));
lprint("D2SBPu0 done");

[Steve Chou]

Dear Roberson

First, thank you for your testing on my code.

I try to use your rewrite code paste to my Matlab but error.

Because I am a new one into Matlab, I only coding with C,C++ before.

So, maybe I am not sure that "Maple code" in my Matlab or not, where I can find it?

Thank you.

[Walter Roberson]

Maple is a third party n Symbolic calculation product https://www.maplesoft.com/

Maple is what the original Symbolic Toolbox was based on, before Mathworks turned to MuPad symbolic engine instead.

The current MuPad symbolic engine has many similarities to Maple, and it probably would not be much work to adapt the code for use in a mupad notebook, using the command

mupad

to open the graphics interface to the MuPad engine. However, that interface will be shut down any release now.

You should be able to rewrite the steps in MATLAB by taking your existing code and splitting the solve(abs(expression) ==1) into two pieces, solve(expression)==1 and solve(expression)==-1 and combine the results.

[Steve Chou]

Dear Roberson

I have already change my code spreat to if Eq2==1 else Eq2==-1

But error, I am stiill not find the problem with my code.

Do I use the wrong coding method, thank you?

for i=0:15
    fiPn0=phi_Pn(i+1);
    fi_Mnw0=phi_Mnw(i+1);
    Pe1w0=Pe1w(i+1);
    fi_Mnz0=phi_Mnz(i+1);
    Pe1z0=Pe1z(i+1);
    
    Lni(i+1)=i;
    
    syms SBPu0
    Eq1=SBPu0/fiPn0+(8/9)*(-SBPu0*Dew/fi_Mnw0+SBPu0*Dez/(fi_Mnz0*(1-SBPu0/Pe1z0)))==1;
    SBPu0=solve(Eq1,SBPu0);
    SBPu0=subs(SBPu0);  SBPu(i+1)=min(abs(SBPu0(:)));
   
        
    syms BPu0
    Eq2=BPu0/fiPn0-BPu0*Dew/(fi_Mnw0*(1-BPu0/Pe1w0))+BPu0*Dez/(fi_Mnz0*(1-BPu0/Pe1z0))==1;
    if Eq2==1
     BPu0=solve(Eq2,BPu0);
    elseif Eq2==-1
     BPu0=solve(Eq2,BPu0);
    end
    BPu0=subs(BPu0); BPu(i+1)=min(abs(BPu0(:)));
    
     syms D2SBPu0
     Eq3=D2SBPu0/fiPn0+(8/9)*(-D2SBPu0*D2ew/fi_Mnw0+D2SBPu0*D2ez/(fi_Mnz0*(1-D2SBPu0/Pe1z0)))==1;
     D2SBPu0=solve(Eq3,D2SBPu0);
     D2SBPu0=subs(D2SBPu0);  D2SBPu(i+1)=min(abs(D2SBPu0(:)));
    
end

[Walter Roberson]

    syms BPu0
    Eq2=BPu0/fiPn0-BPu0*Dew/(fi_Mnw0*(1-BPu0/Pe1w0))+BPu0*Dez/(fi_Mnz0*(1-BPu0/Pe1z0))==1;

That is a symbolic expression containing at least one unresolved variable, BPu0

    if Eq2==1

You cannot use == with an expression with unresolved variables in an if statement. The == creates an equation with an unresolved variable, and if cannot convert that equation into a logical value.

     BPu0=solve(Eq2,BPu0);
    elseif Eq2==-1
     BPu0=solve(Eq2,BPu0);
    end

You do exactly the same thing in both branches, so there is no point having the test.

Furthermore, Eq2 is already an equation, and if it were able to convert to a logical value, then it would convert to symbolic true or symbolic false (not the same as logical true or false), neither of which is equal to 1 or -1, so both branches would fail.

[Steve Chou]

I think matlab is very different with C logic.

So sad with Matlab.

[Walter Roberson]

C defines,

https://en.cppreference.com/w/c/types/boolean

true integer constant 1

false integer constant 0

These are the same values that MATLAB uses. MATLAB also uses exactly the same operator == for comparisons that C uses.

Your Eq2 includes an == operator. In C, if it could be evaluated at all (that is, if BPu0 had a specific numeric value) then the result of the == would be 1 for equality or 0 for inequality, and it is not possible that 1 or 0 can equal -1 for your elseif test.

In C, the only reason to have two tests that executed exactly the same code in either case, would be if the tests had side effects, such as reading from a volatile memory location that was memory mapped to hardware. That would correspond most closely to using a global variable in MATLAB in a situation where a callback might occur between the lines of the tests.

The logic is, in other words, quite similar between MATLAB and C. The logic would not have worked in C and it will not work in MATLAB.

Answer 1

[Walter Roberson]

Like I said, split into two cases and combine the solutions. I even gave a code outline that clearly showed solving the cases individually and combining the results.

for i=0:15
    fiPn0 = phi_Pn(i+1);
    fi_Mnw0 = phi_Mnw(i+1);
    Pe1w0 = Pe1w(i+1);
    fi_Mnz0 = phi_Mnz(i+1);
    Pe1z0 = Pe1z(i+1);
    
    Lni(i+1) = i;
    
    syms SBPu0
    Eq1lhs = SBPu0/fiPn0+(8/9)*(-SBPu0*Dew/fi_Mnw0+SBPu0*Dez/(fi_Mnz0*(1-SBPu0/Pe1z0)));
    SBPu0 = simplify( solve(Eq1lhs == 1, SBPu0) );
    SBPu(i+1) = min(abs(SBPu0(:)));
        
    syms BPu0
    Eq2lhs = BPu0/fiPn0-BPu0*Dew/(fi_Mnw0*(1-BPu0/Pe1w0))+BPu0*Dez/(fi_Mnz0*(1-BPu0/Pe1z0));
    BPu0_plus = simplify( solve(Eq2lhs == 1, BPu0) );
    BPu0_minus = solve(Eq2lhs == -1, BPu0);
    BPu(i+1) = min(abs([BPu0_plus(:); BPu0_minus(:)]));
    
    syms D2SBPu0
    Eq3lhs = D2SBPu0/fiPn0+(8/9)*(-D2SBPu0*D2ew/fi_Mnw0+D2SBPu0*D2ez/(fi_Mnz0*(1-D2SBPu0/Pe1z0)));
    D2SBPu0_plus = simplify(solve(Eq3lhs == 1, D2SBPu0));
    D2SBPu0_minus = simplify(solve(Eq3lhs == -1, D2SBPu0));
    D2SBPu(i+1) = min(abs( [D2SBPu0_plus(:); D2SBPu0_minus(:)] ));
end

[Steve Chou]

Dear Roberson

Thank you for your code.

I paste your code into my matlab but it still get the the same result like before.

I think I really don't know what happen with my code.

[Walter Roberson]

In the code you posted, you have three references to EXPBPu:

%plot(Lni,EXPBPu,'k');

Notice that is commented out

legend('SBPu','BPu','D2SBPu','EXPBPu');

Notice that is as a quoted string and so is not a reference to a variable name. (Notice too that you are giving four legends even though only only plotted three lines)

%double(EXPBPu)

Notice that is commented out.

It is not possible for the code you provided or what I suggested to have "Undefined function or variable 'EXPBPu'"

You have also not given us any equations that we could use to suggest to you what appropriate code might be to calculate EXPBPu .

[Steve Chou]

Dear Roberson

Sorry, my mistake.

I use yo old code and try to fix the problem.

Now, I correct it but still error. It is still happen the same problem with my code.

[Walter Roberson]

Still need to make the adjustment for the bug in R2017a about min() and abs() . I have attached a corrected version.

[Steve Chou]

Dear Roberson

Thank you very much.

I will try it on R2017A, but I am upgrade to R2019A today and it's worked.

It's so glad to meet you.

Answer 2

[Steve Chou]

Ok, I see. Maybe I should be to upgrade my matlab first, thank you very much.