Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica NDSolve Running Infinitely

  1. Dec 8, 2011 #1
    I'm trying to use Mathematica to solve a system of equations based on coordinate transformations and Lagrangian mechanics. The code follows here. However, when the system of equations is run, Mathematica never reaches an answer nor does it give an error, it simply runs infinitely. This is for a class final project that's worth a large portion of my grade.

    (*Definition of various functions which will help us simplify our \
    code later. None of this is physically relevant; it is simply a \
    g[x_, y_, \[Theta]_] := {{Cos[\[Theta]], -Sin[\[Theta]],
    x}, {Sin[\[Theta]], Cos[\[Theta]], y}, {0, 0, 1}}
    unhat[a_] := {a[[1, 3]], a[[2, 3]], a[[2, 1]]}
    hat[a_] := {{0, -a[[3]], a[[1]]}, {a[[3]], 0, a[[2]]}, {0, 0, 0}}
    Hom[point_] := Flatten[{point, 1}]
    Unhom[point_] := {point[[1]], point[[2]]}
    Clear[m1, m2, m3, J1, J2, J3, L1, L2, L3, grav]

    (*The first transformation from an origin of 0,0 to the primary axes \
    of the first bar.*)
    gUpperArm[x_, y_, \[Theta]1_] := g[x, y, 0].g[0, 0, \[Theta]1];

    (*Velocity of this point is determined both in terms of translational \
    velocities as well as rotational. Unhat allows us to calculate this \
    quantity, normally [g^-1][g'], via simple notation.*)
    VUpperArm =
    unhat[Inverse[gUpperArm[x[t], y[t], \[Theta]1[t]]].D[
    gUpperArm[x[t], y[t], \[Theta]1[t]], t]];

    (*Transformation from origin to the primary axes of the second bar \
    using intermediate steps of the first transformation. This links \
    their two motions and allows us to calculate VForearm, a quantity \
    which would otherwise be extremely difficult to describe.*)
    gForearm[x_, y_, \[Theta]1_, \[Theta]2_] :=
    gUpperArm[x, y, \[Theta]1].g[L1/2, 0, 0].g[0, 0, \[Theta]2].g[L2/2,
    0, 0];
    VForearm =
    unhat[Inverse[gForearm[x[t], y[t], \[Theta]1[t], \[Theta]2[t]]] .
    D[gForearm[x[t], y[t], \[Theta]1[t], \[Theta]2[t]], t]];
    gWeapon[x_, y_, \[Theta]1_, \[Theta]2_, \[Theta]3_] :=
    gForearm[x, y, \[Theta]1, \[Theta]2].g[L2/2, 0, 0].g[0,
    0, \[Theta]3].g[L3/4, 0, 0];
    VWeapon =
    gWeapon[x[t], y[t], \[Theta]1[t], \[Theta]2[t], \[Theta]3[t]]] .
    D[gWeapon[x[t], y[t], \[Theta]1[t], \[Theta]2[t], \[Theta]3[t]],

    (*Definitions of inertia for the various components of the system. \
    Notice that the form is x,y,\[Theta].*)
    IUpperArm = DiagonalMatrix[{m1, m1, J1}];
    IForearm = DiagonalMatrix[{m2, m2, J2}];
    IWeapon = DiagonalMatrix[{m3, m3, J3}];

    (*Calculations of Kinetic Energy using previously defined \
    transformations to make this otherwise nearly impossible task much \
    KE = ((1/2)*VUpperArm.IUpperArm.VUpperArm) + ((1/2)*
    VForearm.IForearm.VForearm) + ((1/2)*VWeapon.IWeapon.VWeapon);
    V = m1*grav*(gUpperArm[x[t], y[t], \[Theta]1[t]].{0, 0, 1})[[2]] +
    m2*grav*(gForearm[x[t], y[t], \[Theta]1[t], \[Theta]2[t]].{0, 0,
    1})[[2]] +
    y[t], \[Theta]1[t], \[Theta]2[t], \[Theta]3[t]].{0, 0, 1})[[
    Lag = KE - V;

    EQ1 = D[D[Lag, x'[t]], t] - D[Lag, x[t]];
    EQ2 = D[D[Lag, y'[t]], t] - D[Lag, x[t]];
    EQ3 = D[D[Lag, \[Theta]1'[t]], t] - D[Lag, \[Theta]1[t]];
    EQ4 = D[D[Lag, \[Theta]2'[t]], t] - D[Lag, \[Theta]2[t]];
    EQ5 = D[D[Lag, \[Theta]3'[t]], t] - D[Lag, \[Theta]3[t]];
    EQ = Solve[{EQ1 == 0, EQ2 == 0, EQ2 == 0, EQ4 == 0, EQ5 == 0}, {x[t],
    y[t], \[Theta]1[t], \[Theta]2[t], \[Theta]3[t]}];
  2. jcsd
  3. Dec 8, 2011 #2
    I can give you a general guideline which I think works well: keep stripping down your code until it works, run a simpler version of your code, use explicit values of constants rather than assigned values, cut down the actual differential equations, drop some of the terms, make the terms simpler. Then once you get it working with just a bare-bones version even if it's significantly different than what you want, start building it back up piece by piece and in that way, a small piece by piece, you'll have a better handle on a problem when it occurs because you'll know specifically what you added last that caused the problem and can then focus on that addition and hopefully resolve it.
  4. Dec 8, 2011 #3
    Thanks. I think I managed to get a working version going, but now I'm getting different errors. Hurray for computer programming!
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Similar Threads for Mathematica NDSolve Running
Mathematica Storing Mathematica output
Mathematica Mathematica to MATLAB
Mathematica Cannot do the integral of the Hyper-geometric function?
Mathematica Cannot Plot This Function in Mathematica
Mathematica While Loop in Mathematica