0
$\begingroup$

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

$\endgroup$

1 Answer 1

1
$\begingroup$
$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}]

enter image description here

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]

enter image description here

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]

enter image description here

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]

enter image description here

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

enter image description here

$\endgroup$
8
  • $\begingroup$ Thank you @BobHanlon, but how to plot all the Re[x] and Im[x] roots? $\endgroup$
    – Gallagher
    Commented Feb 8, 2023 at 5:03
  • $\begingroup$ Use the plots to come up with estimates to use in FindRoot. $\endgroup$
    – Bob Hanlon
    Commented Feb 8, 2023 at 5:12
  • $\begingroup$ 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. $\endgroup$
    – Gallagher
    Commented Feb 8, 2023 at 15:55
  • $\begingroup$ Sorry, I have to admit I'm a bit confused, I couldn't plot all the Re[x] and Im[x] roots. $\endgroup$
    – Gallagher
    Commented Feb 9, 2023 at 17:28
  • $\begingroup$ 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? $\endgroup$
    – Gallagher
    Commented Feb 10, 2023 at 2:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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