# plot for complex equation: why is Mathematica plot different from Maple plot?

I have a complex equation:

$2 * B * m + 2 * c * v + 2 * M^2 * v - 4 * M5^2 * v - 4 * k1 * v^3 - 4 * k2 * v^3 = 0$

where $B, m, c, v, M, M5, k1, k2$ are real variables.

I want to solve it for variable $v$:

Solve[2 B m + 2 c v + 2 M^2 v - 4 M5^2 v - 4 k1 v^3 - 4 k2 v^3 == 0, v]


{{v -> (1.63607*10^-14 (-1.33546*10^31 +
5.89443*10^26 M5^2))/(-6.05149*10^46 + Sqrt[
3.66205*10^93 + 4. (-1.33546*10^31 + 5.89443*10^26 M5^2)^3])^(
1/3) - 1.03066*10^-14 (-6.05149*10^46 + Sqrt[
3.66205*10^93 + 4. (-1.33546*10^31 + 5.89443*10^26 M5^2)^3])^(
1/3)}, {v -> -(((8.18035*10^-15 +
1.41688*10^-14 I) (-1.33546*10^31 +
5.89443*10^26 M5^2))/(-6.05149*10^46 + Sqrt[
3.66205*10^93 + 4. (-1.33546*10^31 + 5.89443*10^26 M5^2)^3])^(
1/3)) + (5.1533*10^-15 - 8.92577*10^-15 I) (-6.05149*10^46 +
Sqrt[3.66205*10^93 +
4. (-1.33546*10^31 + 5.89443*10^26 M5^2)^3])^(
1/3)}, {v -> -(((8.18035*10^-15 -
1.41688*10^-14 I) (-1.33546*10^31 +
5.89443*10^26 M5^2))/(-6.05149*10^46 + Sqrt[
3.66205*10^93 + 4. (-1.33546*10^31 + 5.89443*10^26 M5^2)^3])^(
1/3)) + (5.1533*10^-15 + 8.92577*10^-15 I) (-6.05149*10^46 +
Sqrt[3.66205*10^93 +
4. (-1.33546*10^31 + 5.89443*10^26 M5^2)^3])^(1/3)}}


I know the approximate values of the parameters:

k1 = 16.485010961790245; k2 = -13.131344420001051; c = -44687.3983417778; B = 161593.81818181818; m = 5.5; M = 300;


And consequently I may plot real part of it:

Plot[{Re[(1.636069372813509*^-14 (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2))/(-6.051491252967043*^46 + \
\[Sqrt](3.662054638473664*^93 +
4. (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2)^3))^(1/3) -
1.0306591209480217*^-14 (-6.051491252967043*^46 + \
\[Sqrt](3.662054638473664*^93 +
4. (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2)^3))^(1/3)],
Re[-(((8.180346864067545*^-15 +
1.4168776392101725*^-14 I) (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2))/(-6.051491252967043*^46 + \
\[Sqrt](3.662054638473664*^93 +
4. (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2)^3))^(
1/3)) + (5.1532956047401085*^-15 -
8.92576981383125*^-15 I) (-6.051491252967043*^46 + \
\[Sqrt](3.662054638473664*^93 +
4. (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2)^3))^(1/3)],
Re[-(((8.180346864067545*^-15 -
1.4168776392101725*^-14 I) (-1.3354606752754323*^31 +

5.894433894342971*^26 M5^2))/(-6.051491252967043*^46 + \
\[Sqrt](3.662054638473664*^93 +
4. (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2)^3))^(
1/3)) + (5.1532956047401085*^-15 +
8.92576981383125*^-15 I) (-6.051491252967043*^46 + \
\[Sqrt](3.662054638473664*^93 +
4. (-1.3354606752754323*^31 +
5.894433894342971*^26 M5^2)^3))^(1/3)]}, {M5, 0, 400},
PlotStyle -> {{RGBColor[0.461492, 0.563303, 0.0104797], Dashed,
Thickness[0.003]}, {RGBColor[
0.4078757751993936, 0.27540579780035374, 0.780310562001516],
Dashed, Thickness[0.003]}, {RGBColor[0.65, 0., 0.], Dotted,
Thickness[0.003]}, {RGBColor[
0.7055434835026511, 0.5895048315002945, 0.], Dashed,
Thickness[0.003]}},
PlotLegends ->
Placed[{"first solution", "second solution", "third solution"},
Scaled[{0.8, 0.35}]]]


If I do the same thing in Maple:

solve(-4*k1*v^3-4*k2*v^3+4*(-M5)^2*v+2*M^2*v+2*C*v+2*b*m = 0, v)
k1 := 16.48501096179024
C := -44687.3983417778
k2 := -13.13134442000105
b := 161593.8181818181
M := 300
m := 5.5


Stackexchange has a limit for symbols, but in the construction of plot below I will show explicit form solutions from Maple:

plot([6.498027517*10^(-16)*(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)-1.538928540*10^15*(0.9939370219e-1*M5^2-2251.893618)/(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3), -3.249013759*10^(-16)*(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)+7.694642700*10^14*(0.9939370219e-1*M5^2-2251.893618)/(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)+(.8660254038*I)*(6.498027517*10^(-16)*(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)+1.538928540*10^15*(0.9939370219e-1*M5^2-2251.893618)/(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)), -3.249013759*10^(-16)*(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)+7.694642700*10^14*(0.9939370219e-1*M5^2-2251.893618)/(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)-(.8660254038*I)*(6.498027517*10^(-16)*(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3)+1.538928540*10^15*(0.9939370219e-1*M5^2-2251.893618)/(2.414698568*10^50+1.301965372*10^28*sqrt(7.694642700*10^31*M5^6-5.229964194*10^36*M5^4+1.184916421*10^41*M5^2-5.508856818*10^44))^(1/3))], M5 = 0 .. 400)


I get the next plot:

You may see two interesting things:

1) Part of plot that is in Mathematica, isn't in Maple.

2) Solves in Mathematica consist of pieces, but in Maple solves are continuous.

I want to understand where is problem and so my question is: do you know another method for plotting graph for complex equation in Mathematica?

• What makes us confident that Maple is correct, instead of Mathematica?
– ktm
Sep 2, 2017 at 22:42
• @user6014 I don't know what is correct. Sep 2, 2017 at 23:07
• @Nasser I did it. Sep 2, 2017 at 23:09
• Solve[2 B m + 2 c v + 2 M^2 v - 4 M5^2 v - 4 k1 v^3 - 4 k2 v^3 == 0, v, Reals] reproduces the Maple plot, except for the discontinuous branch jumping. In this case, reordering with Solve[2 B m + 2 c v + 2 M^2 v - 4 M5^2 v - 4 k1 v^3 - 4 k2 v^3 == 0 /. v -> -v, v, Reals] and plotting Re[-v] gets the Maple plot exactly. (The full solution set is the one originally produced by Mathematica.) Sep 2, 2017 at 23:29
• I am amazed Mathematica can do something sensible without encountering massive cancellation error due to those awkward exponents. In an example like this I would recommend setting 'Cubics->False' and perhaps also leaving the parameter values unspecified, substituting them after the fact. So something along these lines: solns = v /. (Solve[ 2 B m + 2 c v + 2 M^2 v - 4 M5^2 v - 4 k1 v^3 - 4 k2 v^3 == 0, v, Cubics -> False] /. {k1 -> 16.485010961790245, k2 -> -13.131344420001051, c -> -44687.3983417778, B -> 161593.81818181818, m -> 5.5, M -> 300}). Sep 3, 2017 at 14:44

First, there are three complex roots, counting multiplicities, for all values of the parameters. I don't know Maple nor why the command does not produce all of three. The Mathematica command shows all three (albeit for M5 > 150.2, the real parts of the complex-conjugate pair coincide).

Second, the discontinuities are Sign[x] == x / Abs[x] == x / Sqrt[x^2] type that sometimes arise with fractional powers. The problem is complicated enough that tracking down exactly whether and how it can be fixed seemed more work than might be worth it for the site. One way to see the problem with the first solution returned by Solve[] is to plot its numerator and denominator:

Plot[Flatten[
Re@Through[{Numerator, 30 Denominator[#] &}@
Together[v /. sol[[1]]]]] // Evaluate, {M5, 140, 160}]


One can see an absolute value type vee in the denominator where the numerator and denominator both vanish.

An alternative is Root. However, it sorts the roots and is usually a discontinuous function of the parameters, too. In this particular case, it can be used if we substitute v -> -v:

Block[{
k1 = 16.485010961790245,
k2 = -13.131344420001051,
c = -44687.3983417778,
B = 161593.81818181818,
m = 5.5,
M = 300},
Plot[Re@Table[
-Root[
Evaluate[2 B m + 2 c v + 2 M^2 v - 4 M5^2 v - 4 k1 v^3 - 4 k2 v^3 /. v -> -#] &,
n],
{n, 3}] // Evaluate,
{M5, 0, 400}]
]
`