I have problems solving a cubic equation

Hello,

I have a cubic equation that should yield 2 real and 1 complex solutions.

(Based on the book, from where the equation originates.)

%Original equation with k1 and k2 being constants and R the variable of interest.
syms R V
k1 = 0.3105*V^2;
k2 = 0.0594*V^2;
F = R^6 - k1*R^2 + k2 == 0;

%Substituting R^2 with r
syms r
f = r^3 - k1*r + k2 == 0;

%Solving 
r = simplify(solve(f,r,'maxdegree',3));
R = sqrt(r)

%Plot
figure
plot([0.2:0.05:1],subs(R(1),[0.2:0.05:1]))
hold on
plot([0.2:0.05:1],subs(R(2),[0.2:0.05:1]))

However, when I solve this I get only 1 real and 2 complex solutions.

Maybe my math is off, but the real part of my results are in right order of magnitude...

The Radius curve in the picture below is the desired result. The upper curve would be solution 1 and the lower curve solution 2, or vice versa, respectively. The reactor power is proportional to V...

I use Version R2018a.

Any help is highly appreciated! :D

[Star Strider]

Getting two real and one imaginary roots from a cubic equation does not seem correct to me. (I would expect two complex-conjugate roots and one real root.)

What procedures does the book use to get those results? What are the results in the book?

[Felix Richter]

"This equation gives us three solutions for the square of the square of R (r), but one is not real."

One real and two complex solutions are what I get.

The results are depicted in the figure above. The Radius curve is the result I am looking for. I do not know what you mean with procedures, they just substituted a quadratic variable to get a cubic equation.

[Walter Roberson]

you do not have a cubic you have a sextic, six roots, with closed form solutions available.

[Star Strider]

By ‘procedures’, I am asking what the book does to go from a 6-degree polynomial to a 3-degree polynomial?

[Felix Richter]

Yes initially a sextic. I am not really sure if substitution changes this, or if this just a work around.

Do you have a suggestion what to do about it?

I also tried:

simplify(solve(F,R,'maxdegree',6));

But that doesn't yields the solution I was looking for either.

[Felix Richter]

@Star Strider: Substituting R^2=r. That's it.

[Star Strider]

That does not seem correct to me. With that sort of transformation, I would assume that ‘k1’ and ‘k2’ also change, however they do not.

Something just seems wrong to me about the book’s analysis.

[Felix Richter]

Why would they, k1 and k2 are indepedent from R. And substitution should'nt influence any other than R or r, respectively, shouldn't it?

The book is generally sientifically very correct and I don't think an error like this is probable.

[John D'Errico]

At a very quick glance, what seems to be wrong in the analysis may be just a typo. Just based on units, the way V appears in each of k1 and k2 seem to be a bit strange. Unless V and R are completely unit-less parameters, this would make no physical sense:

k1 = 0.3105*V^2;
k2 = 0.0594*V^2;
F = R^6 - k1*R^2 + k2 == 0;

So my initial conjecture is the book either contains a typo, or Felix has committed a copy-o. Of course, I cannot know if that is true, since we do not see the book, nor do we know the book from whence this was taken.

Looking more carefully at the picture, I see a mention of radius, and I also see a height. Each of those values are apparently measured in meters. Radius is likely to be the variable R. And in that case, R will have units applied to it. Again, if R is measured in some sort of units, then that equation will make no physical sense, unless V also has some units attached. And if V DOES have units attached, then k1 and k2 cannot both have V^2 in them. The units just won't work.

That is, suppose R has units of meters. Then R^6 has units of m^6. To that term you now subtract k1*R^2. So I must now assume that V has units of meters^2, which will give k1 units of meters^4, else you cannot meaningfully add those terms. But then we see k2 must also have units of meters^4. And this makes no sense at all, unless everything is non-dimensional.

(If you need to ask why, then you need to consider what those equations would mean if everything were measured in cm or even feet, instead of m.)

So there is some problem in your equation. Again, this is only a conjecture. Perhaps at some point that we are not shown, the authors have non-dimensionalized the problem. Perhaps R is a variable that has been scaled so it is relative to the height. (This would make R a unit-less variable. And as long as V also has no units attached, then all would make some sense.)

The most probable resolution to this issue would be if we really had something like:

k2 = 0.0594*V^3;

Now all the units would make some sense, as long as V is indeed some sort of area parameter, which would give it units of length squared.

Again, I lack your book to discern the truth here, so my "answer' is only a conjecture, though a fairly strong one. But the formula as written is mathematically meaningless if R and V have units attached to them.

If nothing of what I have said makes sense, you might wish to do some reading here:

https://en.wikipedia.org/wiki/Dimensional_analysis

Read carefully the section about dimensional homogeneity, where you will find the statement:

"The most basic rule of dimensional analysis is that of dimensional homogeneity.[6]1. Only commensurable quantities (physical quantities having the same dimension) may be compared, equated, added, or subtracted."

The concept of non-dimensional parameters can be found in the Buckinham pi theorem.

https://en.wikipedia.org/wiki/Buckingham_π_theorem

And you don't want to know how many years it has been since I was in a class, with my professor discussing the Buckinham pi theorem. I won't tell anyway.

[Star Strider]

John — My sentiments exactly, in your Comment described in appropriate detail.

[Felix Richter]

Hello John, first of all thanks for your extensive analysis of the problem!

In my initial question I "dumbed it down" and on the way the dimension got lost. I am really sorry for leading you on the wrong track...

This is the original formula:

(R_core^2)^3 - (B*V_core/pi^2)^2 * (R_core^2) + (2.405*V_core/pi^2)^2 == 0

%where

k1 = (B/pi^2)^2 * V^2 = 0.3283*V^2
%I changed one parameter a bit, compared to the initial question

k2 = (2.405/pi^2)^2 *V^2 = 0.0594*V^2

%R_core [m]
%V_core [m^3]
%B [1/m]

So, dimensionwise there should definitely be no error.

However, I just re-ran the m.file and I am getting now, for instance the following solution.

 %For 
 P_core = 1.5e9; %[W]

 %it yields
 
 0.63939 - 1.1732e-10i
 0.46169 + 2.7665e-10i
 4.9293e-11 - 0.78866i

The real parts of the first solutions are just that what I was looking.

When I plot it for a core power range of 1.2e9 to 2e9 (W), it got the following result:

That is pretty close to the books result.

Can I just assume the imaginary parts are negligible and go on with my life? :D

I really don't know WHY there are imaginary parts in the solution, but I am studying engineering and not math, so that might be okay...

[David Goodmanson]

Hello Felix,

The first two roots have a small enough imaginary part that you are probably all right ignoring it. Although in these kind of calculations, an imaginary part down around 1e-13 or less is more common and would inspire more confidence that it's just numerical error. Possibly the fact that P_core is so large could be causing scaling problems and reducing the accuracy.

So moving on seems all right, but even though your task is focused on engineering, it doesn't hurt to know more about the roots. Taking the (mostly) lower curve R_core, it looks like a pitchfork bifurcation as shown in wiki, although the middle tine to the right is missing. There are good reasons why the book may not show that one.

How this result depends on P_core is something you have not posted, so the plot above can't be reproduced. There are going to be six roots at all times. Did you just list three of them in the P_core = 1.5e9 example? Given the plot, that pure imaginary root does not seem consistent with the first two real roots.

Note that the legend in the plot above is not quite right.

[Felix Richter]

Hello David,

next time I will add the whole m.file. However, initially I thought it would be an easy fix, because I made something wrong Matlab-wise. That's why I tried to keep the example easy.

Ah okay, yes I used SI-units, that could be the cause of my troubles. How ironic... :D

I substituted

r_core = R_core^2

so i get only three results.

Yes, I know you would not know what refers to what, I posted that as soon as it worked, I will fix this of course!

[David Goodmanson]

Hi Felix,

In terms of r_core this appears to be a cubic equation with real coefficients, but in that case it's not possible to get two real roots and one imaginary root as you have. Something is not quite right.

[Walter Roberson]

If you plot all six branches of the sextic then 2 of them have no real values and 4 of them have real values for large enough V. However this does not mean that when you plot the cubic that one of the branches will have no real values and the other 2 will have real values will have large enough V: you end up with all three having real values for sufficiently large V.

How can this be, that the sextic has branches with no real values at all, but the cubic has real values in all three branches? It is because one of three branches of the cubic is entirely negative, and with R = sqrt(r ) and r pure negative, that gives +/- pure imaginary roots.

[Felix Richter]

@David

For sqrt(r_core) I get three complex solutions. However, the first and second has substantially higher real values than their imaginary parts. For the third it is the opposite. That must be the right solution since it equals nearly exactly the book results.

@Walter

That may be the case, but P_core or V, respectively, are limited, due to real thermal and mechanical boundaries. A third solution at 1 Quintillion Megawatts would not be of interest. :D

[Walter Roberson]

The real solutions start at the same place on all four branches, at V = (44*sqrt(115))/529

[Felix Richter]

Okay, what does this mean in terms of applicability? Is there be a third/fourth combination of R_core and H_core, that is feasible?

[David Goodmanson]

Probably not, but it's hard to say without seeing the code that uses P_core. It appears that r_core satisfies a cubic equation with real coefficients. If that's the case, you can have either a) three real roots or b) one real root and two complex roots that are complex conjugates of each other. Your results for P_core = 1.5e9 don't fit either catagory.

[Walter Roberson]

When you state it in terms of the cubic, you get three real roots. However, one of them is all negative, and in the original R space corresponds to a pair of solutions that are purely complex (with no real component.)

So if you want to use the cubic, you should eliminate the negative root, as it is falsely implying a real-valued root of the original system.

Your plot will start being valid at V = (44*sqrt(115))/529; for smaller V, all roots of the sextic are complex-valued.

[Felix Richter]

@Walter Robertson

Alright, this does make sense!

Thank you all for your help and consider my problem solved.

(Although I might have to step up my math game a bit! :D)