1
$\begingroup$

I created a nice Manipulate box of the magnetic field outside and inside a spherical shell placed in an homogeneous magnetic field. The code is working great, but it has a problem at a limiting value, and the performance is pretty low (slow). I don't mind much about the slow execution, but I need to find a way to prevent the problem at the limiting value.

Here's a reduction of the code which shows the problem :

MagnAxis = {0, 1};
X[s_] := {x[s], z[s]}
r[s_] := Norm[X[s]]

dipField[s_] := 3(MagnAxis.X[s])X[s]/r[s]^5 - MagnAxis/r[s]^3

Coef0[u_, alpha_] := Which[(alpha < 1), 1, (alpha == 1), 0]
Coef1[u_, alpha_] := 3(3 - alpha)
Coef2[u_, alpha_] := - 3alpha u^3
Coef3[u_, alpha_] := 9(1 - alpha) + 2alpha^2(1 - u^3)
Coef4[u_, alpha_] := alpha(3 - alpha)(1 - u^3)

FieldOrientation[s_, u_, alpha_] := 
  If[u < 1,
    Piecewise[{
      {Coef0[u, alpha]MagnAxis, (0 <= r[s] < u)},
      {Normalize[Coef1[u, alpha]MagnAxis + Coef2[u, alpha]dipField[s]], 
         (u <= r[s] < 1)},
      {Normalize[Coef3[u, alpha]MagnAxis + Coef4[u, alpha]dipField[s]], 
         (r[s] >= 1)}}],
    MagnAxis]

Nlines = 40;

MagnCurve[u_, alpha_, n_] := 
  NDSolve[{
   x'[s] == {1, 0}.FieldOrientation[s, u, alpha],
   z'[s] == {0, 1}.FieldOrientation[s, u, alpha],
   x[0] == 0.15(2n - 1 - Nlines)/2,
   z[0] == -3}, 
  {x, z}, {s, 0, 8},
  Method -> BDF, MaxSteps -> 1000000]

MagnGraph[u_, alpha_, n_] := 
  ParametricPlot[
    Evaluate[{x[s], z[s]}/.MagnCurve[u, alpha, n]], {s, 0, 8},
    PlotStyle -> {Directive[Blue]},
    PerformanceGoal -> "Quality"]

Shell[u_] :=  
  Graphics[
    {{Black, Thick, Circle[{0, 0}, 1]},
     {LightGray, Disk[{0, 0}, 1]},
     {Black, Thick, Circle[{0, 0}, u]},
     {White, Disk[{0, 0}, u]}}]

Manipulate[
  Show[
    Shell[u],
    Table[MagnGraph[u, alpha, n], {n, 1, Nlines}],
    PlotRange -> {{-zoom, zoom}, {-zoom, zoom}},
    ImageSize -> 600],
  {{u, 1/3, Style[Subscript[R, 1]/Subscript[R, 2], 10]}, 0, 1, 0.01, 
    ImageSize -> Large, Appearance -> {"Labeled", "Closed"}},
  {{alpha, 0.5, Style["Susceptibility", 10]}, -5, 1, 0.01, 
    ImageSize -> Large, Appearance -> {"Labeled", "Closed"}},
  Delimiter,
  {{zoom, 2, Style["Zoom", 10]}, 1, 3, 0.1, 
    ImageSize -> Large, Appearance -> {"Labeled", "Closed"}}]

Now, set the Susceptibility slider to 1. The field lines are okay (but the output is slow to get). Then if you set the R1/R2 slider to 1, the field lines turn into simple straight vertical lines, which implies that there is no shell anymore (when the inside radius equals the outside radius : R1/R2 = 1). This is normal behavior.

But when you reduce the R1/R2 slider a bit, the Manipulate box just stop and cancels the calculation. Or you may get a message like this:

NDSolve::ndcf : Repeated convergence test failure at s = ....'; unable to continue.

Ideally, the field lines should be of the same type as what you get with Susceptibility set to 1 and R1/R2 set to 0.9999 (directly into the value boxes).

Here's a preview picture of a very thin shell :

thin shell

So what is wrong with my code above? And how can I improve the performance? Using PerformanceGoal -> "Speed" doesn't much change the execution speed, and the field lines get ugly in some places.

$\endgroup$
6
  • $\begingroup$ Your code doesn't work as is, probably because you did not include the definition of Coordonnees[s]. Please make sure that your posted code can run on its own by copying it out of your post and trying to run it in a fresh Mathematica notebook. We can't help without complete running code. $\endgroup$
    – MarcoB
    Feb 29, 2016 at 16:28
  • $\begingroup$ @MarcoB, oops! sorry, this is a small mistake made to define a MWE from the whole project. I fixed that in the question. Try it again, it should work now. $\endgroup$
    – Cham
    Feb 29, 2016 at 18:10
  • 1
    $\begingroup$ You should define FieldOrientation so that it is only evaluated numerically. Your code doesn't work for me at all by the way (v10.1) so I cant really confirm that will help, but I think so. $\endgroup$
    – george2079
    Feb 29, 2016 at 22:43
  • 2
    $\begingroup$ Cham - you're code does not run as it is... please check again using a new kernel... $\endgroup$
    – bill s
    Mar 1, 2016 at 0:19
  • $\begingroup$ I don't understand why it isn't working for you. I tried the code above in a fresh new session of Mathematica, and the code is working. It is pretty slow however. Do you get any error or warning messages ? How to define the ** FieldOrientation** numerically ? $\endgroup$
    – Cham
    Mar 1, 2016 at 1:03

1 Answer 1

3
+50
$\begingroup$

In V10 (and probably V9), the problem is that NDSolve has trouble projecting onto the discontinuity when the annulus is too thin. (It looks like it tries to step across it.) The trick below is to set an event just outside the annulus to restart the integration with a very small step size. This seems to allow NDSolve to properly detect the discontinuity.

MagnAxis = {0, 1};
X[s_] := {x[s], z[s]}
r[s_] := Norm[X[s]]

dipField[s_] := 3 (MagnAxis.X[s]) X[s]/r[s]^5 - MagnAxis/r[s]^3

Coef0[u_, alpha_] := Piecewise[{{1, alpha < 1}}, 0]
Coef1[u_, alpha_] := 3 (3 - alpha)
Coef2[u_, alpha_] := -3 alpha u^3
Coef3[u_, alpha_] := 9 (1 - alpha) + 2 alpha^2 (1 - u^3)
Coef4[u_, alpha_] := alpha (3 - alpha) (1 - u^3)

FieldOrientation[s_, u_, alpha_] := Piecewise[{
   {Piecewise[{
      {Coef0[u, alpha] MagnAxis, (0 <= r[s] < u)},
      {Normalize[Coef1[u, alpha] MagnAxis + Coef2[u, alpha] dipField[s]], (u <= r[s] < 1)}},
     Normalize[Coef3[u, alpha] MagnAxis + Coef4[u, alpha] dipField[s]]
     ], 
    u < 1}},
  MagnAxis]

Nlines = 40;

MagnCurve[u_, alpha_, n_] := NDSolve[{
   x'[s] == {1, 0}.FieldOrientation[s, u, alpha], 
   z'[s] == {0, 1}.FieldOrientation[s, u, alpha],
   x[0] == 0.15 (2 n - 1 - Nlines)/2, z[0] == -3,
   WhenEvent[x[s]^2 + z[s]^2 == 1 + 1*^-8, "RestartIntegration"]},
  {x, z}, {s, 0, 8}, StartingStepSize -> 1*^-8]

MagnGraph[u_, alpha_, n_] := 
 ParametricPlot[
  Evaluate[{x[s], z[s]} /. MagnCurve[u, alpha, n]], {s, 0, 8}, 
  PlotStyle -> {Directive[Blue]}, PerformanceGoal -> "Quality"]

Shell[u_] := 
 Graphics[{{Black, Thick, Circle[{0, 0}, 1]}, {LightGray, 
    Disk[{0, 0}, 1]}, {Black, Thick, Circle[{0, 0}, u]}, {White, 
    Disk[{0, 0}, u]}}]

Output:

Manipulate[
 Show[
  Shell[u],
  Table[MagnGraph[u, alpha, n], {n, 1, Nlines}],
  Frame -> True,
  PlotRange -> {{-zoom, zoom}, {-zoom, zoom}}],

 {{u, 1/3, Style[Subscript[R, 1]/Subscript[R, 2], 10]}, 0, 1, 0.01, 
  Appearance -> {"Labeled", "Closed"}},
 {{alpha, 0.5, Style["Susceptibility", 10]}, -5, 1, 0.01, 
  Appearance -> {"Labeled", "Closed"}},
 Delimiter,
 {{zoom, 2, Style["Zoom", 10]}, 1, 3, 0.1, Appearance -> {"Labeled", "Closed"}}]

Mathematica graphics

Update: Replacing

Table[MagnGraph[u, alpha, n], {n, 1, Nlines}],

with

Graphics[{
  Blue,
  Line@Table[{x["ValuesOnGrid"], z["ValuesOnGrid"]} /. 
      First@MagnCurve[u, alpha, n] // Transpose, {n, 1, Nlines}]
  }]

saves about 0.8 sec. (out of 2.16) on my machine. If I reduce the PrecisionGoal in NDSolve, to

PrecisionGoal -> 4

I save another 0.2 sec., and the plot is still fairly accurate. It still takes over a second to compute and render the plot.

$\endgroup$
14
  • $\begingroup$ It isn't woking well with Mathematica 7, apparently. Is this WhenEvent compatible with v7 ? $\endgroup$
    – Cham
    Mar 1, 2016 at 2:25
  • $\begingroup$ @Cham I no longer have access to V7. The functionality was called EventLocator in V7 and a similar trick should be possible with it. I can't reproduce your problem exactly, since I don't have V7, so I'm not sure exactly what problem V7 has with your system, or whether this approach will address that problem. $\endgroup$
    – Michael E2
    Mar 1, 2016 at 2:32
  • $\begingroup$ The problem is when you set the susceptibility to 1, and R1/R2 to 1 too, then move a bit the R1/R2 slider a bit to the left. The computation may stop and cancels itself. I would like to prevent this to happen. Or in other words, R1/R2 = 1 is a kind of singularity. I would like to contournate it in some way. $\endgroup$
    – Cham
    Mar 1, 2016 at 2:37
  • $\begingroup$ @Cham I understood that, but what I mean is that I can only reproduce something like your problem in V10 and solve it in V10. I can't reproduce in V7 and search for solutions in V7 (because I don't have V7). I only hope my approach might help out in V7. $\endgroup$
    – Michael E2
    Mar 1, 2016 at 2:43
  • 1
    $\begingroup$ @Cham It's the opposite in V10 (PrecisionGoal is for NDSolve, not for ParametricPlot). The V7 docs suggest it's the same in V7. $\endgroup$
    – Michael E2
    Mar 1, 2016 at 14:01

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.