2

The goal is to find the stationary points of the function $$f(x)=\frac{y}{x-4\times10^{-5}}-\frac{1}{x^2}$$ for $y=7000$ using Mathematica. The input used was

f[x_] := (y)/(x - 4*^-5) - 1/x^2
y = 7000
Solve[D[f[x]] == 0, x]
Plot[f[x], {x, 8*^-5, 20*^-5}]

For $1<y<3500$, Solve outputs two real solutions as expected. For larger values of say $y>7000$ then Solve outputs complex-valued solutions, indicating that D[f[x]]=0 has no real solutions. However, the plot shows that there are two real solutions, although only a precise level of zooming will show this.

Is the code incorrect for this purpose? Or is Mathematica just not able to solve the equation?

David Dinh
  • 121
  • 1
  • 2
    f[x_] := (y)/(x - k) - 1/x^2 y = 7000 Solve[D[f[x]] == 0, x] You get two roots (1 - Sqrt[1 - 28000 k])/14000 and (1 + Sqrt[1 - 28000 k])/14000 but if 28000 * k is bigger than 1 so the square roots have negative numbers under them, hence if k is 410^-5 we get Sqrt[1 - 280004*10^-5] which is Sqrt[-(3/25)], and you get complex values. – flinty May 14 '23 at 10:37
  • i.e Mathematica is able to solve it fine. It's the plot that's deceiving you. The more weird thing is if you set y as high as it can go (y = 1/(4k)) while still allowing a real solution. This is 6250. Mathematica then returns 'two' solutions {{x -> 1/12500}, {x -> 1/12500}} which are exactly the same. – flinty May 14 '23 at 10:48
  • 1
    D[f[x]] doesn't do anything. If you want a derivative with respect to $x$, you need to specify it: D[f[x], x]. – Roman May 14 '23 at 11:14

1 Answers1

5

I think it's easier to look at this problem in terms of exact polynomial roots instead of numerical solutions.

The equation $\partial f/\partial x=0$ is a cubic polynomial in $x$ and therefore has three roots (real or complex). By expressing these roots as Root objects we can ensure that we catch them all, for any combination of $y$ and $\epsilon$:

f[ε_, x_, y_] = y/(x - ε) - 1/x^2;
z[ε_, y_] = SolveValues[D[f[ε, x, y], x] == 0, x, Cubics -> False]
(*    {Root[-2 ε^2 + 4 ε #1 - 2 #1^2 + y #1^3 &, 1],
       Root[-2 ε^2 + 4 ε #1 - 2 #1^2 + y #1^3 &, 2],
       Root[-2 ε^2 + 4 ε #1 - 2 #1^2 + y #1^3 &, 3]}    *)

where I've written $\epsilon$ instead of $4\times10^{-5}$.

Now we plot them:

Plot[z[4*^-5, y], {y, 0, 10000}]

enter image description here

we see that there are three real roots for $y\lesssim7500$ and one real root above. This is quite common for cubic polynomials. The exact cutoff is found from the zero of the polynomial discriminant:

Solve[Discriminant[-2 ε^2 + 4 ε #1 - 2 #1^2 + y #1^3 &[x], x] == 0, y]
(*    {{y -> 0}, {y -> 8/(27 ε)}}

8/(27 ε) /. ε -> 4.^-5
(    7407.41    *)
Roman
  • 47,322
  • 2
  • 55
  • 121