# Plot all real and imaginary roots of FindRoot (2)

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 ## 1 Answer $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},

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

• Thank you @BobHanlon, but how to plot all the Re[x] and Im[x] roots? Commented Feb 8, 2023 at 5:03
• The computation time of starting values for specific values of {k,t}, @BobHanlon, is very long . In addition, I did not succeed in plotting the roots. Commented Feb 8, 2023 at 15:55
• Sorry, I have to admit I'm a bit confused, I couldn't plot all the Re[x] and Im[x] roots. Commented Feb 9, 2023 at 17:28
• Thank you so much @Bob Hanlon for this great wonderful job you are doing. Last thing please, how to plot Re[x] as a function of k and Im[x] as a function of k for 0<=k<=100 and for a fixed value of t, for example t=10`? Commented Feb 10, 2023 at 2:27