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
E
as a variable name, since that's used by Mathematica for the constant $e$.) $\endgroup$