1
$\begingroup$

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]

An answer is complex:

{{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}]]]

enter image description here

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: enter image description here

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?

|improve this question
$\endgroup$
  • 1
    $\begingroup$ What makes us confident that Maple is correct, instead of Mathematica? $\endgroup$ – user6014 Sep 2 '17 at 22:42
  • $\begingroup$ @user6014 I don't know what is correct. $\endgroup$ – illuminates Sep 2 '17 at 23:07
  • $\begingroup$ @Nasser I did it. $\endgroup$ – illuminates Sep 2 '17 at 23:09
  • 1
    $\begingroup$ 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.) $\endgroup$ – Michael E2 Sep 2 '17 at 23:29
  • 1
    $\begingroup$ 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}). $\endgroup$ – Daniel Lichtblau Sep 3 '17 at 14:44
3
$\begingroup$

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}]

Mathematica graphics

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}]
 ]

Mathematica graphics

|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.