Manipulate[
(*Nasser M. Abbasi, version 8/11/2014*)
tick;
Module[{a1, a2, a3, v1, v2, v3, h1, a1d, a2d, a3d, z},
Which[state == "init" || state2 == "on",
state2 = "off";
If[state == "init", state = "paused"];
{pe, ke, sol} =
solve[initial\[Theta]1*Pi/180., initial\[Theta]2*Pi/180., initial\[Theta]3*Pi/180., initial\[Theta]1dot, initial\[Theta]2dot,
initial\[Theta]3dot, 0, maxTime, \[Theta]1, \[Theta]2, \[Theta]3, t, m1, m2, m3, L, k1, k2, len];
currentTime = 0,
state == "running" || state == "step",
If[currentTime > maxTime, currentTime = 0]
];
a1 = (\[Theta]1[t] /. sol) /. t -> currentTime; a2 = (\[Theta]2[t] /. sol) /. t -> currentTime;
a3 = (\[Theta]3[t] /. sol) /. t -> currentTime; v1 = (\[Theta]1'[t] /. sol) /. t -> currentTime;
v2 = (\[Theta]2'[t] /. sol) /. t -> currentTime; v3 = (\[Theta]3'[t] /. sol) /. t -> currentTime;
If[state == "running" || state == "step",
currentTime = currentTime + timeStep;
If[state == "running", tick += del]
];
ppe = pe /. {\[Theta]1[t] -> a1, \[Theta]2[t] -> a2, \[Theta]3[t] -> a3};
kke = ke /. {\[Theta]1'[t] -> v1, \[Theta]2'[t] -> v2, \[Theta]3'[t] -> v3};
totalE = kke + Abs[ppe];
kkef = kke/totalE;
ppef = Abs[ppe]/totalE;
{a1d, a2d, a3d} = (z = Mod[#*180/Pi, 360]; Round@If[z > 180, z - 360, z]) & /@ {a1, a2, a3};(*normalize*)
h1 = Style[Grid[{
{"time (sec)",
"P.E. (J)",
"K.E. (J)",
"energy (J)",
Row[{Subscript["\[Theta]", 1], " (deg)"}],
Row[{Subscript["\[Theta]", 2], " (deg)"}],
Row[{Subscript["\[Theta]", 3], " (deg)"}]
},
{padIt2[currentTime, {5, 3}],
padIt1[ppe, {6, 3}],
padIt2[kke, {6, 3}],
padIt1[totalE, {6, 3}],
padIt1[a1d, 3],
padIt1[a2d, 3],
padIt1[a3d, 3]
}
},
Frame -> All, FrameStyle -> Directive[Thickness[.001], Gray], Spacings -> {1, 1.2},
Alignment -> Center]
, 11];
Text@Grid[{
{h1},
{
Framed[Evaluate@Graphics[
{
{EdgeForm[Thin], Opacity[.8], Gray, Rectangle[{-1.5 w, -w/2}, {2 len + 1.5 w, w/2}, RoundingRadius -> 0]},
Rotate[{
{EdgeForm[Thin], Opacity[m1], Red, Rectangle[{-w/2, 0}, {w/2, -L}, RoundingRadius -> 0]},
{EdgeForm[Thin], Blue, Opacity[.3], Disk[{0, 0}, w/2]},
{Thick, Disk[{0, 0}, w/5]},
{Thick, Disk[{0, -L}, w/5]}
}, a1, {0, 0}],
Rotate[{
{EdgeForm[Thin], Opacity[m2], Red, Rectangle[{len - w/2, 0}, {len + w/2, -L}, RoundingRadius -> 0]},
{EdgeForm[Thin], Blue, Opacity[.3], Disk[{len, 0}, w/2]},
{Thick, Disk[{len, 0}, w/5]},
{Thick, Disk[{len, -L}, w/5]}
}, a2, {len, 0}],
Rotate[{
{EdgeForm[Thin], Opacity[m3], Red, Rectangle[{2 len - w/2, 0}, {2 len + w/2, -L}, RoundingRadius -> 0]},
{EdgeForm[Thin], Blue, Opacity[.3], Disk[{2 len, 0}, w/2]},
{Thick, Disk[{2 len, 0}, w/5]},
{Thick, Disk[{2 len, -L}, w/5]}
}, a3, {2 len, 0}],
If[k1 > 0,
{AbsoluteThickness[k1/5], Line[makeSpring[L Sin[a1], -L Cos[a1], len + L Sin[a2], -L Cos[a2], .15]]}
],
If[k2 > 0,
{AbsoluteThickness[k2/5], Line[makeSpring[len + L Sin[a2], -L Cos[a2], 2 len + L Sin[a3], -L Cos[a3], .15]]}
]
},
PlotRange -> {{-.8 L, 2 len + 0.8 L}, {-L, L}}, Frame -> False, Axes -> False, ImageSize -> {450}, ImagePadding -> 15,
Background -> RGBColor[.99, .97, .90]
], FrameStyle -> LightGray]
}}, Spacings -> {.1, .1}]
],
Text@Grid[{
{
Grid[{
{
Button[Text[Style["run", 11]], state = "running"; tick += del, ImageSize -> {50, 35}],
Button[Text[Style["pause", 11]], state = "paused"; tick += del, ImageSize -> {50, 35}]
},
{Button[Text[Style["step", 11]], state = "step"; tick += del, ImageSize -> {50, 35}],
Button[Text[Style["reset", 11]], state = "init"; tick = 0, ImageSize -> {50, 35}
]
}
}, Spacings -> {0.5, .5}, Alignment -> Center
]
,
Grid[{
{
"slow", Manipulator[Dynamic[timeStep, {timeStep = #;} &], {0.001, 0.1, .001}, ImageSize -> Tiny], "fast"
},
{
TraditionalForm@Style[Subscript[k, 1], 11],
Manipulator[Dynamic[k1, {k1 = #; tick += del; state2 = "on"} &], {0, 10, .1}, ImageSize -> Tiny], Dynamic[padIt2[k1, {3, 1}]]
},
{
TraditionalForm@Style[Subscript[k, 2], 11],
Manipulator[Dynamic[k2, {k2 = #; tick += del; state2 = "on"} &], {0, 10, .1}, ImageSize -> Tiny], Dynamic[padIt2[k2, {3, 1}]]
},
{
TraditionalForm@Style[Subscript[m, 1], 11],
Manipulator[Dynamic[m1, {m1 = #; tick += del; state2 = "on"} &], {0.1, 1, .1}, ImageSize -> Tiny],
Dynamic[padIt2[m1, {2, 1}]]
},
{
TraditionalForm@Style[Subscript[m, 2], 11],
Manipulator[Dynamic[m2, {m2 = #; tick += del; state2 = "on"} &], {0.1, 1, .1}, ImageSize -> Tiny],
Dynamic[padIt2[m2, {2, 1}]]
},
{
TraditionalForm@Style[Subscript[m, 3], 11],
Manipulator[Dynamic[m3, {m3 = #; tick += del; state2 = "on"} &], {0.1, 1, .1}, ImageSize -> Tiny],
Dynamic[padIt2[m3, {2, 1}]]
}
}, Frame -> True, Alignment -> Center, Spacings -> {.5, .5}, FrameStyle -> Directive[Thickness[.001], Gray]
]
,
Grid[{
{
Row[{Subscript["\[Theta]", 1], "(0)"}],
Manipulator[Dynamic[initial\[Theta]1, {initial\[Theta]1 = #; tick += del; state2 = "on"} &], {-90, 90, 1},
ImageSize -> Tiny], Row[{Dynamic[padIt1[initial\[Theta]1, 2]], Spacer[1], Degree}]
},
{
Row[{Subscript["\[Theta]", 2], "(0)"}],
Manipulator[Dynamic[initial\[Theta]2, {initial\[Theta]2 = #; tick += del; state2 = "on"} &], {-90, 90, 1},
ImageSize -> Tiny], Row[{Dynamic[padIt1[initial\[Theta]2, 2]], Spacer[1], Degree}]
},
{
Row[{Subscript["\[Theta]", 3], "(0)"}],
Manipulator[Dynamic[initial\[Theta]3, {initial\[Theta]3 = #; tick += del; state2 = "on"} &], {-90, 90, 1},
ImageSize -> Tiny], Row[{Dynamic[padIt1[initial\[Theta]1, 2]], Spacer[1], Degree}]
},
{
Row[{Subscript["\[Theta]", 1]', "(0)"}],
Manipulator[Dynamic[initial\[Theta]1dot, {initial\[Theta]1dot = #; tick += del; state2 = "on"} &], {-2, 2, .1},
ImageSize -> Tiny], Dynamic[padIt1[initial\[Theta]1dot, {2, 1}]]
},
{
Row[{Subscript["\[Theta]", 2]', "(0)"}],
Manipulator[Dynamic[initial\[Theta]2dot, {initial\[Theta]2dot = #; tick += del; state2 = "on"} &], {-2, 2, .1},
ImageSize -> Tiny], Dynamic[padIt1[initial\[Theta]2dot, {2, 1}]]
},
{
Row[{Subscript["\[Theta]", 3]', "(0)"}],
Manipulator[Dynamic[initial\[Theta]3dot, {initial\[Theta]3dot = #; tick += del; state2 = "on"} &], {-2, 2, .1},
ImageSize -> Tiny], Dynamic[padIt1[initial\[Theta]3dot, {2, 1}]]
}
}, Frame -> True, Alignment -> Center, Spacings -> {.5, .5}, FrameStyle -> Directive[Thickness[.001], Gray]
],
Grid[{
{Evaluate@Graphics[
{
{Text[Style["K.E.", 11], {w/2, .1}]},
{Text[Style["P.E.", 11], {2.5 w, .1}]},
{EdgeForm[Thin], Green, Rectangle[{-0.4 w, -1}, {1.4 w, -1 + kkef}, RoundingRadius -> 0]},
{EdgeForm[Thin], Blue, Rectangle[{1.6 w, -1}, {3.4 w, -1 + ppef}, RoundingRadius -> 0]}
}
,
PlotRange -> {{-w/2, 3.5 w}, {-1, 0.1}}, Frame -> False, Axes -> False, ImageSize -> {120}, ImagePadding -> 5,
ImageMargins -> 1
]
}}, Frame -> True, FrameStyle -> Directive[Thickness[.001], Gray]]
}
}, Alignment -> Center, Spacings -> {.5, .5}],
{{del, $MachineEpsilon}, None},
{{tick, 0}, None},
{{maxTime, 20}, None},
{{w, 0.2}, None},
{{totalE, 0}, None},
{{ppe, 0}, None},
{{ppef, 0}, None},
{{kke, 0}, None},
{{kkef, 0}, None},
{{k1, 2.1}, None},
{{k2, 5}, None},
{{m1, 1}, None},
{{m2, 0.5}, None},
{{m3, 1}, None},
{{state, "init"}, None},
{{state2, "off"}, None},
{{currentTime, 0}, None},
{{sol, {}}, None},
{{len, 2.2}, None},
{{L, 1}, None},
{{initial\[Theta]1, 60}, None},
{{initial\[Theta]2, 17}, None},
{{initial\[Theta]3, -10}, None},
{{initial\[Theta]1dot, 1.4}, None},
{{initial\[Theta]2dot, 1.6}, None},
{{initial\[Theta]3dot, -1.3}, None},
{{pe, 0}, None},
{{ke, 0}, None},
{{timeStep, 0.02}, None},
ControlPlacement -> Above,
SynchronousUpdating -> False,
SynchronousInitialization -> False,
ContinuousAction -> False,
Alignment -> Center,
ImageMargins -> 5,
FrameMargins -> 5,
Paneled -> True,
Frame -> False,
AutorunSequencing -> {1},
TrackedSymbols :> {tick},
Initialization :>
(
integerStrictPositive = (IntegerQ[#] && # > 0 &);
integerPositive = (IntegerQ[#] && # >= 0 &);
numericStrictPositive = (Element[#, Reals] && # > 0 &);
numericPositive = (Element[#, Reals] && # >= 0 &);
numericStrictNegative = (Element[#, Reals] && # < 0 &);
numericNegative = (Element[#, Reals] && # <= 0 &);
bool = (Element[#, Booleans] &);
numeric = (Element[#, Reals] &);
integer = (Element[#, Integers] &);
(*--------------------------------------------*)
(* helper function for formatting *)
(*--------------------------------------------*)
padIt1[v_?numeric, f_List] :=
AccountingForm[Chop[v] , f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"}, SignPadding -> True];
padIt1[v_?numeric, f_Integer] :=
AccountingForm[v , f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"}, SignPadding -> True];
(*--------------------------------------------*)
(* helper function for formatting *)
(*--------------------------------------------*)
padIt2[v_?numeric, f_List] :=
AccountingForm[Chop[v] , f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"}, SignPadding -> True];
padIt2[v_?numeric, f_Integer] := AccountingForm[v , f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"}, SignPadding -> True];
(*-------------------------------*)
makeSpring[xFirst_?numeric,
yFirst_?numeric,
xEnd_?numeric,
yEnd_?numeric,
szel_?numeric(*the larger, the wider the spring twists*)] := Module[{hx, veghossz, hossz, hy, dh, tbl},
hx = xEnd - xFirst;
If[Abs[hx] <= $MachineEpsilon, hx = 10^-6];
hy = yEnd - yFirst;
If[Abs[hy] <= $MachineEpsilon, hy = 10^-6];
veghossz = 0.28;
hossz = Sqrt[hx^2 + hy^2];
dh = (hossz - .5*veghossz)/30.;
tbl = Table[If[
OddQ[i], {xFirst + hx*(i*dh + veghossz)/hossz + hy*szel/hossz,
yFirst + hy*(i*dh + veghossz)/hossz - hx*szel/hossz}, {xFirst + hx*(i*dh + veghossz)/hossz - hy*szel/hossz,
yFirst + hy*(i*dh + veghossz)/hossz + hx*szel/hossz}], {i, 2, 18}];
{{xFirst, yFirst}}~Join~{{xFirst + hx*(dh + veghossz)/hossz, yFirst + hy*(dh + veghossz)/hossz}}~Join~tbl~
Join~{{xFirst + hx*(19*dh + veghossz)/hossz, yFirst + hy*(19*dh + veghossz)/hossz}}~Join~{{xEnd, yEnd}}
];
(*-------------------------------*)
solve[initial\[Theta]1_, initial\[Theta]2_, initial\[Theta]3_, initial\[Theta]1dot_, initial\[Theta]2dot_, initial\[Theta]3dot_,
from_, to_, \[Theta]1_, \[Theta]2_, \[Theta]3_, t_, m1_, m2_, m3_, L_, k1_, k2_, len_] :=
Module[{I1, I2, I3, spring1, spring2, pe, ke, lag, eq1, eq2, eq3, sol, ic, g = 9.8},
I1 = 1/3 m1 L^2;
I2 = 1/3 m2 L^2;
I3 = 1/3 m3 L^2;
ke = 1/2 I1 \[Theta]1'[t]^2 + 1/2 I2 \[Theta]2'[t]^2 + 1/2 I3 \[Theta]3'[t]^2;
spring1 = Sqrt[(len - L Sin[\[Theta]1[t]])^2 + L^2 Sin[\[Theta]2[t]]^2 + 2 L Sin[\[Theta]2[t]] (len - L Sin[\[Theta]1[t]]) +
L^2 + L^2 Cos[\[Theta]2[t]]^2 + L^2 (1 - Cos[\[Theta]1[t]])^2 + 2 L^2 Cos[\[Theta]2[t]] (1 - Cos[\[Theta]1[t]]) -
2 L^2 (Cos[\[Theta]2[t]] + 1 - Cos[\[Theta]1[t]])] - len;(*spring 1 extension*)
spring2 = Sqrt[(len - L Sin[\[Theta]2[t]])^2 + L^2 Sin[\[Theta]3[t]]^2 + 2 L Sin[\[Theta]3[t]] (len - L Sin[\[Theta]2[t]]) +
L^2 + L^2 Cos[\[Theta]3[t]]^2 + L^2 (1 - Cos[\[Theta]2[t]])^2 + 2 L^2 Cos[\[Theta]3[t]] (1 - Cos[\[Theta]2[t]]) -
2 L^2 (Cos[\[Theta]3[t]] + 1 - Cos[\[Theta]2[t]])] - len;(*spring 2 extension*)
pe = (1/2) k1 (spring1)^2 + (1/2) k2 (spring2)^2 - (1/2) m1 g L Cos[\[Theta]1[t]] - (1/2) m2 g L Cos[\[Theta]2[t]] - (1/
2) m3 g L Cos[\[Theta]3[t]];
lag = ke - pe;
eq1 = D[D[lag, \[Theta]1'[t]], t] - D[lag, \[Theta]1[t]] == 0;
eq2 = D[D[lag, \[Theta]2'[t]], t] - D[lag, \[Theta]2[t]] == 0;
eq3 = D[D[lag, \[Theta]3'[t]], t] - D[lag, \[Theta]3[t]] == 0;
ic = {\[Theta]1'[0] == initial\[Theta]1dot, \[Theta]1[0] == initial\[Theta]1, \[Theta]2'[0] ==
initial\[Theta]2dot, \[Theta]2[0] == initial\[Theta]2, \[Theta]3'[0] == initial\[Theta]3dot, \[Theta]3[0] ==
initial\[Theta]3};
sol = First@NDSolve[{eq1, eq2, eq3,
ic}, {\[Theta]1[t], \[Theta]2[t], \[Theta]3[t], \[Theta]1'[t], \[Theta]2'[t], \[Theta]3'[t]}, {t, from, to}];
{pe, ke, sol}
]
)
]