Plot all real and imaginary roots of FindRoot (2)

Question

Following the comment of our friend @Bob Hanlon's in my post Plot all real and imaginary roots of FindRoot vrs `k`, which I would like to thank again for the assistance it gave, I'm creating a new post here. As I clarified with @BobHanlon, my problem is to solve a complicated algebraic equation including an integral and then plot all the real Re[x] and imaginary Im[x] parts of the roots x versus k, t and a parameters.

Equation to solve is

eq[x_,k_,t_,a_] := 1 -( k^2/x^2) -((t^2) /(2 BesselK[2, t] k^2))
NIntegrate[s Sqrt[s^2-1] ((1 + (a/x)) - (x /k) Sqrt[s^2/(s^2 - 1)] ((1 + (a/x))^2 - ( 
(s^2 - 1) k^2)/(s^2 x^2))ArcTanh[(Sqrt[(s^2 - 1)/s^2] k)/(x(1 + (a/x )))]) Exp[-t s], 
{s,1, \[Infinity]}]

with 0.001 <= t <= 1000, 0 <= k <= 100 and a = 25

To my knowledge, I suppose that this equation can only be solved by FindRoot (SolveValues gives nothing), which I do not master well, especially the choice of the initial value and precisions of the computations, ...

Here is my code

rootx[k_,t_,a_]:=FindRoot[1 -( k^2/x^2) -((t^2) /(2 BesselK[2, t] k^2)) 
NIntegrate[s Sqrt[s ^2 -1] ((1 + (a/x)) - (x /k) Sqrt[s^2/(s^2 - 1)] ((1 + (a/x))^2 - (( 
s^2 - 1) k^2)/(s^2 x^2)) ArcTanh[(Sqrt[(s^2 - 1)/s^2] k)/(x(1 + (a/x)))]) Exp[-t s], 
{s,1, \[Infinity]}]== 0, {x, 1 + 0.01 I}, WorkingPrecision -> $MachinePrecision]


(* Real root of x *)
rex[k_,t_,a_]:=Re[rootx[k,t,a]]
(* rex gives a one real root for a given (k,t,a) but I'm trying to plot all real roots *)

(* Imaginary root of x *)
imx[k_,t_,a_]:=Im[rootx[k,t,a]]
(* imx gives a one imaginary root for a given (k,t,a) but I'm trying to plot all imaginary roots *)

Please, how can I know if the initial value x=1+0.01 I is the right value who gives the right roots? and how to plot all the Re and Im roots vrs k,t,a ?

THANKS

Answer 1

$Version

(* "13.2.0 for Mac OS X ARM (64-bit) (November 18, 2022)" *)

Clear["Global`*"]

eq[x_?NumericQ, k_?NumericQ, t_?NumericQ,
  a_?NumericQ, wp_ : 10] :=
 1 - (k^2/x^2) - ((t^2)/(2 BesselK[2, t] k^2))*
   NIntegrate[
    s Sqrt[s^2 - 
       1]* ((1 + (a/x)) - (x/k) Sqrt[
         s^2/(s^2 - 1)]*((1 + (a/x))^2 - ((s^2 - 1) k^2)/(s^2 x^2))*
        ArcTanh[(Sqrt[(s^2 - 1)/s^2] k)/(x (1 + (a/x)))]) Exp[-t s], {s, 
     1, ∞}, WorkingPrecision -> wp]

To get estimates of the roots (starting values for use with FindRoot) for specific values of {k, t}, look at the plot of Abs[eq]^2

Manipulate[
 Module[{k, t},
  k = Rationalize[kv, 0];
  t = Rationalize[tv, 0];
  ComplexPlot3D[
   {Abs[eq[x, k, t, 25]]^2, 0},
   {x, -10 - 5*I, 10 + 5*I},
   PlotRange -> {-1/2, 1},
   PlotStyle -> Opacity[0.75],
   AxesLabel -> {"Re", "Im", None},
   ClippingStyle -> None,
   PlotPoints -> 50,
   MaxRecursion -> 3,
   SphericalRegion -> True]],
 {{kv, 3, "k"}, 0.01, 100, 0.5,
  Appearance -> "Labeled"},
 {{tv, 10, "t"}, 0.001, 1000, 5,
  Appearance -> "Labeled"},
 SynchronousUpdating -> False,
 TrackedSymbols :> {kv, tv}]

The roots are where the upper curves intersect the plane at z == 0

EDIT: Manipulating the plots above suggests that the roots are located near x == k and x == -k with small if any imaginary parts.

root[k_?NumericQ, t_?NumericQ, a_?NumericQ] :=
 x /. FindRoot[eq[x, k, t, a, 20] == 0, {x, #}] & /@
  {-k + I/10, k + I/10}

SetOptions[RandomReal, WorkingPrecision -> 20];

SeedRandom[1234];

rgn = ImplicitRegion[10^-3 <= t <= 1000 && 0 < k <= 100, {k, t}];

Generate random points in the {k, t} plane

pts = RandomPoint[rgn, 200];

roots = Transpose[
    Thread[{#[[1]], #[[2]], root[#[[1]], #[[2]], 25]}] & /@ pts] // Quiet;

Show[
 ListPlot3D[Re[roots],
  AxesLabel -> (Style[#, 14] & /@
     {k, t, "Re[roots]"}),
  PlotLabel -> "a = 25"],
 ListPointPlot3D[roots,
  PlotStyle -> Red],
 ImageSize -> Medium]

rootsIm = Map[{#[[1]], #[[2]], Im[#[[3]]]} &, roots, {2}];

Most of the imaginary component of the roots are negligible

MinMax[#[[All, 3]]] & /@ rootsIm

(* {{-0.000172782, 0.000508208}, 
    {-5.41156*10^-14, 5.58909*10^-13}} *)

Show[
 ListPlot3D[rootsIm,
  PlotRange -> All,
  AxesLabel -> (Style[#, 14] & /@
     {k, t, "Im[roots]"}),
  PlotLabel -> "a = 25"],
 ListPointPlot3D[rootsIm,
  PlotStyle -> Red],
 ImageSize -> Medium]

EDIT:

kvalues = Range[0, 100, 2] /. {0 -> 10^-3};

data = With[{t = 10, a = 25}, 
    Transpose@Table[Thread[{k, root[k, t, a]}], 
       {k, kvalues}]] // Quiet;

ListLinePlot[Re@data]

ListLinePlot[data /. Complex[r_, i_] :> i, PlotRange -> All]