Mathematica NDSolve Running Infinitely

  • Context: Mathematica 
  • Thread starter Thread starter AtLastForgot
  • Start date Start date
  • Tags Tags
    Mathematica Running
Click For Summary
SUMMARY

The discussion focuses on troubleshooting an infinite loop issue encountered while using Mathematica's NDSolve function to solve a system of equations related to coordinate transformations and Lagrangian mechanics. The user provided a detailed code snippet that defines transformations and calculates kinetic energy and potential energy. Key recommendations include simplifying the code by stripping down to essential components, using explicit constant values, and gradually rebuilding the code to isolate errors. The user successfully managed to create a working version after following these guidelines.

PREREQUISITES
  • Familiarity with Mathematica 12.0 syntax and functions
  • Understanding of Lagrangian mechanics and its application in physics
  • Knowledge of coordinate transformations in mathematical modeling
  • Experience with debugging and optimizing code in computational environments
NEXT STEPS
  • Explore Mathematica's NDSolve documentation for advanced usage and troubleshooting techniques
  • Learn about Lagrangian mechanics and its mathematical formulations
  • Investigate methods for simplifying complex systems in Mathematica
  • Study debugging strategies in Mathematica to effectively identify and resolve errors
USEFUL FOR

Students in physics or engineering, researchers using Mathematica for computational modeling, and anyone involved in solving complex differential equations in a programming context.

AtLastForgot
Messages
2
Reaction score
0
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 \
convenience.*)
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 =
unhat[Inverse[
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]],
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 \
simpler.*)
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]] +
m3*grav*(gWeapon[x[t],
y[t], \[Theta]1[t], \[Theta]2[t], \[Theta]3[t]].{0, 0, 1})[[
2]];
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]}];
 
Physics news on Phys.org
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.
 
Thanks. I think I managed to get a working version going, but now I'm getting different errors. Hurray for computer programming!
 

Similar threads

  • · Replies 1 ·
Replies
1
Views
3K
  • · Replies 1 ·
Replies
1
Views
3K
  • · Replies 1 ·
Replies
1
Views
3K
Replies
5
Views
2K
  • · Replies 6 ·
Replies
6
Views
7K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 4 ·
Replies
4
Views
2K
  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 8 ·
Replies
8
Views
2K
  • · Replies 4 ·
Replies
4
Views
3K