Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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]
share|improve this question
    
What you mentioned is a warning and not an error. –  Listing Jul 1 '11 at 6:42

2 Answers 2

up vote 0 down vote accepted

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...

share|improve this answer
    
Sorry, it tells me an error that "Encountered a gradient that is effectively zero." I've got help here: groups.google.com/d/topic/comp.soft-sys.math.mathematica/… that I should use FindMinimum[z, {t, 0.01, 0, b trif}]. Maybe this is the problem? –  zsero Jul 12 '11 at 10:58

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.

share|improve this answer
    
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? –  zsero Jul 1 '11 at 15:45

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.