ayuda con un archivo antiguo de Mathematica

Question

Estoy intentando resolver un problema Finalmente encontré un libro que contiene algunos archivos de Mathematica, pero ahora estoy atascado porque no puedo ejecutar los archivos.

Mi problema es que no puedo ejecutar el programa tal y como está escrito en Mathematica 3.0 y no sé qué debería cambiar para hacerlo funcionar bajo las versiones actuales de Mathematica. Este es el error que devuelve.

BuscarMínimo::fmgz: Encountered a gradiente que es efectivamente cero. El resultado de resultado devuelto puede no ser un mínimo; puede ser un máximo o un punto de silla.

Aquí está el código original copiado y pegado (llamado MOTIPOIN.NB en el archivo zip original ):

Off[General::"spell"]
Off[General::"spell1"]

MotiPoin[A_, B_, C0_, r0_, theta0_, b_, alpha_] :=
Module[{q0, trif, K2, T, h, eq},
q0 = (C0 r0 Tan[theta0])/B;
trif = (2 \[Pi] B)/(r0 (A Cos[theta0] + C0 Sin[theta0] Tan[theta0]));
K2 = (B q0)^2 + (C0 r0)^2;
T = 1/2 (B q0^2 + C0 r0^2);
h = Sqrt[(2 T)/K2];
eq1 = Derivative[1][p][t] == ((B - C0) q[t] r[t])/A;
eq2 = Derivative[1][q][t] == ((C0 - A) p[t] r[t])/B;
eq3 = Derivative[1][r][t] == ((A - B) p[t] q[t])/C0;
eq4 = Derivative[1][psi][t] == (Cos[phi[t]] q[t] + p[t] Sin[phi[t]])/Sin[theta[t]];
eq5 = Derivative[1][phi][t] == r[t] - Cot[theta[t]] (Cos[phi[t]] q[t] + p[t] Sin[phi[t]]);
eq6 = Derivative[1][theta][t] == p[t] Cos[phi[t]] - q[t] Sin[phi[t]];
w1 = (Cos[phi[t]] Cos[psi[t]] - Sin[phi[t]] Sin[psi[t]] Cos[theta[t]]) p[t] - (Cos[psi[t]] Sin[phi[t]] + Cos[phi[t]] Cos[theta[t]] Sin[psi[t]]) q[t] + r[t] Sin[psi[t]] Sin[theta[t]];
w2 = (Cos[psi[t]] Cos[theta[t]] Sin[phi[t]] + Cos[phi[t]] Sin[psi[t]]) p[t] + (Cos[phi[t]] Cos[psi[t]] Cos[theta[t]] - Sin[phi[t]] Sin[psi[t]]) q[t] - Cos[psi[t]] r[t] Sin[theta[t]];
w3 = Cos[theta[t]] r[t] + Cos[phi[t]] q[t] Sin[theta[t]] + p[t] Sin[phi[t]] Sin[theta[t]];
sol = NDSolve[{eq1, eq2, eq3, eq4, eq5, eq6, p[0] == 0, q[0] == q0, r[0] == r0, psi[0] == 0, phi[0] == 0, theta[0] == theta0}, {p, q, r, psi, phi, theta}, {t, 0, b trif}];
{x, y} = Flatten[{-((w1 h)/w3), -((w2 h)/w3)} /. sol];
z = x^2 + y^2;
If[A < C0 < B || B < C0 < A, Goto[2], Goto[1]];

Label[1];
m = FindMinimum[z, {t, 0, 0, b trif}];
M = FindMinimum[-z, {t, 0, 0, b trif}];
ra1 = Sqrt[m[[1]]];
ra2 = Sqrt[-M[[1]]];
Print["L'erpoloide è contenuta in una corona circolare"];
Print["avente raggio interno ra1 e raggio esterno ra2"];
Print["ra1=", ra1]; Print["ra2=", ra2];
c1 = ParametricPlot[{ra1 Sin[u], ra1 Cos[u]}, {u, 0, 2 \[Pi]}, AspectRatio -> 1, DisplayFunction -> Identity, PlotStyle -> RGBColor[0.8669, 0.258, 0.227]];
c2 = ParametricPlot[{ra2 Sin[u], ra2 Cos[u]}, {u, 0, 2 \[Pi]}, AspectRatio -> 1, DisplayFunction -> Identity, PlotStyle -> RGBColor[0.925, 0.140, 0.129]];
Plot[Sqrt[z], {t, 0, b trif}, AxesLabel -> {"t", "ra"}];
erp = ParametricPlot[{x, y}, {t, 0, b trif}, AspectRatio -> 1, PlotRange -> All, DisplayFunction -> Identity];
Show[erp, c1, c2, DisplayFunction -> $DisplayFunction];
Goto[3];

Label[2];
Plot[Sqrt[z], {t, 0, b trif}, AxesLabel -> {"t", "ra"}];
erp = ParametricPlot[{x, y}, {t, 0, b trif}, AspectRatio -> 1, PlotRange -> All];

Label[3];
xp = p[t]/Sqrt[2 T] /. sol;
yp = q[t]/Sqrt[2 T] /. sol;
zp = r[t]/Sqrt[2 T] /. sol;
X = (Cos[u] Sin[v])/Sqrt[A];
Y = (Sin[u] Sin[v])/Sqrt[B];
Z = Cos[v]/Sqrt[C0];
el = ParametricPlot3D[{X, Y, Z}, {u, 0, 2 \[Pi]}, {v, 0, alpha}, LightSources -> {{{-1, -1, 3}, GrayLevel[0.999]}}, Boxed -> False, DisplayFunction -> Identity];
pol = ParametricPlot3D[Evaluate[Flatten[{xp, yp, zp}] /. sol], {t, 0, b trif}, PlotPoints -> 200, DisplayFunction -> Identity];
Show[el, pol, DisplayFunction -> $DisplayFunction];
]

MotiPoin[1,1.5,0.5,3,Pi/4,1.5,Pi/4]
MotiPoin[1, 1.5, 0.5, 3, 0.01, 1.5, Pi/100]
MotiPoin[0.5, 1.5, 1, -3, 0.01, 3.5, Pi]
MotiPoin[1, 1, 1.5, 3, Pi/4, 2.5, Pi/2]

Answer 1

Goto[] y Label[] ? ¿Cómo barroco ¡Qué pintoresco! No tengo Mathematica conmigo en este momento, pero ¿podría comprobar primero que NDSolve[] emite un valor distinto de cero InterpolatingFunction[] ? De lo contrario, la sintaxis para FindMinimum[] debería ser correcta (aunque FindMinimum[-z, {t, 0, 0, b trif}] probablemente debería ser FindMaximum[z, {t, 0, 0, b trif}] ).

Probablemente sea mejor reescribirlo todo desde cero; parece un lío sangrante.

Answer 2

Aquí está mi magro intento de mejorar el código:

MotiPoin[A_, B_, C0_, r0_, theta0_, b_, alpha_] := 
  Module[{q0, trif, K2, T, h, p, q, r, psi, phi, theta},
    q0 = (C0 r0 Tan[theta0])/B; 
    trif = (2 Pi B)/(r0 (A Cos[theta0] + C0 Sin[theta0] Tan[theta0]));
    K2 = (B q0)^2 + (C0 r0)^2; T = (B q0^2 + C0 r0^2)/2;
    h = Sqrt[2 T/K2];
    {p, q, r, psi, phi, theta} = 
      First[{p, q, r, psi, phi, theta} /. 
          NDSolve[{p'[t] == ((B - C0) q[t] r[t])/A, 
              q'[t] == ((C0 - A) p[t] r[t])/B, 
              r'[t] == ((A - B) p[t] q[t])/C0, 
              psi'[t] == (Cos[phi[t]] q[t] + p[t] Sin[phi[t]])/Sin[theta[t]], 
              phi'[t] == 
                r[t] - Cot[theta[t]] (Cos[phi[t]] q[t] + p[t] Sin[phi[t]]),
              theta'[t] == p[t] Cos[phi[t]] - q[t] Sin[phi[t]], p[0] == 0, 
              q[0] == q0, r[0] == r0, psi[0] == 0, phi[0] == 0, 
              theta[0] == theta0}, {p, q, r, psi, phi, theta}, {t, 0, 
              b trif}]];
    {x, y} = -h{(Cos[phi[t]] Cos[psi[t]] - 
                    Sin[phi[t]] Sin[psi[t]] Cos[theta[t]]) p[
                  t] - (Cos[psi[t]] Sin[phi[t]] + 
                    Cos[phi[t]] Cos[theta[t]] Sin[psi[t]]) q[t] + 
              r[t] Sin[psi[t]] Sin[
                  theta[t]], (Cos[psi[t]] Cos[theta[t]] Sin[phi[t]] + 
                    Cos[phi[t]] Sin[psi[t]]) p[
                  t] + (Cos[phi[t]] Cos[psi[t]] Cos[theta[t]] - 
                    Sin[phi[t]] Sin[psi[t]]) q[t] - 
              Cos[psi[t]] r[t] Sin[theta[t]]}/(Cos[theta[t]] r[t] + 
              Cos[phi[t]] q[t] Sin[theta[t]] + 
              p[t] Sin[phi[t]] Sin[theta[t]]);
    z = x^2 + y^2;
    Plot[Sqrt[z], {t, 0, b trif}, AxesLabel -> {"t", "ra"}];
    If[A < C0 < B || B < C0 < A,
      ParametricPlot[{x, y}, {t, 0, b trif}, AspectRatio -> 1, 
          PlotRange -> All];,
      ra1 = Sqrt[First[FindMinimum[z, {t, 0, 0, b trif}]]];
      ra2 = Sqrt[First[FindMaximum[z, {t, 0, 0, b trif}]]];
      Print[
        StringForm[
          "The herpolhode is contained in an annulus with inner radius `` and \
outer radius ``.", ra1, ra2]];
      ParametricPlot[{x, y}, {t, 0, b trif}, AspectRatio -> Automatic, 
        Epilog -> {{RGBColor[0.8669, 0.258, 0.227], 
              Circle[{0, 0}, ra1]}, {RGBColor[0.925, 0.140, 0.129], 
              Circle[{0, 0}, ra2]}}, PlotRange -> All];];
    el = ParametricPlot3D[{(Cos[u] Sin[v])/Sqrt[A], (Sin[u] Sin[v])/Sqrt[B], 
          Cos[v]/Sqrt[C0], SurfaceColor[GrayLevel[.75]]}, {u, 0, 2 Pi}, {v, 0,
           alpha}, AmbientLight -> GrayLevel[1], Boxed -> False, 
        DisplayFunction -> Identity, LightSources -> {}];
    pol = 
      ParametricPlot3D[
        Evaluate[
          Append[{p[t], q[t], r[t]}/Sqrt[2 T], {AbsoluteThickness[3], 
              RGBColor[0, 0, 1]}]], {t, 0, b trif}, PlotPoints -> 200, 
        DisplayFunction -> Identity];
    Show[el, pol, DisplayFunction -> $DisplayFunction];]

Sé a ciencia cierta que las ecuaciones diferenciales herpolhode pueden resolverse en términos de funciones elípticas, pero ahora no tengo tiempo para hacer esa derivación...