0
$\begingroup$

I am new in Mathematica, and this is the first time I coding. Following is the code solving differential equations of motions in spherical coordinates of magnetic pendulum and plotting basins of attraction. I got idea from this code:

X[1] = 1;
X[2] = -(1/2);
X[3] = -(1/2);
Y[1] = 0;
Y[2] = Sqrt[3]/2;
Y[3] = -Sqrt[3]/2;
nf = Nearest[{{1, 0}, {-0.5, Sqrt[3]/2}, {-0.5, -Sqrt[3]/2}} -> 
     Automatic] /* First;
getStateData[k_, c_, h_, x0_, y0_] := 
 First@NDSolve`ProcessEquations[{x''[t] + k x'[t] + c x[t] - 
      Sum[(X[i] - 
          x[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i, 
        3}] == 0, 
    y''[t] + k y'[t] + c y[t] - 
      Sum[(Y[i] - 
          y[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i, 
        3}] == 0, x'[0] == 0, y'[0] == 0, x[0] == x0, y[0] == y0}, {x,
     y}, t, Method -> "Adams"]

sd = getStateData[.15, .2, .2, 1, 1];

getBasin2[x0_, y0_] := 
 Module[{state = sd, sol}, 
  state = First@NDSolve`Reinitialize[state, {x[0] == x0, y[0] == y0}];
  NDSolve`Iterate[state, 100.5];
  sol = {x[100], y[100]} /. NDSolve`ProcessSolutions[state];
  nf[sol]]
ArrayPlot[
  ParallelTable[
   getBasin2[xpos, ypos], {xpos, -2, 2, 0.1}, {ypos, -2, 2, 0.1}], 
  ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue}] // AbsoluteTiming

And following is my code.(It's quite long)

mags = {{1.2, 0, 0, 1, 0, 0, -1}, {1.2, Pi/6, 0, 1, -Sin[Pi/6], 
    0, -Cos[Pi/6]}, {1.2, Pi/6, 2*Pi/3, 
    1, -Sin[Pi/6]*Cos[2*Pi/3], -Sin[Pi/6]*
     Sin[2*Pi/3], -Cos[Pi/6]}, {1.2, Pi/6, 4*Pi/3, 
    1, -Sin[Pi/6]*Cos[4*Pi/3], -Sin[Pi/6]*Sin[4*Pi/3], -Cos[Pi/6]}};

nf = Nearest[{{-0.082788, 0}, {0.041394, 
       0.071696}, {0.041394, -0.071696}} -> Automatic] /* First;

getStateData[k_, c_, M_, m_, L_, g_, tht0_, phi0_] := 
  First@NDSolve`ProcessEquations[{tht''[t] == 
      Sin[tht[t]]*Cos[tht[t]]*
        phi'[t]^2 - (Sum[
           c*mags[[i, 
              4]]*((mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] + 
                 mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] - 
                 mags[[i, 7]]*
                  Sin[tht[t]]) + (3*
                   mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Cos[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
                    Sin[phi[t]] - 
                    Cos[mags[[i, 2]]]*
                    Sin[tht[
                    t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 7]]*Cos[tht[t]]) - 
                    mags[[i, 
                    1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]])) - 
                  3*(mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] + 
                    mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] - 
                    mags[[i, 7]]*Sin[tht[t]])*(1 - 
                    mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))/(((1 + 
                    mags[[i, 1]]^2 - 
                    2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))^1) - 
               15*mags[[i, 
                  1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                   Cos[tht[t]]*Cos[phi[t]] + 
                  Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
                   Sin[phi[t]] - 
                  Cos[mags[[i, 2]]]*
                   Sin[tht[t]])*((mags[[i, 5]]*Sin[tht[t]]*
                    Cos[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 7]]*Cos[tht[t]]) - 
                  mags[[i, 
                    1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 - 
                   mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]]))/((1 + 
                    mags[[i, 1]]^2 - 
                    2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))^2)/((1 + 
                 mags[[i, 1]]^2 - 
                 2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))^1.5, {i, 
            4}] + (m + M/2)*g*L*Sin[tht[t]] + 
          k*(L^2)*tht'[t]/3)/((M/3 + m)*L^2), 
     phi''[t] == -Sum[
          c*mags[[i, 
             4]]*((-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] + 
                mags[[i, 6]]*Sin[tht[t]]*
                 Cos[phi[t]]) + (3*
                  mags[[i, 
                    1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Sin[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Cos[phi[
                    t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 7]]*Cos[tht[t]]) - 
                    mags[[i, 
                    1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]])) - 
                 3*(-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Cos[phi[t]])*(1 - 
                    mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))/((1 + 
                   mags[[i, 1]]^2 - 
                   2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]]))^1) - 
              15*mags[[i, 
                 1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                  Sin[tht[t]]*Sin[phi[t]] + 
                 Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*

                  Cos[phi[t]])*((mags[[i, 5]]*Sin[tht[t]]*
                    Cos[phi[t]] + 
                   mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                   mags[[i, 7]]*Cos[tht[t]]) - 
                 mags[[i, 
                    1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 - 
                  mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + Cos[mags[[i, 2]]]*Cos[tht[t]]))/(1 +
                    mags[[i, 1]]^2 - 
                   2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]]))^2)/(1 + 
               mags[[i, 1]]^2 - 
               2*mags[[i, 
                  1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                   Sin[tht[t]]*Cos[phi[t]] + 
                  Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                   Sin[phi[t]] + 
                  Cos[mags[[i, 2]]]*Cos[tht[t]]))^1.5, {i, 
           4}]/(((M/3 + m)*L^2)*Sin[tht[t]]^2) - 
       2*tht'[t]*phi'[t]/Tan[tht[t]] - 
       k*(L^2)*phi'[t]/(3*((M/3 + m)*L^2)), tht'[0] == 0, 
     phi'[0] == 0, tht[0] == tht0, phi[0] == phi0}, {tht, phi}, t, 
    Method -> "Adams"];

sd = getStateData[0.1, 0.0004, 0.05, 0.01, 0.3, 9.8, 1, 3];

getBaSin[tht0_, phi0_] := 
 Module[{state = sd, sol}, 
  state = First@
    NDSolve`Reinitialize[state, {tht[0] == tht0, phi[0] == phi0}];
  NDSolve`Iterate[state, 110.];
  sol = {tht[100], phi[100]} /. NDSolve`ProcessSolutions[state];
  nf[{0.3*Sin[sol[[1]]]*Cos[sol[[2]]], 
    0.3*Sin[sol[[1]]]*Sin[sol[[2]]]}]]

ArrayPlot[
  ParallelTable[
   getBaSin[tht, phi], {tht, 0.1, 1.5, 0.1}, {phi, 0.1, 6.2, 0.1}], 
  ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue}] // AbsoluteTiming

First I want to know what the red lines mean, which says such like:

NDSolveReinitialize::dvnoarg : -- Message text not found -- (0.1),

since I expected the result as picture colored evenly with three colors.

And second I want to know how the result will differ if I change the NDSolve method(ex)Adams -> ExplicitRungeKutta), and which method will be most appropriate.

$\endgroup$
2
$\begingroup$

I fixed the code for my taste. There's only a typo when calling getBaSin[tht, phi]. You can not use the parameters tht, phi. We can reduce the time to t = 20. Fig.1 is different from what is obtained on the plane. In my opinion it is even more beautiful.

mags = {{1.2, 0, 0, 1, 0, 0, -1}, {1.2, Pi/6, 0, 1, -Sin[Pi/6], 
    0, -Cos[Pi/6]}, {1.2, Pi/6, 2*Pi/3, 
    1, -Sin[Pi/6]*Cos[2*Pi/3], -Sin[Pi/6]*
     Sin[2*Pi/3], -Cos[Pi/6]}, {1.2, Pi/6, 4*Pi/3, 
    1, -Sin[Pi/6]*Cos[4*Pi/3], -Sin[Pi/6]*Sin[4*Pi/3], -Cos[Pi/6]}};

nf = Nearest[{{-0.082788, 0}, {0.041394, 
       0.071696}, {0.041394, -0.071696}} -> Automatic] /* First;
eq = {tht''[
      t] - (Sin[tht[t]]*Cos[tht[t]]*
        phi'[t]^2 - (Sum[
           c*mags[[i, 
             4]]*((mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] + 
                 mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] - 
                 mags[[i, 7]]*
                  Sin[tht[t]]) + (3*
                   mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Cos[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
                    Sin[phi[t]] - 
                    Cos[mags[[i, 2]]]*
                    Sin[tht[
                    t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 7]]*Cos[tht[t]]) - 
                    mags[[i, 
                    1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]])) - 
                  3*(mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] + 
                    mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] - 
                    mags[[i, 7]]*Sin[tht[t]])*(1 - 
                    mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))/(((1 + 
                    mags[[i, 1]]^2 - 
                    2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))^1) - 
               15*mags[[i, 
                 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Cos[tht[t]]*
                   Cos[phi[t]] + 
                  Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
                   Sin[phi[t]] - 
                  Cos[mags[[i, 2]]]*
                   Sin[tht[t]])*((mags[[i, 5]]*Sin[tht[t]]*
                    Cos[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 7]]*Cos[tht[t]]) - 
                  mags[[i, 
                    1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 - 
                   mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]]))/((1 + 
                    mags[[i, 1]]^2 - 
                    2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))^2)/((1 + 
                 mags[[i, 1]]^2 - 
                 2*mags[[i, 
                   1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))^1.5, {i, 
            4}] + (m + M/2)*g*L*Sin[tht[t]] + 
          k*(L^2)*tht'[t]/3)/((M/3 + m)*L^2)) == 0, 
   phi''[t] - (-Sum[
          c*mags[[i, 
            4]]*((-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] + 
                mags[[i, 6]]*Sin[tht[t]]*
                 Cos[phi[t]]) + (3*
                  mags[[i, 
                   1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Sin[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*

                    Cos[phi[
                    t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 7]]*Cos[tht[t]]) - 
                    mags[[i, 
                    1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]])) - 
                 3*(-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] + 
                    mags[[i, 6]]*Sin[tht[t]]*Cos[phi[t]])*(1 - 
                    mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]])))/((1 + 
                   mags[[i, 1]]^2 - 
                   2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]]))^1) - 
              15*mags[[i, 
                1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Sin[tht[t]]*
                  Sin[phi[t]] + 
                 Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                  Cos[phi[t]])*((mags[[i, 5]]*Sin[tht[t]]*
                    Cos[phi[t]] + 
                   mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] + 
                   mags[[i, 7]]*Cos[tht[t]]) - 
                 mags[[i, 
                   1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
                    Cos[mags[[i, 3]]] + 
                    mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
                     mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 - 
                  mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + Cos[mags[[i, 2]]]*Cos[tht[t]]))/(1 +
                    mags[[i, 1]]^2 - 
                   2*mags[[i, 
                    1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
                    Sin[tht[t]]*Cos[phi[t]] + 
                    Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                    Sin[phi[t]] + 
                    Cos[mags[[i, 2]]]*Cos[tht[t]]))^2)/(1 + 
               mags[[i, 1]]^2 - 
               2*mags[[i, 
                 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Sin[tht[t]]*
                   Cos[phi[t]] + 
                  Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
                   Sin[phi[t]] + 
                  Cos[mags[[i, 2]]]*Cos[tht[t]]))^1.5, {i, 
           4}]/(((M/3 + m)*L^2)*Sin[tht[t]]^2) - 
       2*tht'[t]*phi'[t]/Tan[tht[t]] - 
       k*(L^2)*phi'[t]/(3*((M/3 + m)*L^2))) == 0, tht'[0] == 0, 
   phi'[0] == 0, tht[0] == tht0, phi[0] == phi0};

getStateData[op_] := 
  First@NDSolve`ProcessEquations[eq /. op, {tht, phi}, t, 
    Method -> "Adams"];
op = {k -> .1, c -> .0004, M -> .05, m -> .01, L -> .3, g -> 9.8, 
   tht0 -> 1, phi0 -> 3};
sd = getStateData[op]




getBaSin[x_, y_] := 
 Module[{state = sd, sol}, 
  state = First[NDSolve`Reinitialize[sd, {tht[0] == x, phi[0] == y}]];
  NDSolve`Iterate[state, 21];
  sol = {tht[20], phi[20]} /. NDSolve`ProcessSolutions[state];
  nf[{0.3*Sin[sol[[1]]]*Cos[sol[[2]]], 
    0.3*Sin[sol[[1]]]*Sin[sol[[2]]]}]]

ArrayPlot[
  ParallelTable[
   getBaSin[x, y], {x, 0.1, 1.5, 0.01}, {y, 0.1, 6.2, 0.025}], 
  ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue}] // AbsoluteTiming

fig1

$\endgroup$
  • $\begingroup$ Thank you so much! May I ask you a few more questions? First, what is the principle of this code solving differential equations numerically? I mean I thought there should be somewhere defining the time step, which is not shown in my code. Second, does NDSolve`Iterate[state, 21] means that function works about 20 times? Is it related to the time step? Lastly, is Adams' method usually better than RK4 method? $\endgroup$ – 최준호 May 10 at 16:44
  • $\begingroup$ Sorry, I accidentally edited my question... I hope you are able to see my edited comment $\endgroup$ – 최준호 May 10 at 17:19
  • $\begingroup$ 1) Read the tutorial section Advanced Numerical Differential Equation Solving in the Wolfram Language -> Components and Data Structures. 2) 21 is a time of integration (t,0,21}, it is related to the time step rough as 21/number of steps. 3) Adams method is 20% faster than Automatic whereas with the Method -> "ExplicitRungeKutta" you can not find a solution. $\endgroup$ – Alex Trounev May 10 at 17:26

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.