I am trying to plot some results of chaotic behavior in a system based on a paper for one of my classes. I am a physics undergrad. I found some sample code online springy.nb which I have been trying to manipulate to fit my situation. Here is my code so far. I think the problem is my non-understanding of getting the x to iterate. Any help would be appreciated. Thank you!

springalt[{c_, a_, x_}, {rs_, vs_}, 
  tmax_] := {ϕ[t], ω[t], Mod[θ[t], 2*π]} /. 
  NDSolve[{ϕ'[t] == ω[t], ω'[
       t] == -c (((
            Sqrt[(Cos[ϕ[t]] + a)^2 + (Sin[ϕ[t]])^2] - 3)/
            Sqrt[(Cos[ϕ[t]] + a)^2 + (Sin[ϕ[t]])^2])*a*
             t]] + (((√((x0 - 
                    Cos[ϕ[t]])^2 + (0.6*Sin[(2*π)/3] - 
                    Sin[ϕ[t]])^2) - 
                3)/(√((x - 
                    Cos[ϕ[t]])^2 + (0.6*Sin[(2*π)/3] - 
                    Sin[ϕ[t]])^2)))*((0.6*Sin[(2*π)/3] - 
                 Sin[ϕ[t]])^2*Cos[ϕ[t]] - 
              x*Sin[ϕ[t]]))) - ω[t], θ'[t] == (
      2*π)/3, ϕ[0] == rs, ω[0] == vs}, {ϕ[
      t], ω[t], θ[t]}, {t, 0, 1 tmax}][[1]]
poincarealt[{c_, a_, x_}, ics_, tmax_] := 
 Block[{sol = springalt[{c, a, x}, ics, tmax]}, 
  Print[ParametricPlot3D[{sol[[1]] Cos[sol[[3]]], 
     sol[[1]] Sin[sol[[3]]], sol[[2]]}, {t, 0, tmax}]]; 
  Print[Plot[sol, {t, 0, tmax}]]; 
      Table[{t, sol}, {t, 0, tmax, .001}], #[[2, 3]] < 10^-2 &], 
     SameTest -> (Abs[First[#1] - First[#2]] < .1 &)] PlotStyle -> 
  • $\begingroup$ What are good values for the parameters? Probably all you need to do is call poincarealt[{c_, a_, x_}, ics_, tmax_] with appropriate values. $\endgroup$
    – bill s
    Apr 14 '16 at 20:52
  • $\begingroup$ c=10, a=6 I know the attracting x0 is around 6.38 $\endgroup$ Apr 14 '16 at 20:58
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Apr 16 '16 at 21:21

I corrected errors to get you started, but the physical equations themselves may not be accurate.

springalt[c_, a_, x_, rs_, vs_, tmax_] := NDSolve[{ϕ'[t] == ω[t], ω'[t] == 
    -c (((Sqrt[(Cos[ϕ[t]] + a)^2 + (Sin[ϕ[t]])^2] - 3)/
    Sqrt[(Cos[ϕ[t]] + a)^2 + (Sin[ϕ[t]])^2])*a*                     
    Sin[ϕ[t]] + (((\[Sqrt]((x - Cos[ϕ[t]])^2 + (0.6*Sin[(2*π)/3] - 
    Sin[ϕ[t]])^2) - 3)/(\[Sqrt]((x - 
    Cos[ϕ[t]])^2 + (0.6*Sin[(2*π)/3] - Sin[ϕ[t]])^2)))*((0.6*Sin[(2*π)/3] - 
    Sin[ϕ[t]])^2*Cos[ϕ[t]] - x*Sin[ϕ[t]]))) - ω[t], θ'[t] == (2*π)/3,
    ϕ[0] == rs, ω[0] == vs, θ[0] == 0}, {ϕ[t], ω[t], θ[t]}, {t, 0, tmax}][[1]]

tmax = 10; sol = springalt[10, 6, 6.38, 1, 1, tmax];
Plot[Evaluate[{ϕ[t], ω[t], Mod[θ[t], 2*π]} /. sol], {t, 0, tmax}]
ParametricPlot3D[{ϕ[t] Cos[θ[t]], ϕ[t] Sin[θ[t]], ω[t]} /. sol, {t, 0, tmax}]

enter image description here

enter image description here

Note that no initial condition was provided for θ[t], so I chose θ[0] == 0. Union[...] probably is meant to be part of a ListPlot calling sequence to produce a return map, but insufficient information is provided to correct it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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