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.
– Michael E2
Aug 10 at 13:41

  

 

How could I deal with that?
– usumdelphini
Aug 10 at 13:42

=================

2 Answers
2

=================

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.
– usumdelphini
Aug 10 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.♦
Aug 10 at 16:18

  

 

@J.M. nope, but looks like a good read for an amateur like me – thanks!
– Chris K
Aug 10 at 16:46

  

 

No problem. From one amateur to another… 😉
– J. M.♦
Aug 10 at 16:47