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.
