Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

4. NonLinear Oscillating Systems

Introduction

Off[General::"spell"]; 
Off[General::"spell1"];

4.1 Nonlinear Pendulum

Overview of the Nonlinear Pendulum

Initialization of User-Defined Functions:

1. User-defined procedure for the Pendulum's angle-time graph

anglePend[{θ0_, dθ0_}, {γ_, 
            f0_, ω0_, ωd_}, tmax_, step_,  opts___]:=
Module[{temp},
temp=NDSolve[
{dθ'[
                                t] + γ  dθ[
                                  t] + ω0^2  Sin[θ[
                                    t]] == f0 Cos[t ωd] ,
θ'[t] ==dθ[t],θ[0] == θ0,dθ[0] == dθ0},
{dθ[t], θ[t]},
{t, 0, tmax},MaxSteps -> step 
]//Flatten ;

Plot[{ θ[t]  } /. temp   // 
                  Evaluate, 
{t, 0, 
                  tmax}, 
opts, 
DisplayFunction -> Identity] 

]

2. User-defined procedure for the Pendulum's Phase diagram

phasePend[{θ0_, dθ0_}, {γ_, 
          f0_, ω0_, ωd_}, tmax_, step_,  opts___]:=
Module[{temp},
temp=NDSolve[
{dθ'[
                              t] + γ  dθ[
                                t] + ω0^2  Sin[θ[t]] == 
                          f0 Cos[t ωd] ,  
θ'[t] ==dθ[t],θ[0] == θ0,dθ[0] == dθ0},
{dθ[t], θ[t]},
{t, 0, tmax},MaxSteps -> step 
]//Flatten ;

numgraph = 
ParametricPlot[
{ θ[t] , 
                      dθ[t]} /. temp  // 
                  Evaluate, 
{t, 0, 
                  tmax}, 
opts, 
DisplayFunction -> Identity] 
]

3. User-defined procedure for the Pendulum's Poincaré diagram

poincarePend[{θ0_, dθ0_}, {γ_, 
          f0_, ω0_, ωd_},  start_, tmax_, step_,  opts___]:=
Module[{temp1, temp2},
temp1=NDSolve[
{dθ'[
                              t] + γ  dθ[
                                t] + ω0^2  Sin[θ[t]] == 
                          f0 Cos[t ωd],  
θ'[t] ==dθ[t],θ[0] == θ0,dθ[0] == dθ0},
{dθ[t], θ[t]},
{t, 0, tmax},
MaxSteps -> step 
]//Flatten ;

temp2 =  
              Table[ { θ[t] , dθ[t]} /. temp1, {t, start, 
                  tmax, 2 π/ωd}];
ListPlot[temp2,  
              PlotStyle -> PointSize[0.02], opts]
  ]

4. User-defined procedure to map the pendulum to the interval -π  and  π

reduce[x_] = Mod[x ,  
        2  π, -π]

Plot[reduce[x], {x, -4 π , 4 π}, 
        Ticks -> {Table[i π, {i, -4, 4}], None}];

5. User-defined procedure for the Pendulum's reduced angle-time graph

anglePendR[{θ0_, dθ0_}, {γ_, 
          f0_, ω0_, ωd_}, tmax_, step_,  opts___]:=Module[{temp},temp=NDSolve[
{dθ'[
                              t] + γ  dθ[
                                t] + ω0^2  Sin[θ[t]] == 
                          f0 Cos[t ωd] ,  
θ'[t] ==dθ[t],θ[0] == θ0,dθ[0] == dθ0},
{dθ[t], θ[t]},
{t, 0, tmax},
MaxSteps -> step 
]//Flatten ;

Plot[
              reduce[θ[t] /. temp]  // 
                Evaluate, 
{t, 0, 
                tmax}, 
opts, 
DisplayFunction -> Identity]
]

6. User-defined procedure for the Pendulum's reduced phase diagram

phasePendR[{θ0_, dθ0_}, {γ_, 
          f0_, ω0_, ωd_}, tmax_, step_,  opts___]:=Module[{temp},temp=NDSolve[
{dθ'[
                              t] + γ  dθ[
                                t] + ω0^2  Sin[θ[t]] == 
                          f0 Cos[t ωd] ,  
θ'[t] ==dθ[t],θ[0] == θ0,dθ[0] == dθ0},
{dθ[t], θ[t]},
{t, 0, tmax},
MaxSteps -> step 
]//Flatten ;

ParametricPlot[{ 
                  reduce[θ[t] /. temp ], dθ[t] /. temp} // 
                Evaluate, {t, 0, 
                tmax}, 
opts, 
DisplayFunction -> Identity]
]

7. User-defined procedure for the Pendulum's reduced Poincaré diagram

poincarePendR[{θ0_, dθ0_}, {γ_, 
          f0_, ω0_, ωd_},  start_, tmax_, step_,  opts___]:=Module[{temp1, temp2},
temp1=
NDSolve[
{dθ'[
                              t] + γ  dθ[
                                t] + ω0^2  Sin[θ[t]] == 
                          f0 Cos[t ωd] ,  
θ'[t] ==dθ[t],θ[0] == θ0,dθ[0] == dθ0},
{dθ[t], θ[t]},
{t,  0, tmax},
MaxSteps -> step 
]//Flatten ;

temp2 =  
Table[ {reduce[ θ[
                        t] /. temp1] , 
                  dθ[t] /. temp1} , 
{t, start, 
                  tmax, 2 π/ωd}];

ListPlot[temp2, 
              PlotStyle -> PointSize[0.02], opts]
  ]

8. Protect User-defined procedures

Protect[anglePend,phasePend,poincarePend,
reduce,anglePendR,phasePendR,poincarePendR];

Off[Clear::wrsym]; 
Off[SetDelayed::write];

Problem 1: Analytic Solution for the Planar Pendulum

Remarks and Outline:

Solution

Clear["Global`*"]

Part a

NLeq = θ''[t] + g/L Sin[θ[t]] == 0

eq1 = 
      En == (θ'[t] NLeq[[1]]  // ExpandAll // 
            Integrate[#, t] &)

Erule = {En -> -g/ L  Cos[θ0]};

En /. 
        Erule/.{θ0 -> {0, π}}

\(\(pt = 
        ContourPlot[
          eq1[\([2]\)] /. {g -> 1, L -> 1} // 
            Evaluate, \ 
{θ[t], \(-4\) π, \ 
            4\ \ π}, 
{\(θ'\)[t], \(-3\), 
            3}\ , 
Contours -> 
            Range[\(-1\), 3, 
              0.5], \ 
FrameTicks -> {π\ \ Range[\(-4\), 4\ \ ], 
              Automatic}, 
ColorFunction -> 
            Hue, 
PlotPoints -> 
            100, \ 
FrameLabel -> {"θ", "θ'"}, 
DisplayFunction -> 
            Identity\ ];\)\)

Show[
        pt,  
Graphics[{Hue[ .2], PointSize[ .02], 
            Map[Point, {  π (2 # + 1), 0 }  & /@ 
                Range[-2, 1]]}], 
Graphics[{Hue[ .8], 
            PointSize[ .03], 
Map[
              Point, {2 π  #, 0 }  & /@ 
                Range[-2, 
                  2]]}], 
DisplayFunction -> $DisplayFunction];

cylind =  
        ParametricPlot3D[
{v, Sin[θ], 
            Cos[θ]}, 
{θ, -π, π}, 
{v, -3, 3}, 
DisplayFunction -> 
            Identity  ];

    
Show[cylind, DisplayFunction -> $DisplayFunction];

Show[  cylind, 
Graphics3D[
{AbsoluteThickness[
              3], 
Cases[Graphics[pt], _Line, 
                Infinity] /. 
{θ_ , v_} :> {v, 
                  1.02  Cos[θ ], 
                  1.02  Sin[θ ] }}], 
 
        ViewPoint -> {1, 1, 
            1}, 
DisplayFunction -> $DisplayFunction
];

Part b

eq2 = Solve[eq1, θ'[t]][[2]]  /. Erule // 
          Simplify // PowerExpand

eq3 =   t == 
            Integrate[  (1/(θ'))[t] //. eq2 , {θ[t], 
                0, θ}]  // Simplify // PowerExpand

eq4 = Solve[eq3, θ, 
          InverseFunctions -> 
            True ]//Flatten

NLθsol[
        t1_, θ1_, ω_ ]  = 
θ  /. Solve[eq3, θ, 
                InverseFunctions -> True][[1]]  /. {t -> 
              t1, θ0 -> θ1, 
            g ->   L (ω^2) } //PowerExpand

NLeq[[1]] /. θ -> ((NLθsol[#, θ0 , ω  ]) &) /. g -> (ω^2) L // Simplify

initial = 
{θ[0], θ'[
                0]} == 
{ 
              NLθsol[0, θ0, ω] , 
              D[NLθsol[t, θ0, ω], t]} /. t -> 0 // 
        Thread

list =  
        Table[ NLθsol[t , i π, 1 ], {i, 0.1 , 0.9 , 
            0.2}] ;

\(\(Plot[
        list // Evaluate, 
{t, 0, 
          16}, 
AxesLabel -> {"t", "θ"}, 
PlotStyle -> \((\({Thickness[#/100], 
                    Hue[#]} &\)\  /@ \ Range[0.1, 0.9, 0.2] // 
              Evaluate)\)
];\)\)

Part c

eq3[[2]]  /.  {θ -> θ0}

period =  4 eq3[[2]] /. {θ -> θ0}

Plot[ period  /. {g -> 1, L -> 1}  // Re  // 
          Evaluate, 
{θ0, 
          0.0 , π }, 
GridLines -> 
          Automatic, 
AxesLabel -> {"θ0", "\
Period"}, 
PlotStyle -> {Hue[ .6], 
            Thickness[ .01]}, 
Ticks -> {π/4 
               Range[0, 4, 1], π Range[0, 10, 
                2]}, 
PlotRange -> {{0, π}, {0, 
              10 π}}, 
Epilog ->  {Hue[ .8], 
            Thickness[0.01], 
            Line[{ {0, 2 π}, {π, 
                  2  π}} ]}
];

Part d

Leq = θ''[t] + g/L θ[t] == 0;

Lθsol[
        t1_, θo_, ω_ ] = θ[
              t] /. DSolve[{Leq, θ[0] == 
                    0, θ'[0] == 
                    2 ω Sin[θ0/2]}, θ[t], 
                t][[1]] /. {t -> 
              t1, θ0 -> θo , 
            g ->   L ω^2   } // PowerExpand

pt1 = 
        Plot[{  Lθsol[t, π/2, 1],  
              NLθsol[t, π/2, 1]}  // 
            Evaluate, 
{t, 0, 
            4 π}, 
PlotStyle -> {
{Dashing[{1.00, 0.00}], 
                Thickness[0.01]}, {Dashing[{0.06, 0.03}], 
                Thickness[0.001]}}, 
Frame -> 
            True,  
DisplayFunction -> 
            Identity
];

Lpts =  {#,  
                    Lθsol[#, π/2, 1]} &   /@   Range[0, 
                  4 π,  .5] // Flatten // Partition[#, 2] & // 
          Chop;

    
NLpts =  {#,  NLθsol[#, π/2, 1]} & /@ 
                Range[0, 4 π,  .5] // Flatten // Partition[#, 2] & // 
          Chop;

Show[
pt1, 
Epilog -> {Hue[ .6], PointSize[ .02], 
            Map[ Point, 
              Join[NLpts , 
                Lpts ]]}, 
DisplayFunction -> $DisplayFunction
];

Part e

?Fourier

fourierdataNL = 
        Table[NLθsol[t, 0.99  π , 1]  // N, {t, 0, 100, 
            100/511}] ;

ptNL = 
        ListPlot[
          Fourier[fourierdataNL]  // 
            Abs, 
PlotJoined -> 
            True, 
PlotStyle -> 
            Hue[Random[]], 
PlotRange -> 
            All
];

Show[ptNL, PlotRange -> {{0, 45}, {0, 8 }}];

fourierdataL = 
        Table[Lθsol[t,  .99  π , 1]  // N, {t, 0, 100, 
            100/511}] ;

ptL = 
        ListPlot[
          Fourier[fourierdataL]  // 
            Abs, 
PlotJoined -> 
            True, 
PlotStyle -> 
            Thickness[ .01],  
PlotRange -> {{0, 
                45}, {0, 10 }}, 
Epilog -> 
            ptNL[[1]]
];

Problem 2: Damped Pendulum

Solution

Clear["Global`*"]

Part a

eq1={θ''[t]+γ  θ'[t]+ω0^2 Sin[θ[t]]==0,
θ[
                    0] ==  .9  π,θ'[0]==0}/.{γ -> 0.2, ω0 -> 1}//Flatten;

solDP = 
        NDSolve[ eq1,  θ[t] , {t, 0, 40}] // Flatten;

\(\(Plot[\ \ θ[t] /. solDP\ \ \  // 
          Evaluate, 
{t, 0, 
          40}, 
AxesLabel -> {"t", "θ"}, 
PlotRange -> 
          All, 
PlotStyle -> {\ Thickness[0.01], 
            Hue[Random[]]}
];\)\)

angle[{θ0_, v0_}, {γ_, ω0_}, tmax_,  
        opts___]:=
Module[{numsol},

numsol=
NDSolve[
{θ''[t]+γ  θ'[t]+ω0^2 Sin[θ[t]]==0,
θ[
                          0] == θ0,θ'[0]==v0},θ[t],
{t, 0, tmax}];

numgraph = 
Plot[
θ[t] /. numsol, 
{t, 0, 
                  tmax}, 
opts, 
DisplayFunction -> Identity] 
]

\(\(\(grlist = \ \ 
\(angle[\ {0, #}\ , {0.3, 1}, 23, 
              AxesLabel -> {"t", "θ"}, 
GridLines -> 
                Automatic, 
PlotStyle -> {Thickness[ .01], Hue[Random[]]\ }, 
Ticks -> {π\ Range[0, 8, 
                      2], 2\ π\ \ Range[\(-2\), 
                      2\ ]}, 
PlotRange -> 
                All\ , 
DisplayFunction -> 
                Identity
] &\)\ \  /@ \ \ {\(-6\), \(-3\), \(-1\), 1, 3, 6}\ ;\)\(\ \ \)\)\)

Show[grlist,  
        DisplayFunction -> $DisplayFunction];

Part b

ParametricPlot[
{θ[t] /. solDP, 
            D[ θ[t] /. solDP , t] }  // 
          Evaluate, 
{t, 0, 
          40}, 
AxesLabel -> {"θ", "\
θ'"}, 
Ticks -> {π Range[-2, 2, 
                1/4], Automatic}, 
PlotRange -> 
          All, 
PlotStyle -> {Hue[ .8], 
            Thickness[ .009] }
];

phase[{θ0_, v0_}, {γ_, ω0_}, tmax_,  
        opts___]:=
Module[{numsol},

numsol=
NDSolve[
{v'[
                              t] + γ  v[
                                t] + ω0^2  Sin[θ[t]] == 0 ,  
θ'[t] ==v[t],
 θ[0] == θ0,v[0] == v0},
{v[t], θ[t]}, 
{t, 0, tmax}]//Flatten;

numgraph = 
              ParametricPlot[
{ θ[t] , v[t]} /. 
                    numsol  // Evaluate, 
{t, 0, 
                  tmax}, 
opts, 
DisplayFunction -> Identity]
]

\(\(pt1 = \ 
phase[{0, 4}, { .2, 1}, 
          40, 
PlotRange -> 
            All, 
AxesLabel -> {"θ", "θ'"}, 
PlotStyle -> {Thickness[ .009], 
              Hue[ .7]}, 
Ticks -> {\ π\ \ Range[
                  0, 4\ ], 
              Range[\(-2\), 4\ ]}, 
GridLines -> 
            Automatic, \ 
          PlotPoints -> 
            200, 
DisplayFunction -> $DisplayFunction
];\)\)

\(\(phlist = \ \ 
\(phase[{0, #\ }, { .3, 1}, 
              23, 
AxesLabel -> {"θ", "θ'"}, 
PlotStyle -> {Thickness[ .008], 
                  Hue[Random[]]\ }, 
Ticks -> {π\ Range[\(-8\), 8, 2], 
                  Range[\(-6\), 6\ ]}, 
PlotRange -> 
                All\ \ 
] &\)\ \  /@ \ \ {\(-6\), \(-3\), \(-1\), 1, 3, 6}\ ;\)\)

pt2 = 
        Show[  
phlist , 
GridLines -> Automatic, 
          DisplayFunction -> $DisplayFunction
];

list = {0, 3 + #} & /@ 
          Range[-0.51 , 0.51 , 0.2];

\(\(grlist = \ 
\(phase[#, { .2, 1}, 
              155, 
PlotPoints -> 
                100, 
AxesLabel -> {"θ", "θ'"}, 
PlotStyle -> {Thickness[ .005], 
                  Hue[ .9]}\ \ 
] &\)\  /@ \ 
          list\ ;\)\)

pt3 = 
        Show[
grlist, 
Graphics[ Point[{2 π #, 0} ] & /@ 
              Range[0, 2] ], 
PlotRange -> All, 
          Ticks -> {π Range[0, 4 ],   π /4 Range[0, 4 ]}, 
          DisplayFunction -> $DisplayFunction
];

Part c

Needs["Graphics`PlotField`"]

pt4 = 
        PlotVectorField[
{ 
            y, - .2 y  - 1  Sin[x]   }, 
{x, 0,  
            5  π}, 
{y, -1,  
            5}, 
PlotPoints -> 
            20, 
ScaleFunction -> ( .4  # &), 
ScaleFactor -> 
            None , 
DisplayFunction -> 
            Identity 
];

Show[ pt4, pt1, pt3, 
        DisplayFunction -> $DisplayFunction];

Problem 3: Periodic Solutions for the Driven Pendulum

Solution

Clear["Global`*"]

Part a

\(\(pt1 = \(anglePend[{\ 0, #\ \ \ π\ }, {\  .5\ \ , 1\ , 1, 2/3}, 
              100\ , 5000\ , 
AxesLabel -> {"θ", "θ'"}, 
PlotRange -> 
                All, 
Ticks -> {π\ Range[0, 20\ , 
                      5], π\ Range[0, 10\ , 
                      2\ ]}, \ 
PlotStyle -> 
                Hue[Random[]]\ 
] &\)\ \  /@ \ \ {1/2, 3/2, 
            2, 3};\)\)

Show[pt1, DisplayFunction -> $DisplayFunction] ;

Part b

\(\(pt2 = \ \(phasePend[{\ 0, #\ π}, {\  .5\ \ , 1\ , 1, 2/3}, 400\ , 
              5000\ , \ 
              AxesLabel -> {"θ", "θ'"}, 
PlotRange -> 
                All, 
Ticks -> {π\ Range[0, 
                      6\ \ \ ], \ Automatic}, 
              PlotPoints -> 150, 
PlotStyle -> 
                Hue[Random[]]\ \ \ 
] &\)\ \  /@ \ \ {1/2, 
            3/2, 2, 3};\)\)

Show[pt2, 
        DisplayFunction -> $DisplayFunction];

Part c

pt3 =  
poincarePend[
{ 
                0, # π}, {  .5  , 1 , 1, 2/3}, 100,  400 , 
              5000 , 
PlotRange -> {{-25, 
                    25}, {-3, 
                    3}} , 
DisplayFunction -> 
                Identity 
] &  /@   {1/2, 3/2, 2, 
            3};

Show[ 
GraphicsArray[Partition[pt3, 2]] ,  
        DisplayFunction -> $DisplayFunction
];

Show[ pt3,  DisplayFunction -> $DisplayFunction];

Problem 4: Looping Solutions for the Driven Pendulum

Clear["Global`*"]

Solution

Part a

\(\(\(eq1 = \ 
anglePend[\ {0, 0}\ , {\ 1/2, 1.44\ , 1, 
            2/3}, 50\ , 10000\ , AxesLabel -> {"t", "θ"}, 
          Ticks -> {Automatic, \(-π\)\ Range[0, 10, 
                  2\ \ ]\ }, \ 
PlotStyle -> \ {\ \ Thickness[ .01]}, \ 
          DisplayFunction -> $DisplayFunction\ \ 
]\ ;\)\(\ \ \)\)\)

eq2 =   
anglePendR[ {0, 0}  , { 1/2, 1.44, 
            1, 2/3}, 50, 6000 , 
          Ticks -> {π Range[0, 20 ,  5], 
              Range[-30, 0 , 
                5 ]}, 
PlotStyle ->   Thickness[ .01] 
] ;

Show[  eq1, eq2, 
        DisplayFunction -> $DisplayFunction];

Part b

eq3 =  
 #[{0, 0}  , 
{ 
                1/2, 1.44  , 1, 2/3}, 100 , 6000 , 
              PlotRange -> All, 
Ticks -> None, 
              PlotPoints -> 150,  
PlotStyle -> 
                Hue[Random[]]
]  &  /@   {phasePend, 
            phasePendR} ;

Show[GraphicsArray[eq3]];

Show[ 
eq3[[2]], 
PlotRange -> {{-1, 0 }, {-1, 1}}, 
        DisplayFunction -> $DisplayFunction
];

Part c

poincarePend[{0, 0}  , 
{ 1/2, 1.44 , 1, 2/3}, 
        400, 500 , 6000  ] ;

poincarePendR[{0, 0} , 
{ 1/2, 1.44 , 1, 2/3}, 
        400, 500 , 6000  ] ;

Problem 5: Chaotic motion for the Driven Pendulum

Solution

Clear["Global`*"]

Part a

eq1 = anglePend[
{ π/2, 
                0}, 
{ 0.2,  #, 1, 0.694}, 50 , 10000 , 
              PlotRange -> 
                All  
] &  /@   Range[0, 0.3, 
            0.1] ;

Show[
        GraphicsArray[
          Partition[eq1, 
            2]] , 
DisplayFunction -> $DisplayFunction];

eq2 =  phasePend[
{ π/2, 
                0}, 
{ 0.2, #, 1, 0.694} ,   600 , 
              10000 ,  
PlotRange -> 
                All , 
Ticks  -> {π Range[-1/2, 1/2, 1/4], None}, 
PlotPoints -> 
                100 
] &  /@   Range[0, 0.3, 
            0.1];

Show[GraphicsArray[Partition[ eq2, 2] ] , 
        DisplayFunction -> $DisplayFunction];

eq3 = poincarePend[
{ π/2, 
                0}, 
{ 0.2, #, 1, 0.694} ,  
              500,   600 , 
              10000 , 
Ticks  -> 
                False  , 
PlotRange -> {{-1 , 
                    1}, {-1, 
                    1}} , 
DisplayFunction -> 
                Identity
] &   /@   Range[0, 0.3, 
            0.1];

Show[GraphicsArray[Partition[ eq3, 2] ] , 
        DisplayFunction -> $DisplayFunction];

Part b

anglePend[
{ π/2, 0}, 
{ 
          0.2,  0.69, 1, 0.694}, 200 , 
        10000 , 
Ticks -> {Automatic, π Range[0, 
                2, 1/2 ]} , 
PlotPoints -> 
          150, 
DisplayFunction -> $DisplayFunction 
];

phasePend[
{ π/2, 0}, 
{ 
          0.2,  0.69, 1, 0.694},   700 , 
        8000 ,  
PlotRange -> 
          All,  
PlotPoints -> 
          200, 
DisplayFunction -> $DisplayFunction
];

poincarePend[
{ π/2, 0}, 
{ 
          0.2,  0.69, 1, 0.694},   200, 700 , 
        10000 , 
PlotRange -> {{-5 , 
              10}, {-3, 
              3}} ,   
DisplayFunction -> $DisplayFunction
];

Part c

\(\(anglePend[
{\ π/2, 0}, 
{\ 
          0.2, \ 0.74\ , 1, 0.694}, 300\ , 
        12000\ , 
AxesLabel -> {"t", "θ"}, 
Ticks -> {Automatic, π\ Range[\(-5\), 5, 1\ \ ]}\ , \ 
        DisplayFunction -> $DisplayFunction\ 
]\ ;\)\)

phasePendR[
{ π/2, 0}, 
{ 
          0.2 ,  0.83 , 1, 0.694},   300 , 
        10000 ,  
PlotRange -> 
          All,  
DisplayFunction -> $DisplayFunction
];

poincarePendR[
{ π/2, 
          0}, 
{ 0.2,  0.8, 1, 0.694}, 100,   300 , 
        10000 ,  
PlotRange -> 
          All,  
DisplayFunction -> $DisplayFunction
];

Part d

\(\(\(anglePend[
{\ π/2, 0}, 
{\ 
          0.2, \ 0.982\ , 1, 0.694}, 50\ , 
        12000\ , 
AxesLabel -> {"t", "θ"}, 
Ticks -> {Automatic, π\ Range[\(-1\), 1, 1/2\ \ ]}\ , 
        DisplayFunction -> $DisplayFunction\ 
];\)\(\ \)\)\)

\(\(\(anglePendR[
{\ π/2, 0}, 
{\ 
          0.2, \ 0.982\ , 1, 0.694}, 50\ , 
        12000\ , 
AxesLabel -> {"t", "θ"}, 
Ticks -> {Automatic, π\ Range[\(-1\), 1, 1/2\ \ ]}\ , 
        DisplayFunction -> $DisplayFunction\ \ 
]\ ;\)\(\ \ \ \)\)\)

phasePend[
{ π/2, 0}, 
{ 
          0.2,  1 , 0.982 , 0.694},   300 , 
        12000 , 
PlotPoints -> 
          150, 
DisplayFunction -> $DisplayFunction
];

pt1 =  
        phasePendR[
{ π/2, 0}, 
{ 
            0.2,  1 , 0.982 , 0.694},   300 , 
          12000 , 
PlotPoints -> 
            150, 
DisplayFunction -> $DisplayFunction
];

Show[ 
pt1, 
PlotRange -> {{ 1, 3 }, {-1, 
              1}}, 
DisplayFunction -> $DisplayFunction
] ;

poincarePendR[
{ π/2, 
          0}, 
{ 0.2,  1 , 0.985, 0.694}, 200, 300 , 
        12000  ,  
PlotRange -> {{-3 , 
              3}, {-3, 
              3}} ,  
DisplayFunction -> $DisplayFunction
];

4.2 Duffing Equation

Overview of the Duffing Equation

Initializations for the Duffing Equation

1. User-defined procedure to plot the Duffing displacement motion

dispDuff[{q0_, p0_}, {a_, b_, c_, f0_, wd_}, tmax_, step_,  
          opts___]:=Module[{temp},
temp=
NDSolve[
{p'[t] - a q[t] + b q[t]^3 + 
                              c  p[t]  == f0  Cos[t wd],  
q'[t] ==p[t],
q[0] == q0,
p[0] == p0},
{p[t], q[t]},
{t, 0, tmax},
MaxSteps -> step 
]//Flatten ;

numgraph =  
Plot[ 
q[t] /. temp  // Evaluate, 
{t, 0, 
                    tmax}, 
opts, 
DisplayFunction -> Identity
]

]

2. User-defined procedure to plot the Duffing phase

phaseDuff[{q0_, p0_}, {a_, b_, c_, f0_, ωd_}, tmax_, 
        step_,  opts___]:=Module[{temp},
temp=
NDSolve[
{p'[t] - a  q[t] + b q[t]^3  + 
                            c p[t] == f0 Cos[t ωd] ,
q'[t] ==p[t],
q[0] == q0,
p[0] == p0},
{p[t], q[t]},
{t, 0, tmax},
MaxSteps -> step 
]//Flatten ;

ParametricPlot[
{ 
                    q[t] , p[t]} /. temp  // 
                Evaluate, 
{t, 0, 
                tmax}, 
opts, 
DisplayFunction -> Identity] 

]

3. User-defined procedure to plot the Duffing Poincaré map

poincareDuff[{q0_, p0_}, {a_, b_, c_, f0_, ωd_},  
        start_, tmax_, step_,  opts___]:=
Module[{temp1, temp2},
temp1=
NDSolve[
{p'[t] - a  q[t] + b q[t]^3 + 
                            c  p[t]  == f0 Cos[t ωd],
q'[t] ==p[t],
q[0] == q0,
p[0] == p0},
{p[t], q[t]},
{t, 0, tmax},
MaxSteps -> step 
]//Flatten ; 

temp2 =  
Table[ { 
                    q[t] , p[t]} /. temp1  , 
{t, start, 
                  tmax, 2 π/ωd}
];

ListPlot[temp2, 
              PlotStyle -> PointSize[0.02], opts]

]

4. Protect User-defined procedures

Protect[dispDuff,phaseDuff,poincareDuff];

Off[Clear::wrsym]

Problem 1:  Potential and Phase Diagrams for the Duffing Oscillator

Remarks and outline

Solution

Clear["Global`*"]

Part a

Acc=q''[t];

Flinear =  a q[t]; 
Fnonlinear = - b q[t]^3;

eqMotion =  Acc - Flinear - Fnonlinear == 0

energy=Integrate[ q'[t] ( 
            Acc - Flinear - Fnonlinear) , t]

\(\(\({T = \ Integrate[\ \(q'\)[t] Acc\ , t], 

      V0 = \(-Integrate[\ \(q'\)[t] Flinear, t]\), 

      V1 = \(-Integrate[\ \(q'\)[t] Fnonlinear\ , t]\)}\)\(
\)
    \(\ \)\)\)

\(\(Plot3D[
        V0 + V1 //. {q[t] -> q, b -> 0.05}\  // 
          Evaluate, 
{q, \(-15\), 
          15}, 
{a, \(-1\), 
          4}, \ 
AxesLabel -> {"q", "a", "Potential"}, 
BoxRatios -> {1, 1, 
            1}, 
PlotRange -> {\(-100\), 
            50}, 
ViewPoint -> {0, 3, 
            1}
];\)\)

Needs["Graphics`Legend`"]

\(\(Plot[\(V0 + V1 //. {q[t] -> q, b -> 0.05}\) /. \ 
            a -> {\(-1\), 4}\  // 
          Evaluate, 
{q, \(-15\), 
          15}, \ 
Ticks -> 
          None, 
PlotStyle -> {
{Dashing[\ {0.01}\ ], Hue[0.6]}, 
{Dashing[\ {0.03}\ ], 
              Hue[0.9]}}, 
PlotLegend -> {" a=-1", " a=4"}, \ \ \ \ 
LegendPosition -> {0.1, 
            0.3\ }]\ ;\)\)

Part b

values={q[t] -> q,q'[t]->v,b -> 0.05};

contour =  
ContourPlot[energy //. 
                    values //. { a ->  # }  // 
                Evaluate, 
{q, - 15,  
                15}, 
{v, -15, 
                15}, 
Contours -> 
                10 ,  
PlotPoints -> 
                20, 
ColorFunction -> 
                Hue,    
DisplayFunction -> 
                Identity
] &   /@   {-1, 
            4};

Show[GraphicsArray[contour]];

pt1 =  
        ContourPlot[energy //. values //. { a ->  4 }  // 
            Evaluate, 
{q, -15, 
            15}, 
{v, -15, 
            15}, 
Contours -> 
            10 , 
DisplayFunction -> 
            Identity 
];

pt2 =   Plot3D[energy //. values //. { a ->  4  }  // 
            Evaluate, 
{q, -15, 
            15}, 
{v, -15, 
            15},   
PlotPoints -> 
            20,  
ColorFunction -> (If[# < 0.226 , 
                  Hue[0.9], 
                  Hue[0.1]] &), 
DisplayFunction -> 
            Identity
 ];

Show[ 
pt2, 
Graphics3D[
{AbsoluteThickness[
              3], 
Cases[Graphics[pt1] , _Line, 
                Infinity] /. 
{ q_Real , v_} :> {q, v, 
                  15  - 2 (q^2) + 
                    1/2  (v^2)  +  .0125 (q^4)}}
], 
ViewPoint -> {0, 1, 
            0.8}, 
DisplayFunction -> $DisplayFunction
];

Part c

eq1={eqMotion/.{q''[t]->p'[t]},q'[t]==p[t]}

initial = {p[0] == 0.001, q[0] == 0};

eq2 = Join[eq1, initial] ;

ndsol1 = 
NDSolve[
          eq2 //. {b -> 0.05, a -> -1}, {q[t], p[t]}, {t, 0, 
            30}] // Flatten

q1[t_] = q[t] //. ndsol1; 
p1[t_] = p[t] //. ndsol1;

pt[var_, opts___]   :=  
      Plot[var[t], {t, 0, 30}, opts,  
        DisplayFunction -> Identity]

Show[
GraphicsArray[
{pt[q1, 
              PlotStyle -> Hue[0.7] ] , 
pt[p1, 
              PlotStyle -> 
                Hue[0.9] ] }
]
];

Part d

ndsol2 = 
        NDSolve[eq2 //. {b -> 0.05, a -> +4}, {q[t], p[t]}, {t, 
              0, 30}] // Flatten;

q2[t_] = q[t] //. ndsol2; 
p2[t_] = p[t] //. ndsol2;

Show[
        GraphicsArray[
{pt[q2, 
              PlotStyle -> Hue[ .7] ] , 
pt[p2, 
              PlotStyle -> 
                Hue[ .9] ] }
]
];

\(\(ParametricPlot[
{\ q2[t]\ , p2[t]}\ \ \  // 
          Evaluate, 
{t, 0, 
          30}, 
AxesLabel -> {"q", "p"}\ , 
AspectRatio -> 
          Automatic
];\)\)

Part e

oscMovie[f_, tMax_, frames_, range_:  Automatic] := 
        Module[{t, 
            i}, 
Do[
t = 
              i tMax/frames; 
plot[
                i] = 
Show[
Graphics[
{PointSize[0.05], 
                    Point[{N[f[t]], 
                        0}]}, 
PlotRange -> {range, {-1, 1}}], 
Axes -> 
                  True
], 
{i, 1, 
              frames}]
];

frames = 1;
oscMovie[q2, 30, frames, {-15, 15}];

Problem 2:  Phase Diagram and Orbits for the Damped,  Duffing Equation

Solution

Clear["Global`*"]

Part a

Needs["Graphics`PlotField`"]

lista = 
PlotVectorField[ 
{p, # q -   .05  q^3 -  .15  p  }, 
{q, -13, 
                13}, 
{p, -20  , 
                20 } , 
PlotPoints -> 
                10, 
ScaleFunction -> ( .2  # &), 
ScaleFactor -> 
                None, 
DisplayFunction -> 
                Identity 
] &   /@   { 
            4 , -1};

pt1 = Show[GraphicsArray[lista]];

Part b

values = {{-4, 4 }, {4, 4} , {4, 
            23}, {-4, -23}};

listb = 
dispDuff[ 
#  , 
{4,  .05,  .1, 0, 0}, 30, 
              1000 , 
PlotStyle -> 
                Hue[Random[]], 
PlotRange -> 
                All
] &   /@   values;

Show[
listb[[{1, 
              2}]], 
PlotRange -> {{0, 
              30}, {-14, 
              14}}, 
DisplayFunction -> $DisplayFunction 
];

Show[
listb[[{3, 
              4}]], 
PlotRange -> {{0, 
              30}, {-14, 
              14}}, 
DisplayFunction -> $DisplayFunction 
];

listc = 
phaseDuff[ 
#  , 
{4,  .05,  .1, 0, 0}, 30, 
              1000 , 
PlotStyle -> 
                Hue[ .6], 
PlotRange -> 
                All
] &   /@   values;

Show[GraphicsArray[Partition[ listc, 2]], 
        DisplayFunction -> $DisplayFunction ];

Show[
{  listc[[4]] , 
          lista[[1]] }, 
PlotRange -> {{-15, 
              15}, {-20, 
              20}}, 
DisplayFunction -> $DisplayFunction
];

Problem 3: Driven Duffing Orbits with No Damping

Solution

Clear["Global`*"]

Part a

dispDuff[
{2.038, -1.5 }, 
{ -1 , 1 , 0 , 1/4, 2 }, 100, 
        10000 , 
PlotPoints -> 
          200, 
DisplayFunction  -> $DisplayFunction 
];

pt1 =  
phaseDuff[
{2.038, -1.5 }, 
{ -1 , 1 , 0 , 1/4, 2 }, 100 , 
          10000 , 
PlotPoints -> 
            300, 
DisplayFunction  -> $DisplayFunction 
];

phaseDuff[
{2.038, -1.5 }, 
{-1 , 1 , 0 , 1/4, 2}, 100, 
        10000 , 
PlotPoints -> 
          300, 
PlotRange -> {{1.93, 2.2}, {-3 , 
              3 }}, 
DisplayFunction -> $DisplayFunction 
];

list =  
phaseDuff[
{2.038, -1.5 }, 
{ -1 , 1 , 0 , 1/4, 2 }, #, 
              10000 , 
PlotPoints -> 
                300, 
PlotRange -> {{1.93, 
                    2.2}, {-3 , 
                    3 }}  
] &   /@   {30, 
            50};

Show[GraphicsArray[ list ]];

pt2 =  
poincareDuff[
{2.038, -1.5 }, 
{ -1 , 1 , 0 , 1/4, 2 }, 0 , 300, 
          10000 , 
DisplayFunction  -> $DisplayFunction 
];

Show[pt1, pt2, 
        PlotRange -> {{1.8, 2.4}, {-3 , 3 }} ];

Part b

dispDuff[
{ .99 , 0 }, 
{1/5 , 
          1 , 0 , 1/4, 2}, 200, 
        10000 , 
PlotPoints -> 
          200, 
DisplayFunction  -> $DisplayFunction 
];

pt3 =  
phaseDuff[
{ .99  , 
            0 }, 
{  1/5 , 1 , 0 , 1/4, 2 }, 300, 
          10000 , 
PlotPoints -> 200, 
          DisplayFunction  -> $DisplayFunction 
];

pt4 =  
poincareDuff[ 
{ .99  , 0 }, 
{  1/5 , 1 , 0 , 1/4, 2 }, 0 , 350, 
          10000  , 
          DisplayFunction -> $DisplayFunction 
];

Problem 4: Two Well Driven Duffing Oscillator with Damping

Solution

Clear["Global`*"]

Part a

lista  = 
dispDuff[
{ - .13, - .7}, 
{  .5,  .5,  .3, #, 1 }, 50 , 10000, 
              PlotStyle -> 
                Hue[Random[]], 
PlotPoints -> 
                200   
] &   /@    { .05  ,  .24 };

Show[GraphicsArray[Partition[lista, 2]]];

lista = 
phaseDuff[
{ - .13, - .7}, 
{  .5,  .5,  .3, #, 1 }, 300 , 
              10000, 
PlotPoints -> 
                200   
] &   /@   { .05  ,  .24 };

Show[GraphicsArray[ lista ]];

listb = 
poincareDuff[
{ - .13, - .7}, 
{  .5,  .5,  .3, #, 1 },  100, 400 , 
              10000, PlotRange -> {{-2, 2}, {-1, 
                    1}}, 
Ticks -> 
                False,  
DisplayFunction -> 
                Identity 
] &   /@   { .05 ,  .24 };

Show[GraphicsArray[ listb]  ];

Part b

listc  = 
dispDuff[
{ -0.13, -0.7}, 
{ 0.5, 0.5, 0.3, #, 1 }, 50 , 10000, 
              PlotStyle -> 
                Hue[Random[]], 
PlotPoints -> 
                200   
] &   /@   { 0.29 , 
            0.32};

Show[GraphicsArray[ listc ]];

listd = 
phaseDuff[
{ -0.13, -0.7}, 
{ 0.5, 0.5, 0.3, #, 1 }, 300 , 
              10000, 
PlotPoints -> 
                200   
] &   /@   { 0.29 , 
            0.32};

Show[GraphicsArray[ listd ]];

liste = 
poincareDuff[
{ -0.13, -0.7}, 
{ 0.5, 0.5, 0.3, #, 1 },  100, 
              400 , 10000, 
              PlotRange -> {{-2, 2}, {-1, 
                    1}}, 
Ticks -> 
                False, 
 
              DisplayFunction -> 
                Identity 
] &   /@   { 0.29, 
            0.32};

Show[GraphicsArray[ liste ]];

4.3 Exercises

Exercise 1:  van der Pol  Oscillator and Limiting Cycles   (Solved in the Appendix)

Exercise 2:  Springs

Exercise 3: Buckling column system

Exercise 4: Nonlinear Equation

Exercise 5:   Inverted Pendulum

Exercise 6: Driven Nonlinear Equation