3
$\begingroup$

The following stripped down code works, but I'm experiencing some strange warnings: from an old document made with Mathematica 7.0 (but used in version 13.2), I'm getting several warning messages about "Event location failed to converge to requested accuracy or precision...". The output is okay, though. When I copy paste the exact same code in a new document made with Mathematica 13.2, I don't get any warning/error message, but the StoppingTest is drawn in red. So the issue lies there (indicated with the comment (* ISSUE IS HERE! *) in the code below).

MomMagn = {1, 1, 1};  (* norm of each magnetic moment *)
thetaMu = {45, 45, 0}Pi/180;  (* inclinaison relative to Z axis *)
phiMu = {90, 270, 0}Pi/180;  (* azimutal angle around Z axis *)

source[k_] := {{3, 0, 0}, {-8, 0, 0}, {0, 8, 0}}[[k]] (* position of each magnetic sphere *)

rotationX[k_] := {{1, 0, 0}, {0, Cos[thetaMu[[k]]], -Sin[thetaMu[[k]]]}, {0, Sin[thetaMu[[k]]], Cos[thetaMu[[k]]]}} (* rotation matrix around X axis *)
rotationZ[k_] := {{Cos[phiMu[[k]]], -Sin[phiMu[[k]]], 0}, {Sin[phiMu[[k]]], Cos[phiMu[[k]]], 0}, {0, 0, 1}}  (* rotation matrix around Z axis *)
mu[k_] := MomMagn[[k]]rotationZ[k].rotationX[k].{0, 0, 1} (* magnetic moments *)

Nspheres = Length[MomMagn]; (* number of magnetic spheres *)

(* magnetic field definition: *)
xyz[s_] := {x[s], y[s], z[s]}
r[s_, k_] := Norm[xyz[s] - source[k]]
dipole[s_, k_] := 3 mu[k].(xyz[s] - source[k])(xyz[s] - source[k])/r[s, k]^5 - mu[k]/r[s, k]^3
magneticfield[s_] := Sum[dipole[s, k], {k, 1, Nspheres}]

(* magnetic field lines definition: *)
Nlines = 15;

theta[n_] := Which[n == -2, 20, n == -1, 0, n == 0, 0, n == 1, 20, n == 2, 27.5]Pi/180
phi[n_, m_] := If[n > 0, m 2 Pi/Nlines, 0]
sphericalstart[n_, m_] := If[n < 0, -1, 1] {Sin[theta[n]] Cos[phi[n, m]], Sin[theta[n]] Sin[phi[n, m]], Cos[theta[n]]}
linestart[n_, m_, k_] := rotationZ[k].rotationX[k].sphericalstart[n, m] + source[k]

magneticCurves[n_, m_, k_] := NDSolve[{
    x'[s] == If[n < 0, -1, 1]{1, 0, 0}.magneticfield[s]/Norm[magneticfield[s]],
    y'[s] == If[n < 0, -1, 1]{0, 1, 0}.magneticfield[s]/Norm[magneticfield[s]],
    z'[s] == If[n < 0, -1, 1]{0, 0, 1}.magneticfield[s]/Norm[magneticfield[s]],

    x[0] == {1, 0, 0}.linestart[n, m, k],
    y[0] == {0, 1, 0}.linestart[n, m, k],
    z[0] == {0, 0, 1}.linestart[n, m, k]
    }, {x, y, z}, {s, 0, 500},
    Method -> Automatic, MaxSteps -> Automatic,
    StoppingTest -> (Apply[Or, Table[EuclideanDistance[xyz[s], source[kk]] < 1, {kk, 1, Nspheres}]] || r[s, k] > 25)] (* ISSUE IS HERE! *)

Do[magneticCurves[n, m, k], {k, 1, Nspheres}, {n, -2, 2}, {m, 1, Nlines}]

Smin[n_, m_, k_] := (x/.magneticCurves[n, m, k])[[1]][[1]][[1]][[1]] (* start of each field line *)
Smax[n_, m_, k_] := (x/.magneticCurves[n, m, k])[[1]][[1]][[1]][[2]] (* end of each field line *)

magneticGraph[n_, m_, k_] := ParametricPlot3D[Evaluate[xyz[s]/.magneticCurves[n, m, k]], {s, Smin[n, m, k], Smax[n, m, k]}]

Show[Table[magneticGraph[n, m, k], {n, -2, 2}, {m, 1, Nlines}, {k, 1, Nspheres}], PlotRange -> All]

Preview of what the code is doing:

enter image description here

So what should I use instead of StoppingTest? If I remove completely this part, the output isn't right at all. The field lines output is right when I still use the red StoppingTest line. I don't understand what's going on here.

Note: This code is slow (on a Silicon Mac M2 Pro with 32GB ram). This is NOT a concern to me. The issue isn't about the compilation speed and I don't need to optimize the whole code.

$\endgroup$
0

1 Answer 1

6
$\begingroup$

StoppingTest isn't available anymore in the newer versions of Mathematica. You should use WhenEvent or "EventLocator" method instead.

For example:

magneticCurves[n_, m_, k_] := 
 With[{conds = Or@@Table[EuclideanDistance[xyz[s], source[kk]] < 1,
   {kk, 1, Nspheres}] || r[s, k] > 25}, 
  NDSolve[{x'[s] == 
     If[n < 0, -1, 1] {1, 0, 0} . magneticfield[s]/Norm[magneticfield[s]], 
    y'[s] == 
     If[n < 0, -1, 1] {0, 1, 0} . magneticfield[s]/Norm[magneticfield[s]], 
    z'[s] == 
     If[n < 0, -1, 1] {0, 0, 1} . magneticfield[s]/Norm[magneticfield[s]], 
    x[0] == {1, 0, 0} . linestart[n, m, k],
    y[0] == {0, 1, 0} . linestart[n, m, k], 
    z[0] == {0, 0, 1} . linestart[n, m, k],
    WhenEvent[conds, "StopIntegration"]
  }, 
  {x[s], y[s], z[s]}, 
  {s, 0, 500}, 
  Method -> Automatic, MaxSteps -> Automatic]]
$\endgroup$
10
  • $\begingroup$ Why is the cond event declared inside a With environment, and not the other conditions? $\endgroup$
    – Cham
    Jul 28, 2023 at 17:42
  • $\begingroup$ @Cham With applies to the whole block $\endgroup$ Jul 28, 2023 at 17:43
  • $\begingroup$ @Cham, because r[s, k] has to be evaluated to Sqrt[...] > 25 before it is fed into NDSolve. $\endgroup$
    – Domen
    Jul 28, 2023 at 17:43
  • $\begingroup$ @Domen, then what about the r[s, k] < 1 condition? Why not put it alongside with >25? I don't understand this part $\endgroup$
    – Cham
    Jul 28, 2023 at 17:46
  • $\begingroup$ @Cham, I have simplified the code to make it more clear. Let me know if there is still anything unclear. $\endgroup$
    – Domen
    Jul 28, 2023 at 18:11

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.