0
$\begingroup$

My aim

I want to find the solution of a cubic equation as a function of a parameter. In my code, the cubic equation is to be solved for the variable th with Pdyn as the parameter. Then, I want to plot the solutions (there will be three) with respect to varying Pdyn. Below is what I have tried.

Code and outputs

k1 = 10; 
k2 = 2;   
S = 100;
c1 = 2*Pi;
c2 = 1; 
alpha = 0;
e = 2/10; 
c = 5;  

eqncub = 
  k1*th + k2*th^3 - Pdyn*S*e*c*(c1*(alpha + th) - c2*(alpha + th)^3) == 0;
theq = Solve[eqncub, th][[All, All, 2]];
th2 = theq[[2]][[1]]
Plot[th2, {Pdyn, 1/20/Pi, 20}]

Note the below result.

Solve[th2 == 0, Pdyn]
{{Pdyn -> 1/(20*Pi)}}

However, my plot, which I specify that Pdyn should run from $ \frac{1}{20\pi}$ to $20$, doesn't even start from zero! Below is the plot.

Second solution, th2, of cubic equation plotted versus parameter Pdyn

A few more notes

The solution calculated by Mathematica is correct. I have verified it using hand calculation.

The solution, th2, is clearly only real for Pdyn $\geq \frac{1}{20\pi}$. However, Mathematica seems to have no problem giving me plots for any values of Pdyn, even negative values. This is weird.

Edit/Update

Just about 10 minutes after posting this, I tried another thing and it makes me believe the problem lies solely when plotting the function.

f[Pdyn_] := Evaluate[th2];
Plot[f[Pdyn], {Pdyn, 1/20/Pi, 20}]

So here, f[1/20/Pi] gives me zero but the plot again starts nowhere near zero.

$\endgroup$
4
  • 2
    $\begingroup$ Use AxesOrigin -> {0, 0} in your plot options. $\endgroup$
    – Ben Izd
    Feb 21 at 12:36
  • 1
    $\begingroup$ Or: Plot[th2, {Pdyn, 1/20/Pi, 20}, PlotRange -> All] $\endgroup$ Feb 21 at 12:37
  • $\begingroup$ @BenyIzd It works! Why do I need to do this though? $\endgroup$
    – ModCon
    Feb 21 at 12:45
  • $\begingroup$ @MariuszIwaniuk This also works! But for some reason, the y-axis crosses x-axis at x=0.05 and so my plot crosses y-axis, which can be misleading. Any idea why this happens? $\endgroup$
    – ModCon
    Feb 21 at 12:47
3
$\begingroup$

Generally the answer is AxesOrigin -> {0, 0}, but sometimes PlotRange -> All as Mariusz Iwaniuk mentioned might works as well.

Regard your second question about Why Mathematica changes the plot range automatically and how it determines?, is above my knowledge to answer it comprehensively but this is what I find:

AxesOrigin documentation states:

In 2D graphics, AxesOrigin -> Automatic uses an internal algorithm to determine where the axes should cross. If the point {0,0} is within, or close to, the plotting region, then it is usually chosen as the axes origin.

Mathematica automatically determined how to show your plot even sometimes it might cross your ranges. It will do its best to show the best picture, instead of showing drastic changes, it might crop a little and show a beautiful curve like yours 😉.

Some plots start from $n$, so Mathematica usually crop the range to not waste any space in the plot if is away from (0,0):

ListLinePlot[Table[{x, x^2}, {x, 5, 10, .1}]]

enter image description here

It almost plots from (5,0), the plot with (0,0) origin is this:

enter image description here

But if you're close to (0,0) it almost doesn't change:

Plot[x^15, {x, 0, 1}]

enter image description here

Plot[x^15, {x, 0.5, 1}, AxesOrigin -> {0, 0}]

enter image description here

But for some reason with the function below, internal algorithm determines to change the origin:

Plot[Surd[x, 10], {x, 0, 50}]

enter image description here

Plot[Surd[x, 10], {x, 0, 50}, AxesOrigin -> {0, 0}]

enter image description here

$\endgroup$
1
  • $\begingroup$ Makes sense. Thank you very much for your answer! $\endgroup$
    – ModCon
    Feb 24 at 11:23
2
$\begingroup$

Plot will default to a PlotRange that shows what it believes is the likely region of interest. If the default is not your region of interest, then you must override the default by specifying the PlotRange.

Clear["Global`*"]

k1 = 10;
k2 = 2;
S = 100;
c1 = 2*Pi;
c2 = 1;
alpha = 0;
e = 2/10;
c = 5;

eqncub = k1*th + k2*th^3 - 
    Pdyn*S*e*c*(c1*(alpha + th) - c2*(alpha + th)^3) == 0;
theq = Solve[eqncub, th] // Simplify

(* {{th -> 0}, {th -> -(Sqrt[-5 + 100 Pdyn π]/Sqrt[1 + 50 Pdyn])}, {th -> 
   Sqrt[-5 + 100 Pdyn π]/Sqrt[1 + 50 Pdyn]}} *)

Verifying the solutions,

eqncub /. theq // Simplify

(* {True, True, True} *)

For Pdyn < -1/50 the argument of both radicals are negative and the ratio is real, i.e., the imaginary part of the ratio is zero.

Assuming[Pdyn < -1/50, 
 Im[th /. theq] // ComplexExpand[#, TargetFunctions -> {Re, Im}] & // 
  Simplify]

(* {0, 0, 0} *)

Plotting the solutions,

Plot[Evaluate[th /. theq], {Pdyn, -1/2, 1/2},
 PlotLegends -> "Expressions",
 Epilog -> {Red, Thin,
   InfiniteLine[{{-1/50, -1}, {-1/50, 1}}],
   InfiniteLine[{{1/20/Pi, -1}, {1/20/Pi, 1}}]}]

enter image description here

$\endgroup$
1
  • $\begingroup$ I completely missed that Pdyn < -1/50 gives a real number! Thank you very much for your entire answer. $\endgroup$
    – ModCon
    Feb 24 at 11:23

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