Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Some fellow students and I are trying to calculate and plot the Lyapunov spectrum for the Duffing oscillator (currently a test approach). Each of us has Mathematica knowledge based purely on past 'need to know' experiences and are having troubles with getting the following code to work.
Can anyone point out where we have gone wrong (probably in multiple places)?

solDuffingOscillator[A_, B_, C_, x0_, x1_, maxt_] :=   
 NDSolve[{D[x[t], {t, 2}] + 2 A D[x[t], {t, 1}] + B x[t] + x[t]^3 ==   
    C Cos[ t], x[0] == x0, x'[0] == x1},  
  {x[t], x'[t]},  
  {t, 0, maxt},   
  MaxSteps -> Infinity]  

Manipulate[  
 Module[{sol1, sol2, X0, V0, X1, V1},  
  sol1 = solDuffingOscillator[A, B, C, x0, v0, tlim];  
  X0[t_] = First[x[t] /. sol];  
  V0[t_] = First[x'[t] /. sol1];  

  sol2[t_] =   
   solDuffingOscillator[A, B, C, (X0[t] + dx), (V0[t]), tlim];  
  X1[t_] = First[x[t] /. sol2[t]];  
  V1[t_] = First[x'[t] /. sol2[t]];  

  d = First[Evaluate[{X0}]] - First[Evaluate[{X1}]];  
  L = 1/dt*Log[d/dx];  

  Tbl1 = Table[{t, d}, {t, dt, tlim, dt}];  
  ListPlot[Tbl1, PlotRange -> {{0, tlim }, {-2 Pi, 2 Pi}},  
   ImageSize -> {1200, 800}]  
  ],  
 {{A, 0.5, "A"}, 0, 10, 0.001, Appearance -> "Labeled"},  
 {{B, -8.54, "B"}, -10, 0, 0.01, Appearance -> "Labeled"}, 
 {{C, 14, "C"}, 0, 20, 0.05, Appearance -> "Labeled"},  
 {{dt, 1, "dt"}, 0.001, 1, 0.0001, Appearance -> "Labeled"},  
 {{x0, 1, "x0"}, -2, 2, 0.01, Appearance -> "Labeled"},  
 {{v0, 0, "v0"}, -1, 1, 0.0001, Appearance -> "Labeled"},  
 {{dx, 0.5, "dx"}, -0.5, 0.5, 0.001, Appearance -> "Labeled"},  
 {{tlim, 1000, "t Cutoff"}, 100, 4000, 5, Appearance -> "Labeled"}]  
share|improve this question
    
Aside from a missing 1 from sol in the line X0[t_] = First[x[t] /. sol];, thing go bad in the 2nd call to solDuffingOscillator (to get sol2[t]). Try it first for fixed values of the parameters: subs = {a -> 0.5, b -> -8.54, c -> 14, dt -> 1, x0 -> 1, v0 -> 0, dx -> 0.5, tlim -> 1000};, sol1 = solDuffingOscillator[a, b, c, x0, v0, tlim] /. subs; X0[t_] = First[x[t] /. sol1], V0[t_] = First[x'[t] /. sol1], then solDuffingOscillator[a, b, c (X0[t] + dx), (V0[t]), tlim] /. subs`. –  murray Apr 18 at 20:09

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.