NDSolve-Problem (Calculating field lines)

Question

I want to calculate the field lines of a quite nasty vectorfield. In general the method works fine - except, one of the components of the VF gets 0 - then NDSolve complains about

NDSolve::mxst: Maximum number of 300 steps reached at the point s == 0.06617737920434953`

The weird thing is, that the same routine worked fine in Mathematica 7, but starts complaining in newer Mathematica versions.

The routine I use, is the following:

IntegrateFieldLines[{Ex_, Ez_}, {x0_, z0_}, length_] :=
   Module[{Exs, Ezs},
          Exs[s_] := Ex /. {x -> x[s], z -> z[s]};
          Ezs[s_] := Ez /. {x -> x[s], z -> z[s]};
          NDSolveValue[
            {x'[s] Ezs[s] == z'[s] Exs[s],
             x'[s]^2 + z'[s]^2 == 1,
             x[0] == x0,
             z[0] == z0},
            {x, z}, {s, 0, length}]]

Thus, I want to solve $E_x(x,z)/E_y(x,z) = dx/dy$, where x and z are paramatrised by $s$.

The Vectorfield I want to compute is quite nasty

$$\vec E(x,y,z)=Re\left(\frac{e^{i (k r - \omega t)}}{r} \left(k \vec n \times (\vec p \times \vec n )+\frac{(1-i k r)}{r ^2}\left(3 \vec n (\vec p \vec n) - \vec p)\right)\right)\right),$$ where $\omega t=\pi$, $k=2 \pi /5, r=\sqrt{x^2+y^2+z^2}$.

If you want to play around with it, just copy the Mathematica formula below,

    Evec[x_, y_, z_] := Re[Exp[I (k r - \[Omega] t)]/
r (k^2 n\[Cross](p\[Cross]n) + (1 - I k r)/
   r^2 (3 n Dot[p, n] - p))] /. {k -> 2 \[Pi]/5, \[Omega] t -> \[Pi], p -> {0, 0, 0.2}, n -> {x, y, z}/Sqrt[x^2 + y^2 + z^2], r -> Sqrt[x^2 + y^2 + z^2]}

And I only want to get the field lines in the xz-plane, therefore I call

IntegrateFieldLines[{Evec[x, 0, z][[1]], Evec[x, 0, z][[3]]}, {4, 0}, 10]

if I want to calculate the field line starting at (4,0).

Like I said before, it works fine until NDSolve reaches MaxSteps. In the following picture, one can see the calculated field lines (in red) and a StreamLinePlot - but I can't manage to calculate the full field line.

I have found at, that the point s, where NDSolve stops working, is where $E_z(x(s),z(s)) == 0$:

For solving my problem, I have already tried out a smaller number in AccuarcyGoal, as well as a higher number in MaxSteps... but neither worked. Maybe one of you has an idea, how I get NDSolve to integrate further and why Mathematica 7 could do the job, but all newer versions fail.

Greetings, Jürgen

Answer 1

If you look at the output of your original data, you'll notice that you're getting small imaginary parts in your results. This implies to me that Mathematica is trying to take square roots, which is why it's failing when $E_z$ is going to zero; and the main culprit, I believe, is your constraint that $x'(s)^2 + z'(s)^2 = 1$. Let's differentiate that constraint explicitly to eliminate the possibility of square roots:

$$ x''(s) x'(s) + z''(s) z'(s) = 0 $$

This equation implies that $x'(s)^2 + z'(s)^2$ is a constant $C$, but not necessarily 1. We can add in an initial condition that this constant be 1 to ensure this: $$ x'(0)^2 + z'(0)^2 = 1 $$

So we rewrite your code as follows:

IntegrateFieldLines[{Ex_, Ez_}, {x0_, z0_}, length_] := 
 Module[{Exs, Ezs}, Exs[s_] := Ex /. {x -> x[s], z -> z[s]};
  Ezs[s_] := Ez /. {x -> x[s], z -> z[s]};
  NDSolveValue[{x'[s] Ezs[s] == z'[s] Exs[s], x''[s] x'[s] + z''[s] z'[s] == 0, 
                 x'[0]^2 + z'[0]^2 == 1, x[0] == x0, z[0] == z0},
               {x, z}, {s, 0, length}]]

Running through the same procedure, the equations can now be integrated for the full time:

My guess as to why your original code worked in previous versions of Mathematica but not the present one is that the previous versions defaulted to a method that didn't have this problem with these particular equations, but that the newer versions pick a different method that fails. This would imply (to me) that there's problem some Method option that would allow your original set of equations to be integrated as well. I played around with a few of them, but I couldn't get any of them to work; you might dig around in the NDSolve documentation if you're interested in getting your original equations to work again.