# Using FindRoot to track moving Bessel Functions zeros

Setting up the problem

I want to study the common zeros of these two functions:

r[x_, y_] := Sqrt[x^2 + y^2];
Θ[x_, y_] := ArcTan[x, y];

f[x_, y_, t_] := 1/2 - BesselK[0, r[x, y]]/(2 BesselK[0, 1]) + ϵ1[t] BesselK[0, r[x, y]]/BesselK[0, 1] + D1[t] BesselK[1, r[x, y]]/(2 BesselK[1, 1])    Cos[Θ[x, y] + ζ1[t]]

g[x_, y_, t_] :=  Sin[Θ[x, y] + ζ2[t]]/  r[x, y] + ϵ2/2 r[x, y]^(-1/2) Cos[Θ[x, y]/2]


I assign values to the time-dependent parameters:

ϵ1[t_] = -0.5 + 10^-6 t;
ζ2[t_] = π/3 + 10^-4 t;
ζ1[t_] = π/3 + 10^-6 t;
ϵ2 = 0.5;
D1[t_] = 0.5 - 10^-3 t;


And I start with t=0. A CountourPlot shows that there are two zeros:

ContourPlot[{g[x, y, 0] == 0, f[x, y, 0] == 0}, {x, -10, 10}, {y, -10,10}]


Their positions are:

  xG0 = x /. FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, 0.1}, {y, -1},        AccuracyGoal -> 5]
yG0 = y /. FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, 0.1}, {y, -1},        AccuracyGoal -> 5]

(* 0.378089 *)
(* -1.26037 *)


And:

  xG0b =  x /. FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, -1}, {y, 1},       AccuracyGoal -> 5]
yG0b = y /. FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, -1}, {y, 1},       AccuracyGoal -> 5]

(* -1.02897 *)
(* 1.31148 *)


Tracking the zeros:

Now, I want to be able to track their positions as the variable t (time) increases.

In order to do this, I have devised a simple For loop with discretised time-stepping, which I Break when the two zeros get close enough.

At each time step, I look for the i-th zero in a neighbourhood of the previous one, so that each individual zero is not mistaken with the other one.

(* number of steps *)
Nf=100;

(* Final time *)
Tf = 10000;

(* Time Step *)
dt=Tf/Nf;

(* arrays with the positions and distance *)
xGc = Table[0, {i, Nf + 1}];
yGc = Table[0, {i, Nf + 1}];
xGb = Table[0, {i, Nf + 1}];
yGb = Table[0, {i, Nf + 1}];
PosC = Table[{0, 0}, {i, Nf + 1}];
PosB = Table[{0, 0}, {i, Nf + 1}];
Distanza = Table[0, {i, Nf + 1}];

(* I assign the first elements - initial positions *)
xGb[[1]] = xG0b;
yGb[[1]] = yG0b;
xGc[[1]] = xG0;
yGc[[1]] = yG0;
PosB[[1]] = {xG0b, yG0b};
PosC[[1]] = {xG0, yG0};
Distanza[[1]] = ((yGb[[1]] - yGc[[1]])^2 + (xGb[[1]] - xGc[[1]])^2)^(1/2);

(* Search region tolerance*)
xtol=0.1;
ytol=xtol;

(* For Loop *)

For[i = 2, i <= Nf + 1, i = i + 1, t = (i - 1) dt;

(* I look for the i-th zero in a neighbourhood of the previous one *)

xGc[[i]] =
x /. FindRoot[{f[x, y, t] == 0, g[x, y, t] == 0}, {x,
xGc[[i - 1]], xGc[[i - 1]] - xtol, xGc[[i - 1]] + xtol}, {y,
yGc[[i - 1]], yGc[[i - 1]] - ytol, yGc[[i - 1]] + ytol},
AccuracyGoal -> 2, PrecisionGoal -> 2];
yGc[[i]] =
y /. FindRoot[{f[x, y, t] == 0, g[x, y, t] == 0}, {x,
xGc[[i - 1]], xGc[[i - 1]] - xtol, xGc[[i - 1]] + xtol}, {y,
yGc[[i - 1]], yGc[[i - 1]] - ytol, yGc[[i - 1]] + ytol},
AccuracyGoal -> 2, PrecisionGoal -> 2];

xGb[[i]] =
x /. FindRoot[{f[x, y, t] == 0, g[x, y, t] == 0}, {x,
xGb[[i - 1]], xGb[[i - 1]] - xtol, xGb[[i - 1]] + xtol}, {y,
yGb[[i - 1]], yGb[[i - 1]] - ytol, yGb[[i - 1]] + ytol},
AccuracyGoal -> 2, PrecisionGoal -> 2];
yGb[[i]] =
y /. FindRoot[{f[x, y, t] == 0, g[x, y, t] == 0}, {x,
xGb[[i - 1]], xGb[[i - 1]] - xtol, xGb[[i - 1]] + xtol}, {y,
yGb[[i - 1]], yGb[[i - 1]] - ytol, yGb[[i - 1]] + ytol},
AccuracyGoal -> 2, PrecisionGoal -> 2];

Distanza[[i]] =
N[((yGb[[i]] - yGc[[i]])^2 + (xGb[[i]] - xGc[[i]])^2)^(1/2)];

PosB[[i]] = {xGb[[i]], yGb[[i]]};
PosC[[i]] = {xGc[[i]], yGc[[i]]};

If[Distanza[[i]] <= 0.5,
Tf = t;
Nf = i;
xGc[[Nf + 1]] = xGc[[i]];
xGb[[Nf + 1]] = xGb[[i]];
yGb[[Nf + 1]] = yGb[[i]];
yGc[[Nf + 1]] = yGc[[i]];
Break[]
]
]


I get various errors:

FindRoot::reged: The point {-0.350794,0.695484} is at the edge of the search region {0.695484,0.895484} in coordinate 2 and the computed search direction points outside the region. >>

FindRoot::reged: The point {-0.350794,0.695484} is at the edge of the search region {0.695484,0.895484} in coordinate 2 and the computed search direction points outside the region. >>

FindRoot::reged: The point {-0.290278,0.595484} is at the edge of the search region {0.595484,0.795484} in coordinate 2 and the computed search direction points outside the region. >>

General::stop: Further output of FindRoot::reged will be suppressed during this calculation. >>

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

General::stop: Further output of FindRoot::lstol will be suppressed during this calculation. >>


And when I plot some of the position variables (here Distanza and xGb), I see some suspect jumps:

Questions:

1) Why do I see the errors, and how could I improve the situation?

2) Is there a way to implement an adaptive Search-region size, so that if I get an error I increase xtol and ytol until I find the zero within the right accuracy?

3) Would something like the extrapolation of a guess from the two previous zeros help more? How would I go about implementing that?

• There's a singularity (some kind of node) that develops near t == 1744 or so. Aug 10, 2016 at 13:41
• How could I deal with that? Aug 10, 2016 at 13:42

I found three roots for 0 < t < ~1745 and only one root after that. I use NDSolve[] and DAEs to trace the equations. There are several examples on the site.

ics1 = {x[0], y[0]} == ({x, y} /.
FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, 0.1}, {y, -1}]);
ics2 = {x[0], y[0]} == ({x, y} /.
FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, -1}, {y, 1}]);
ics3 = {x[0], y[0]} == ({x, y} /.
FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, 0.01}, {y, -0.01}]);
dae = {Simplify[                                    (* remove discontinuity of ArcTan *)
{f[x[t], y[t], t], g[x[t], y[t], t]} /.
Cos[θ_/2] :> Sqrt[(1 + Cos[θ])/2] /. tr_Cos | tr_Sin :> TrigExpand@tr
] == {0, 0}, s'[t] == 1, s[0] == 0};

{sol1} = NDSolve[{dae, ics1}, {x, y}, {t, 0, 10000},
"ExtrapolationHandler" -> {Nothing &, "WarningMessage" -> False}];
{sol2} = NDSolve[{dae, ics2}, {x, y}, {t, 0, 10000},
"ExtrapolationHandler" -> {Nothing &, "WarningMessage" -> False}];
{sol3} = NDSolve[{dae, ics3}, {x, y}, {t, 0, 10000},
"ExtrapolationHandler" -> {Nothing &, "WarningMessage" -> False}];


NDSolve::ndsz: At t == 1744.8567484346531, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 1744.8567484340776, step size is effectively zero; singularity or stiff system suspected.

The warning are the side-effect of two roots coalescing and disappearing.

Manipulate[
Show[
ContourPlot[{g[x, y, t] == 0, f[x, y, t] == 0},
{x, -2.5, 2.5}, {y, -2.5, 2.5}],
Graphics[{Red, PointSize[Medium],
Point[{x[t], y[t]} /. {sol1, sol2, sol3} /. {} -> Nothing]}]
],
{t, 0, 10000}
]


Note that the apparent intersection at the origin is actually a pole.

Plot3D[f[x, y, 600] // TrigExpand // Simplify //
Evaluate, {x, -0.00002, 0.00002}, {y, -0.00002, 0.00002},
MeshFunctions -> {Function[{x, y, z}, g[x, y, 600]]}, Mesh -> {{0}},
MeshShading -> {Lighter@Blue, Lighter@Red}, WorkingPrecision -> 32,
PlotLabel -> "Plot of f at t = 600, colored by the sign of g"]


• Thanks, this was very helpful. I am now trying with another set functions (the ones posted here were a minimal example) and I get a series of errors, but above all: NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended. Aug 10, 2016 at 20:49

The pseudo-arclength continuation function TrackRootPAL I hacked together here works on this problem. First, define TrackRootPAL from that link. Then start at your two initial points to get two tracks of roots:

tr = TrackRootPAL[{f[x, y, t], g[x, y, t]}, {x, y}, {t, 0, 10000},
0, {xG0, yG0}, SMax -> 10000];
trb = TrackRootPAL[{f[x, y, t], g[x, y, t]}, {x, y}, {t, 0, 10000},
0, {xG0b, yG0b}, SMax -> 10000];


Note I had to increase the range of arclength (SMax) up to 10000.

Plot the results:

Plot[Evaluate[{x[t], y[t]} /. tr], {t, 0, 10000}]
Plot[Evaluate[{x[t], y[t]} /. trb], {t, 0, 10000}]
Show[
ParametricPlot3D[{t, x[t], y[t]} /. tr, {t, 0, 10000}, PlotRange -> All],
ParametricPlot3D[{t, x[t], y[t]} /. trb, {t, 0, 1744}, PlotRange -> All],
BoxRatios -> {1, 0.25, 0.25}, ImageSize -> Large
]


The second track trb folds back on itself around t=1744 as noted by @MichaelE2.

• Have you seen Allgower/Georg, by any chance? Aug 10, 2016 at 16:18
• @J.M. nope, but looks like a good read for an amateur like me - thanks! Aug 10, 2016 at 16:46
• No problem. From one amateur to another... ;) Aug 10, 2016 at 16:47