NDSolve producing non-closed curves for a vector field

Question

I am solving a system of ODEs in the form

$$ \begin{aligned} \dot{x} (t) &= B_x (x, z) \\ \dot{z} (t) &= B_y (x, z) \end{aligned} $$ where $\vec{B}$ is defined as follows $$ \vec{B} (x, z) = \int_0^{2 \pi} \mathrm{d} \phi \frac{(z \cos \phi, 1 - x \cos \left( \phi \right)}{\left[ \left( x - \cos \phi \right)^2 + \sin^2 \phi + z^2 \right]^{3/2}} $$

To set this thing up in Mathematica I defined

Bx[r_?NumericQ, φ_?NumericQ, z_?NumericQ] := 
     NIntegrate[(
      z Cos[ϕ])/((r Cos[φ] - 
          Cos[ϕ])^2 + (r Sin[φ] - Sin[ϕ])^2 + z^2)^(
      3/2), {ϕ, 0, 2 π}]
    Bz[r_?NumericQ, φ_?NumericQ, z_?NumericQ] := 
     NIntegrate[(
      1 - r Cos[ϕ - φ])/((r Cos[φ] - 
          Cos[ϕ])^2 + (r Sin[φ] - Sin[ϕ])^2 + z^2)^(
      3/2), {ϕ, 0, 2 π}]

(note that this also contains angle $\varphi$ which comes from polar coordinates, however, the problem is rotationally symmetric around $z$-axis, so in reality, everything what's interesting in this is $(r, z)$ behaviour for some constant $\varphi$. When $\varphi = 0$ then $r = x$)

Then I proceeded with

b = 5; x0 = 0.8;
sol = First@
  NDSolve[{x'[t] == Bx[Abs[x[t]], 0, z[t]], 
    z'[t] == Bz[Abs[x[t]], 0, z[t]], x[0] == x0, z[0] == 0, 
    WhenEvent[x[t] < 0 || x[t] > b || Abs[z[t]] > b, 
     "StopIntegration"]}, {x, z}, {t, -100, 100}]

Please note that the integrals for $\vec{B}$ are not expressible in terms of elementary functions. But even if it were, in future I will surely use functions that are not expressible in terms of elementary functions.

The problem is:

As you can see, the curve is not closed. That result is nonsensical as it would indicate that at the point of intersection, the field is multivalued ($\vec{B}$ corresponds to two different vectors). I think this is why all lines of vector fields are either closed or they run from infinity to infinity, or they end up in some "sink" point in space. But they just cannot intersect like this. In fact, this particular curve should be closed, but due to some numerics, it isn't. My goal is to find the culprit and fix the numerics to obtain closed curves (or curves that vanish because of running out of image) for various starting points.

Answer 1

The key is setting a finite AccuracyGoal. (When the initial integral is zero or very close to zero, neither the accuracy goal of Infinity nor any precision goal can be met. A lot of time can be spent in needless refinement.) Another suggestion is to use the Trapezoidal method since the integrand is periodic and the interval is a period.

Bx[r_?NumericQ, φ_?NumericQ, z_?NumericQ] := 
  NIntegrate[(z Cos[ϕ])/((r Cos[φ] - Cos[ϕ])^2 + (r Sin[φ] - Sin[ϕ])^2 + z^2)^(3/2),
   {ϕ, 0, 2 π}, Method -> "Trapezoidal", AccuracyGoal -> 12];
Bz[r_?NumericQ, φ_?NumericQ, z_?NumericQ] := 
  NIntegrate[(1 - r Cos[ϕ-φ])/((r Cos[φ] - Cos[ϕ])^2 + (r Sin[φ] - Sin[ϕ])^2 + z^2)^(3/2),
   {ϕ, 0, 2 π}, Method -> "Trapezoidal", AccuracyGoal -> 12];

PrintTemporary@Dynamic@{foo, Clock[Infinity]};
Block[{b = 5, x0 = 0.8},
 sol = First@
   NDSolve[{x'[t] == Bx[Abs[x[t]], 0, z[t]], 
     z'[t] == Bz[Abs[x[t]], 0, z[t]], x[0] == x0, z[0] == 0, 
     WhenEvent[x[t] < 0 || x[t] > b || Abs[z[t]] > b, 
      "StopIntegration"]}, {x, z}, {t, -10, 10}, 
    StepMonitor :> (foo = {t, Abs[z[t]]})]
 ]

ListLinePlot[{x, z} /. sol]

P.S. Using phi and curly-phi seems likely to be error-prone. But to each one's own.

Update: As can be seen, the curves are closed:

Update 2: I'm still not sure what's being asked. If one cycle only is desired, then the trouble is how WhenEvent works. It seems that the event detection criterion (first argument of WhenEvent) is applied only at the steps. Once an event is detected, it is then located. So at the point of detection, the return of the solution will be very roughly approximately equal to its initial value. After the event is located, it will be much closer.

NDSolve[{x'[t] == Bx[Abs[x[t]], 0, z[t]], 
  z'[t] == Bz[Abs[x[t]], 0, z[t]], x[0] == x0, z[0] == 0,
  WhenEvent[z[t] == 0 && echo@Norm[{x[t], z[t]} - {x0, 0}] < 10^-1, 
   Print[{t, x[t], z[t], Norm[{x[t], z[t]} - {x0, 0}]}]; 
   If[Norm[{x[t], z[t]} - {x0, 0}] < 10^-7, "StopIntegration"]]},
 {x, z}, {t, 0, 3}]

Concerning "how to raise accuracy so that the point it returns to is within the machine precision the same numerically it was at the start?": You have to have enough extra precision and sufficiently high precision and accuracy goals so that the truncation and rounding errors give a results with an error of 1 ulp. (I'm not sure what error is satisfactory for z[t], since the initial value was 0. The smallest error is either $MinMachineNumber or in V11.3, the smallest subnormal number 5.*10^-324.) In any case, WorkingPrecision -> 34 seems to give an error of 0., when the solution is compared to {0.8, 0}.