2
$\begingroup$

I am trying to solve a problem, for what I finally found a book what has some Mathematica files supplied, but I am stuck now as I cannot run the files.

My problem is that I cannot run the program as it is written in Mathematica 3.0 and I don't know what should I change to make it run under today's Mathematica versions. Here is the error it returns.

FindMinimum::fmgz: Encountered a gradient that is effectively zero. The result returned may not be a minimum; it may be a maximum or a saddle point.

Here is the original code copy and pasted (called MOTIPOIN.NB in the original zip):

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]
$\endgroup$
1
  • $\begingroup$ What you mentioned is a warning and not an error. $\endgroup$
    – Listing
    Commented Jul 1, 2011 at 6:42

2 Answers 2

2
$\begingroup$

Goto[] and Label[]? How baroque quaint! I don't have Mathematica with me at the moment, but could you first check that NDSolve[] is indeed outputting a nonzero InterpolatingFunction[]? Otherwise, the syntax for FindMinimum[] ought to be correct (though FindMinimum[-z, {t, 0, 0, b trif}] probably ought to be FindMaximum[z, {t, 0, 0, b trif}]).

Probably a better bet to rewrite the whole thing from scratch; it's a bleeding mess, it looks.

$\endgroup$
1
  • $\begingroup$ Thanks Jerry! I never used NDSolve and InterpolatingFunction and stuff like this in Mathematica, thus I really don't know what function should output what value. Can you have a look at it later? $\endgroup$
    – hyperknot
    Commented Jul 1, 2011 at 15:45
0
$\begingroup$

Here's my meagre attempt at improving the code:

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];]

I do know for a fact that the herpolhode differential equations can be solved in terms of elliptic functions, but I don't have the time to do that derivation now...

$\endgroup$
1

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .