0
$\begingroup$

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: Solution

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.

$\endgroup$
2
$\begingroup$

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]

Mathematica graphics

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:

Mathematica graphics

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}.

$\endgroup$
  • $\begingroup$ Perhaps the OP is more interested in the WhenEvent than the integration time? I'm not sure... $\endgroup$ – Michael E2 Apr 14 '18 at 1:01
  • $\begingroup$ That's nice, but the curves are still not closed (see my edit for reference). $\endgroup$ – user16320 Apr 14 '18 at 1:30
  • $\begingroup$ @user16320 It's unclear what the curve is you plotted. $\endgroup$ – Michael E2 Apr 14 '18 at 1:53
  • $\begingroup$ What version are you running? After executing your code, the curve is closer to be closed, but it is still a bit off on my machine. Also, I had to add a WhenEvent to stop integration after it returns to the starting point, but it copies the curve twice anyway (in each direction). How to add a condition so that it ends whenever it visits the vicinity of the same point he visited before? My solution ((x[t] - x0)^2 + z[t] < 10^-4 && Abs@t > 0.1) does not work universally (when the curve is turning far away, it ends in the middle of nowhere)... $\endgroup$ – user16320 Apr 14 '18 at 6:16
  • $\begingroup$ To put it simply: how to raise accuracy so that the point it returns to is within the machine precision the same numerically it was at the start? $\endgroup$ – user16320 Apr 14 '18 at 6:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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