Mathematica关于混沌摆的简单模拟

来源：互联网发布：2016网络情歌对唱大全编辑：程序博客网时间：2024/06/10 19:33

Clear["Global`*"]m = 1; l = 2; g = 10; tmax = 30; num = 200;func = {p1'[ t] == -(1/ 2) m l^2 (\[Theta]1'[t] \[Theta]2'[ t] Sin[\[Theta]1[t] - \[Theta]2[t]] + 3 g/l Sin[\[Theta]1[t]]), p2'[t] == -(1/ 2) m l^2 (-\[Theta]1'[t] \[Theta]2'[ t] Sin[\[Theta]1[t] - \[Theta]2[t]] + g/l Sin[\[Theta]2[t]]), \[Theta]1'[t] == 6/(m l^2) (2 p1[t] - 3 Cos[\[Theta]1[t] - \[Theta]2[t]] p2[t])/(16 - 9 Cos[\[Theta]1[t] - \[Theta]2[t]]^2), \[Theta]2'[t] == 6/(m l^2) (8 p2[t] - 3 Cos[\[Theta]1[t] - \[Theta]2[t]] p1[t])/(16 - 9 Cos[\[Theta]1[t] - \[Theta]2[t]]^2)};ini = {\[Theta]1[0] == Pi/2, \[Theta]2[0] == Pi/2, p1[0] == 0, p2[0] == 0};sol = First@ NDSolve[{func, ini}, {p1, p2, \[Theta]1, \[Theta]2}, {t, 0, tmax}];point1[t_] := Evaluate[{l Sin[\[Theta]1[t]], -l Cos[\[Theta]1[t]]} /. sol]point2[t_] := Evaluate[{l (Sin[\[Theta]1[t]] + Sin[\[Theta]2[t]]), -l (Cos[\[Theta]1[t]] + Cos[\[Theta]2[t]])} /. sol]data = Table[ Show[Graphics[{Red, Thick, Line[{{{0, 0}, point1[t]}, {point1[t], point2[t]}}]}, PlotRange -> {{-4.1, 4.1}, {-4.1, 1.2}}, ImageSize -> 300], ParametricPlot[point2[t0], {t0, 0, t}]], {t, tmax/num, tmax, tmax/num}];Animate[data[[i]], {i, 1, Length[data], 1}, AnimationRate -> num/20]Export["Doublependulum.gif", data]

0 0