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]};
            {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.

enter image description here

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

enter image description here

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

  • $\begingroup$ This doesn't address the Mathematica issues here, but it looks like you're calculating field line diagrams for dipole radiation in the $xz$-plane. It turns out that these field lines are expressible as the level sets of another function in the $xz$-plane. See Zangwill's Modern Electrodynamics, §20.5.2, & references therein. $\endgroup$ Jun 22, 2015 at 16:03
  • $\begingroup$ Mathematica seems to do just fine with the field lines of a static dipole, integrating past the points where $E_z = 0$ without any problem. Can you post your code for the electric field so that potential answerers can toy around with it? (Also, as an aside: it's bad practice to use E as a variable name, since that's used by Mathematica for the constant $e$.) $\endgroup$ Jun 22, 2015 at 17:03
  • $\begingroup$ @Michael: Thanks for the replies. I also tried out the field lines of a static dipole, and it worked fine for me either! But for the field I posted above, the integration wont get beyond a certain parameter s. Thus I thought it is more Mathematica related than Physics related, because I guess I just have to pass another option to NDSolve in order to get the thing working (and also because it already did work in a previous Mathematica version). And as you have suggested, I have also posted the mathematica expression for the field equation above. $\endgroup$
    – juergen
    Jun 23, 2015 at 9:31

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

enter image description here

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.

  • $\begingroup$ Thank you! This solved the problem indeed! I have quit often encountered situations, where Mathematica was quite picky with squareroots, but I didn't have thought that the parametrization condition would be a problem! Thanks again! $\endgroup$
    – juergen
    Jun 24, 2015 at 10:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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