# 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[] = xG0b;
yGb[] = yG0b;
xGc[] = xG0;
yGc[] = yG0;
PosB[] = {xG0b, yG0b};
PosC[] = {xG0, yG0};
Distanza[] = ((yGb[] - yGc[])^2 + (xGb[] - xGc[])^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. – Michael E2 Aug 10 '16 at 13:41
• How could I deal with that? – usumdelphini Aug 10 '16 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, y} == ({x, y} /.
FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, 0.1}, {y, -1}]);
ics2 = {x, y} == ({x, y} /.
FindRoot[{f[x, y, 0] == 0, g[x, y, 0] == 0}, {x, -1}, {y, 1}]);
ics3 = {x, y} == ({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};

{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. – usumdelphini Aug 10 '16 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? – J. M. will be back soon Aug 10 '16 at 16:18
• @J.M. nope, but looks like a good read for an amateur like me - thanks! – Chris K Aug 10 '16 at 16:46
• No problem. From one amateur to another... ;) – J. M. will be back soon Aug 10 '16 at 16:47