Solve using PolyGamma function

Question

I am trying to solve for the value of y.

My code is:

x = 0.165; 

f = -Log[y] - PolyGamma[0.5 + 0.2*(x/y)] + PolyGamma[0.5]; 

Solve[f == 0, y]

Running gave me

Solve::inex: Solve was unable to solve the system with inexact coefficients or >the system obtained by direct rationalization of inexact numbers present in the >system. Since many of the methods used by Solve require exact input, providing >Solve with an exact version of the system may help.

Using MATLAB gave me a value of y.

Answer 1

Solve is a symbolic solver of algebraic as well as various types of transcendental equations and consequently we should use exact numbers.

x = Rationalize[0.165];
f[y_] := -Log[y] - PolyGamma[1/2 + 1/5 (x/y)] + PolyGamma[1/2]

Next we should restrict the domain of y, a natural assumption is y > 0, nevertheless one should also add a bound from above (this is not always necessary). Taking a look at e.g. Plot[f[y], {y, 0, 5}], we should be satisfied with 0 < y < 5, and so Solve provides an exact solution in terms of the Root object:

y /. First @ FullSimplify[ Solve[ f[y] == 0 && 0 < y < 5, y]]

Root[{EulerGamma + Log[4] + Log[#1] + PolyGamma[0, 1/2 + 33/(1000*#1)]&,
      0.8324921392608209131327584577663713541314494518823960255129`30
      .102999566398122}]

This solution can be further transformed symbolically for more detailed purposes. Using a bound for y we could use NSolve as well without rationalizing inexact coefficients to get a numerical solution, however an exact solution is a full information we can get from the given equation. Our solution can be represented in terms of HarmonicNumber if we appropriately simplify f, e.g.

fs = FullSimplify[ f[y], y > 0]

-HarmonicNumber[-(1/2) + 33/(1000 y)] - Log[4 y]

y /. First @ Solve[fs == 0 && 0 < y < 5, y]

Root[{HarmonicNumber[-1/2 + 33/(1000*#1)] + Log[4*#1]&,
      0.8324921392608209131327584577663713541314494518823960255129`
      30.102999566398122}]

The assumption that y > 0 yields a real solution. There are also complex solutions, without specifying them we can get an insight where they can be found with ContourPlot, e.g.

ContourPlot[{ Re[f[u + I w]] == 0, Im[f[u + I w]] == 0}, 
            { u, -0.095, 0.005}, {w, -0.07, 0.07}, 
            PerformanceGoal -> "Quality", ContourStyle -> Thick, 
            AspectRatio -> Automatic]

and slightly more detailed neighbourhood of the singular point (it takes a while)

ContourPlot[{ Re[f[u + I w]] == 0, Im[f[u + I w]] == 0}, 
            {u, -0.023, -0.003}, {w, -0.008, 0.008}, 
            PlotPoints -> 50, MaxRecursion -> 4, PerformanceGoal -> "Quality", 
            ContourStyle -> Thick, AspectRatio -> Automatic]