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. 

 *)

3. Oscillating Systems

Introduction

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

3.1 Linear Oscillations

Overview of Linear Oscillations

Initialization of User-Defined Functions:

1. phasePlot for one-dimensional system

\(phasePlot[q_, p_, tList_List, 
        opts___]\ \  := 
ParametricPlot[\ {p, q}\  // 
          Evaluate, 
tList, 
        opts, 
PlotRange -> 
          All, 
PlotStyle -> {Hue[Random[]], 
            Thickness[ .01]}, 
Axes -> 
          False, 
PlotLabel -> "Phase Trajectory", \ \ 
Frame -> 
          True, 
FrameLabel -> {"p", "q"}\ ]\)

phasePlot[Sin[2 π t], D[Sin[2 π  t], t], {t, 0, 3}, 
        GridLines -> Automatic] ;

2. timePhasePlot: Time behavior of phase plot for a one-dimensional system

\(timePhasePlot[q_, p_, tList_List, 
        opts___] := 
ParametricPlot3D[\ {p, q, 
            tList[\([1]\)], {\ Hue[ .8], Thickness[ .003]}}\  // 
          Evaluate, 
tList, 
        opts, 
PlotRange -> 
          All, 
PlotLabel -> "     Phase[t]", 
AxesLabel -> {"p[t]", "q[t]", "t "}, 
PlotPoints -> 
          300, 
BoxRatios -> {1, 1, 1}]\)

timePhasePlot[Sin[2 π t], 
        D[Sin[2 π  t], t], {t, 0, 3}];

3. doublePlot: Combination of the phase plot and the time behavior of the phase plot

doublePlot[q_, p_, 
          tList_List] := 
Show[
GraphicsArray[ { 
phasePlot[q, p, tList, DisplayFunction -> Identity, 
                FrameTicks -> None] , 
              timePhasePlot[q, p, tList, DisplayFunction -> Identity, 
                Ticks -> 
                  None] } 
] 
];

doublePlot[Sin[2 π t], 
        D[Sin[2 π t], t], {t, 0, 3}];

4.  Protect Commands

Protect[ phasePlot,timePhasePlot,doublePlot];

Problem 1: Linear Oscillator

Remarks and outline

Solution

Clear["Global`*"]

Part a

eq1=x''[t]+ω(0^2)  x[t]==0 ;

initial1={x[0] == x0,x'[0]==v0};

dSol1 =  DSolve[ {eq1 , initial1} // Flatten, x[t], 
          t]//Flatten

Needs["Graphics`Arrow`"]

?Graphics`Arrow`*

?Arrow

arrow =  {
{Hue[0.2], 
            Arrow[{0, 0}, {0.5, 0.5}]}, 
{Hue[0.4], 
            Arrow[{0, 0}, {2, 0}]}, 
{Hue[0.6], 
            Arrow[{0, 0}, {0, 
                0.5}]}
};

Plot[x[t] /. dSol1 /. {x0 -> 0, ω0 -> 1,  
            v0 -> 1}, {t, -2 π, 
          2 π}, 
Epilog -> arrow];

Part b

initial2 = {x[0] == 
            A Cos[δ], x'[
              0] == -A ω0 Sin[δ]};

dSol2 = DSolve[{eq1, initial2} // Flatten, x[t], t] // Flatten  // 
        Simplify

Part c

lapTx = LaplaceTransform[eq1, t, s]

lapSol=Solve[lapTx, 
          LaplaceTransform[x[t], t, s]][[1]]

InverseLaplaceTransform[
        LaplaceTransform[x[t], t, s] /. lapSol, s, 
        t]//Apart

Problem 2: Series expansion Solution

Clear["Global`*"]

    
diffSeriesOne[eq_, z_, t_, initialx_, initialv_, order_]  :=  
      Module[ {tempser, tempeq ,  tempvar,  tempsol, 
          tempinit }, 
tempser =  
          z -> Series[
              z, {t, 0, 
                order}]; 
tempeq = ( (eq[[1]] - 
                    eq[[2]])  //. tempser) == 0 ; 
        tempvar =  (Table[D[z, {t, j}], {j, 2, order + 2}] /. t -> 0) // 
            Flatten; 
tempsol = 
          Solve[tempeq,  tempvar] // 
            Flatten; 
tempinit =  {Head[z][0] -> 
              initialx, Head[z]'[0] -> 
              initialv}; 
z /. tempser /. tempsol /. 
          tempinit ]

Protect[diffSeriesOne ];

Solution

Clear["Global`*"]

Part a

Clear["Global`*"]

    
diffSeriesOne2[eq_, z_, t_, initialx_, initialv_, order_]  :=  
      Module[ { }, 
tempser =  
          z -> Series[
              z, {t, 0, 
                order}]; 
tempeq = ( (eq[[1]] - 
                    eq[[2]])  //. tempser) == 0 ; 
        tempvar =  (Table[D[z, {t, j}], {j, 2, order + 2}] /. t -> 0) // 
            Flatten; 
tempsol = 
          Solve[tempeq,  tempvar] // 
            Flatten; 
tempinit =  {Head[z][0] -> 
              initialx, Head[z]'[0] -> 
              initialv}; 
z /. tempser /. tempsol /. 
          tempinit ]

eq1 = y''[t] + γ y'[t]^2 +  w^2 y[t] == 
          a y'[t]^2  y[t];

diffSeriesOne2[eq1, y[t], t,   y0,  v0, 2]

tempser

tempeq//Short[#, 3] &

tempvar

tempsol // Short[#, 7] &

y[t] /. tempser /. tempsol /. 
      tempinit

diffSeriesOne[eq1, y[t], t,   y0,  v0, 2]

Part b

Clear["Global`*"]

eq2=x''[t]+ω(0^2)  x[t]==0 ;

approxsol = diffSeriesOne[eq2, x[t], t,   x0,  v0 , 
          9]//Collect[#1, {x0, v0}] &

termX[n_]= (((((-1))^n) (((t ω0))^(2 n)))/((2 n)!));
termV[n_]= (((((-1))^n) (((t ω0))^(2 n + 1)))/(ω0 (2 n + 1)!));

finitesum = 
        Sum[(x0 termX[n] + v0 termV[n]), {n, 0, 
            4}];

finitesum  == approxsol//Simplify

exactSum = Sum[(x0 termX[n] + v0 termV[n]), {n, 0, Infinity}] // 
          PowerExpand  // Apart

Problem 3: Potential and Phase Diagrams for the Linear Oscillator

Remarks and outline

Solution

Clear["Global`*"]

Part a

Force=m x''[t];
Spring =   -k x[
            t];
eq1 =  
      Force  ==  Spring

energy=Integrate[ x'[t]( m x''[t] +k x[t]) ,t]

T =  Integrate[ x'[t] Force , t];
V = -Integrate[ x'[t] Spring , t];
{T, V}

energy = T + V

Part b

Needs["Graphics`FilledPlot`"]

?FilledPlot

\(\(FilledPlot[
        V /. {x[t] -> x, k -> 1, m -> 1}, {x, \(-1\), 
          1}, 
PlotStyle -> {Thickness[ .02]\ }, 
Fills -> {Hue[ .8]}, 
Frame -> 
          True\ , 
RotateLabel -> 
          False\ , 
PlotLabel -> "V", 
FrameLabel -> {"x-axis", "Potential"}\ \ \ ];\)\)

Part c

values={m -> 1,k -> 1,x[t] -> x,x'[t]->v};

pt1 = 
        ContourPlot[
          Evaluate[energy /. values], {x, -1, 1}, {v, -1, 
            1}, 
ColorFunction -> 
            Hue,   
Contours -> 
            30, 
PlotPoints -> 
            30, 
ContourStyle  -> { 
              Thickness[0.006]}, 
Axes -> 
            True,  
AxesLabel -> {"x[t]", "v[t]\
"}];

Needs["Graphics`Graphics3D`"]

?ShadowPlot3D

\(\(ShadowPlot3D[
        energy /. values, {x, \(-1\), 1}, {v, \(-1\), 
          1}, 
PlotLabel -> "Energy Surface", 
Axes -> 
          True, 
AxesLabel -> {"x[t]", "v[t]", "energy"}, 
ShadowPosition -> \ 1];\)\)

Problem 4: Damped Linear Oscillator

Remarks and outline

Solution

Clear["Global`*"]

Part a

eq1=x''[t]+γ x'[t]+ ω0^2  x[t]==0;

xsol=x[t]/.DSolve[{eq1,x[0] == x0,x'[0]==v0},x[t],t][[1]] //Simplify

ωRule= Solve[(Sqrt[ω(0^2) - ((γ^2)/4)])==  ω,ω0][[2]]

xsol = xsol  /. ωRule // ExpToTrig  // 
        FullSimplify

xsol /. γ -> 0 // PowerExpand // 
      ExpandAll

Limit[xsol, γ -> 2  ω0]

Part b

values = {x0 -> 1,  ω0 -> 1, 
          v0 -> 0, ω -> (Sqrt[ω(0^2) - ((γ^2)/4)])} ;

pt[γ1_ /; 
          0 <= γ1 < 
            2]  := 
Plot[xsol //. 
              values /. {γ -> γ1 } // 
          Evaluate, 
{t, 0, 
          15} , 
PlotStyle -> {Hue[Random[]], 
            Thickness[ .03]}, 
PlotRange -> 
          All, 
DisplayFunction -> Identity ]

Show[
GraphicsArray[
{{pt[0], 
              pt[ .2],  pt[ .4]}, 
{ pt[ .8], pt[1.4], 
              pt[1.9]}}
]
];

Plot[ (xsol //. values) /. γ -> 
              Range[2.0001, 6.0001 , 0.4]  // 
          Evaluate, 
{t, 0, 
          15} , 
PlotStyle -> 
            Table[ Hue[i], {i, 0, 1, 1/12}] // 
          Evaluate
];

Part c

{x, 
        v} =  {xsol, D[xsol, t]} //. 
          values /. {γ -> 0.2}

doublePlot[x, v, {t, 0.01, 30}];

Part d

values2 = {γ -> 0.1, ω0 -> 1, x0 -> 1, 
          v0 -> 
            0, ω -> (Sqrt[ω(0^2) - ((γ^2)/4)])};

traj  = 
        ParametricPlot[{ xsol, D[xsol, t] } //. values2  // 
            Evaluate, 
{t, 0, 
            60}, 
PlotStyle -> {Thickness[ .008], 
              GrayLevel[0]}, 
DisplayFunction -> 
            Identity];

Needs["Graphics`PlotField`"]

Clear[x, v];

    
 pt1 = 
        PlotVectorField[{ v,  -0.5 v   -  x }  // 
            Evaluate, 
{x, -1, 1}, {v, -1 , 
            1}, 
DefaultColor -> 
            Hue[0.6], 
PlotPoints -> 
            20, 
AspectRatio -> 
            Automatic, 
Axes -> 
            Automatic, 
AxesOrigin -> {0, 
              0}, 
ScaleFunction -> (1 &), 
Epilog -> traj[[1]]];

Part e

eq1

lapTx = LaplaceTransform[eq1, t, s]

lapSol = Solve[lapTx, 
            LaplaceTransform[x[t], t, 
              s]][[1]]//Simplify

sol =  
      InverseLaplaceTransform[LaplaceTransform[x[t], t, s] /. lapSol, s, 
        t]

Problem 5: Damped Harmonic Oscillator and Driving Forces

Remarks and outline

Solution

Clear["Global`*"]

Part a

eq1=x''[t]+γ  x'[t]+ω(0^2)  x[t]==Q0  Cos[ωd t];

eq2 = x[t] /. DSolve[eq1, x[t], t] // First // Simplify

eq3 = eq2 /. Q0 -> 0

eq4 = (eq2 - eq3) // FullSimplify

Part b

cRule=Solve[{(Subscript[c, 1])==(1/2) (A0 + I B0),(Subscript[c, 2])==(1/2) (A0 - I B0)},{C[1], C[2]}][[1]]

ωRule= Solve[(Sqrt[ω(0^2) - ((γ^2)/4)])==  ω,ω0][[2]]

homogSol = eq3  //. cRule  /. ωRule  // 
            ComplexExpand // Simplify  // PowerExpand

eq5 =  A Cos[ωd t + ϕ] // TrigExpand

eq6 = eq4 - eq5 == 0

eq7={
eq6 /. {Sin[_ t] -> 0, 
            Cos[_ t] -> 1}, 
eq6 /. {Cos[_ t] -> 0, 
            Sin[_ t] -> 1}}

eq8 =  Solve[eq7, {A, ϕ} ][[4]] // 
          Simplify//PowerExpand

rule={(cos^(-1))(
                x_) -> (tan^(-1))(((Sqrt[1 - (x^2)])/x))};

solAϕ =  eq8 /. rule //  Simplify  // 
        PowerExpand

partSol =   A Cos[ωd t + ϕ] /. solAϕ

genSol = partSol+homogSol

Part c

values={A -> 3, Q0 -> 1, ωd -> 1, 
            A0 -> 0, 
            B0 -> 
              1, ω -> (Sqrt[ω(0^2) - ((γ^2)/4)])};

\(\(Plot3D[\(genSol //. values\) //. {γ -> 0.1} // 
          Evaluate, 
{t, 0, 30}, {ω0, 0.5, 
          1.5}, 
PlotPoints -> 
          40, 
AxesLabel -> {"Time", "ω0", "x[t]"}, 
FaceGrids -> {{0, 1, 
              0}, {0, \(-1\), 0}}, 
ViewPoint -> {1, 
            0.4, 0.6}];\)\)

\(\(\(
\)\(Plot3D[\ \(A /. solAϕ\) /. values\  // 
          Evaluate, 
{γ, \  .001,  .5}, {ω0,  .5, 1.5}, 
PlotPoints -> 
          25, 
PlotRange -> {0, 
            15}, 
AxesLabel -> {"γ", "ω0", "Amplitude"}, 
FaceGrids -> {{\(-1\), 0, 0}}, 
ViewPoint -> {\  .5, \ 1, 
            0.4}\ \ ];\)\)\)

Part d

pt[name_, 
          opts___] := 
Plot[name //. 
                values //. {γ -> 0.1, ω0 -> 2}  // 
            Evaluate, 
{t, 0, 
            80}, 
opts,  
Ticks -> 
            None, 
DisplayFunction -> 
            Identity];

Show[
Graphics[{
Rectangle[{0,   0}, { 1,   1}, 
pt[genSol, 
                PlotStyle -> Hue[0.8],   Ticks -> 
                  Automatic ]], 
Rectangle[{0.6, 0.7}, {1, 
                1}, 
pt[partSol, 
                PlotStyle -> Hue[0.4],   Background -> 
                  GrayLevel[0.9]]], 
Rectangle[{0.6, 
                0}, {1, 0.3}, 
pt[homogSol, 
                PlotStyle -> Hue[0.8], 
                Background -> 
                  GrayLevel[
                    0.9]]]
}]
];

?Rectangle

q = genSol //. 
            values //. {γ -> 0.1, ω0 -> 
              2};
p = D[q, t];
doublePlot[q, p, {t, 0, 80} ];

phasePlot[q, p, {t, 200, 250}, 
        AspectRatio -> Automatic];

3.2 Small Oscillations

Overview of  Small Oscillations and Normal Modes

Initialization of User-Defined Funtions for Small Oscillations and Normal Modes

1 Eigenvalues and eigenvectors for small oscillating systems

smallOsc[T_, V_, coord_List, equlib_List, 
        t_:  t] := 
Module[{q, p, dim, equlibRule, Vij, 
          Tij, matrix, ω2Sol, freq, vector, eqs, 
          vSol}, 
q = coord; 
p = 
          D[coord, t]; 
        dim = Length[q]; 
equlibRule = 
          Thread[coord -> equlib]; 
Vij = 
          Outer[D[V, #1, #2] &, q, q ]; 
Tij = 
          Outer[D[T, #1, #2] &,  p, p ]; 
matrix = 
          Vij - ω2 Tij /. 
            equlibRule; 
ω2Sol = 
          Simplify[
            Solve[Det[matrix] == 
                0, ω2]]; 
freq = ω2 /. ω2Sol; 
        vector = Array[a, dim]; 
eqs = 
          Simplify[
            ExpandAll[
              Thread[0 == matrix . vector]]]; 
vSol = 
          Table[
Solve[
                Join[eqs /. ω2Sol[[i]], {vector . vector == 
                      1}]
, vector][[2]], {i, 1, 
              dim}]; 
{freq, vSol} ]

Protect[smallOsc ];

T=(1/2) m ((q1'(t))^2)+(1/2) m ((q2'(t))^2);
V=(1/2) k1 (((q2(t) - q1(t)))^2);

smallOsc[T, V, {q1[t], q2[t]}, {0, 0}, t]

Problem 1: Two Coupled Oscillators along a Line

Remarks and outline

Solution

Clear["Global`*"]

Part a

T= (1/2) m x[1]'[t]^2 + (1/2) m x[2]'[t]^2;
V= (1/2) k2 x[1][t]^2 + (1/2) k2 x[2][t]^2 + 
            (1/2) k1 (((x[2][t] - x[1][t]))^2);

q = {x[1][t], x[2][t]};   
p = D[q, t];

Tij = Outer[D[T, #1, #2] &, p, p]
Vij = Outer[D[V, #1, #2] &, q, q]

(matrix = 
          Vij - ω2 Tij)//MatrixForm

ω2Sol = Simplify[Solve[Det[matrix] == 0, ω2]]

vector = Array[a, 2];
eq = Thread[0 == matrix . vector]

eq1 = Simplify[ExpandAll[eq /. ω2Sol[[1]]]]

eqNorm = {vector . vector == 1};

ev1rule = Solve[Join[eq1, eqNorm], vector][[2]]

eq2 = Simplify[ExpandAll[eq /. ω2Sol[[2]]]]

ev2rule = Solve[Join[eq2, eqNorm], vector][[2]]

Part b

eq3 = Thread[0 == matrix . vector]

(sol =  Reduce[
              Flatten[{eq3, eqNorm }], {ω2, a[1], a[2]}][[{8, 
                10}]])//ColumnForm

Part c

eq4 = smallOsc[T, V, {x[1][t], x[2][t]}, {0, 0}, t]

Part d

eigenmatrix = vector //. {ev1rule, ev2rule}

xRule = {x1[t], x2[t]} -> 
            Transpose[eigenmatrix] . {c1 Cos[ω1 t + ϕ1], 
                c2 Cos[ω2 t + ϕ2]} // Thread  // 
        Simplify

freqrule =  
      Thread[{ω1, ω2} -> (Sqrt[ω2] /. ω2Sol)]

values = {k1 -> 1, k2 -> 1, 
          m -> 1, ϕ1 -> 0, ϕ2 -> 0};

allRules = Join[xRule, freqrule, values];

pt1a = 
        Plot[{x1[t] + 1, x2[t] - 1} //. allRules //. {c1 -> 1, 
                c2 -> 0}  // Evaluate, 
{t, 0, 
            20}, 
PlotStyle -> {{Thickness[0.010], 
                Hue[ .6]}, {Thickness[0.005], 
                Hue[ .9]}} 
];

pt1b = 
        ParametricPlot[{x1[t], x2[t]} //. allRules //. {c1 -> 1, 
                c2 -> 0}  // Evaluate, 
{t, 0, 
            20}, 
AxesLabel -> {x1, x2}];

pt2a = 
        Plot[{x1[t] + 1, x2[t] - 1} //. allRules //. {c1 -> 0, 
                c2 -> 1}  // Evaluate, 
{t, 0, 
            20},  
PlotStyle -> {{Thickness[0.010], 
                Hue[ .6]}, {Thickness[0.005], 
                Hue[ .9]}}, 
DisplayFunction -> 
            Identity];

pt2b = 
        ParametricPlot[{x1[t], x2[t]} //. allRules //. {c1 -> 0, 
                c2 -> 1}  // Evaluate, 
{t, 0, 
            20}, 
DisplayFunction -> 
            Identity] ;

Show[GraphicsArray[{pt2a, pt2b}]];

pt3a = 
        Plot[{x1[t] + 1, x2[t] - 1} //. allRules //. {c1 -> 1, 
                c2 -> 1}  // Evaluate, 
{t, 0, 
            20},  
PlotStyle -> {{Thickness[0.010], 
                Hue[ .6]}, {Thickness[0.005], Hue[ .9]}} ];

pt3b = 
        ParametricPlot[{x1[t], x2[t]} //. allRules //. {c1 -> 1, 
                c2 -> 1}  // Evaluate, {t, 0, 50} ];

Problem 2: Three Coupled Oscillators along a Line

Remarks and outline

Solution

Clear["Global`*"]

Part a

T=(1/2) m ((((x[1])'[t])^2)+(((x[3])'[t])^2))+(1/2) M (((x[2])'[t])^2) ;
V = 
        (1/2) k ((((x[2][t] - x[1][t]))^2) + (((x[3][t] - x[2][t]))^2)) ;

solution=smallOsc[T, 
        V, {x[1][t], x[2][t], x[3][t]}, {0, 0, 0}, 
        t]

Part b

eigenmatrix = Array[a, 
          3]/.solution[[2]]

xRule = {x1[t], x2[t], x3[t]} -> 
          Transpose[eigenmatrix] . {c1 t + ϕ1, 
              c2 Cos[ω2 t + ϕ2], 
              c3 Cos[ω3 t + ϕ3]}  // Thread

freqrule = {ω1, ω2, ω3} -> (Sqrt[solution])[[1]]  // Thread

values = {k -> 1, m -> 1, 
          M -> 3, ϕ1 -> 0, ϕ2 -> 0, ϕ3 -> 
            0};

allRules = Join[xRule, freqrule, values];

xPlot[C1_, C2_, C3_] := 
        Plot[{x1[t] + 2, x2[t], x3[t] - 2} //. 
                allRules //. {c1 -> C1, c2 -> C2, 
                c3 -> C3}  // Evaluate, 
{t, 0, 
            20}, 
PlotStyle -> {
{Thickness[0.010], 
                Hue[ .3]}, 
{Thickness[0.009], 
                Hue[ .6]}, 
{Thickness[0.007], 
                Hue[ .9]}}, 
AxesLabel -> {t , 
              x}, 
DisplayFunction -> 
            Identity];

xLPlot[C1_, C2_, C3_, opts___] := 
        ParametricPlot3D[{x1[t], x2[t], x3[t]} //. 
                allRules //. {c1 -> C1, c2 -> C2, 
                c3 -> C3}  // Evaluate, 
{t, 0, 
            20}, 
BoxRatios -> {1, 1, 
              1}, 
AxesLabel -> {x1, x2, 
              x3}, 
opts, 
DisplayFunction -> Identity];

Show[
        GraphicsArray[{xPlot[1, 0, 0],  
            xLPlot[1, 0, 0  , Ticks -> None ]}]];

Show[
        GraphicsArray[{xPlot[0, 1, 0],  
            xLPlot[0, 1, 0 , Ticks -> None, 
              ViewPoint -> {1, 1, 
                  1}]}
]
];

Show[
        GraphicsArray[{xPlot[0, 0, 1], 
            xLPlot[0, 0, 1 , Ticks -> None ]}]];

Show[  xPlot[0, 1, 1] , 
        DisplayFunction -> $DisplayFunction];

Show[
        xLPlot[0, 1, 1, 
          ViewPoint -> {1 , 1, 
              1 }], 
DisplayFunction -> $DisplayFunction];

Part c

number = 1     ;

Do[
        ListPlot[
{{x1[t] - 2, 0}, {x2[t], 
                    0}, {x3[t] + 2, 0}} //. 
                allRules
 //. {c1 -> 0, 
                c2 -> 1, c3 -> 1} 
 // 
            Evaluate, 
PlotStyle -> {PointSize[
                0.05]}, 
PlotRange -> {{-4, 
                4}, {-1, 1}}], 
{t, 0, 1 - 1/number,  
          1/number}];

Problem 3: Three Coupled Oscillators along a Circle

Solution

Clear["Global`*"]

Part a

T=(1/2)    r^2( m   (((θ[1])'[t])^2)+m     (((θ[2])'[t])^2)+m1    (((θ[3])'[t])^2)) ;
V = 
        (1/2) k  r^2 ((((θ[2][t] - θ[1][t]))^2) +   (((θ[3][t] - θ[2][t]))^2)  +   (((θ[1][t] - θ[3][t]))^2) );

eq1 = smallOsc[T, 
            V, {θ[1][t], θ[2][t], θ[3][
                t]}, {0, 0, 0}, t]  // Simplify // PowerExpand

solution = Limit[ eq1, m1 -> m  ] // PowerExpand

Part b

eigenmatrix = Array[a, 3] /. solution[[2]]

solrule = {θ1[t], θ2[t], θ3[t]} -> 
            Transpose[eigenmatrix] . {c1 t + phase1, 
                c2 Cos[freq2 t + phase2], c3 Cos[freq3 t + phase3]}  // 
          Thread;

freqrule =  {freq1, freq2, freq3} ->  Sqrt[solution[[1]]]  // 
        Thread

valuerule = {k -> 1, m -> 1 , phase1 -> 0, 
          phase2 -> 0, phase3 -> 0} ;

solrule

xPlot[C1_, C2_, 
          C3_] := 
Plot[{θ1[t] + 2, θ2[
                        t], θ3[t] - 2}  /. solrule /. freqrule //. 
                valuerule //. {c1 -> C1, c2 -> C2, 
                c3 -> C3}  // Evaluate, 
{t, 0, 
            20}, 
PlotStyle -> {
{Thickness[0.012]}, 
{Thickness[
                  0.008]}, 
{Thickness[
                  0.004]}}, 
DisplayFunction -> 
            Identity];

Show[
        GraphicsArray[{ xPlot[0, 1, 0], xPlot[0, 0, 1] }]];

Problem 4:  Double Pendulum

Remarks and outline

Solution

Clear["Global`*"]

Part a

x1 = L1 Sin[θ1[#] ] &;
y1 = L1 Cos[θ1[#] ] &;
x2 = x1[#] + L2 Sin[θ2[#]] &;
y2 = y1[#] + L2 Cos[θ2[#]] & ;

L1 = L2 = len;

T1 = 1/2   m1 ( x1'[t]^2 + y1'[t]^2);
\(\(T2 = 1/2\ \ \ m2 \((\ \(x2'\)[t]^2 + \(y2'\)[t]^2)\)	;\)\ \)
V1 = - m1 g y1[t];
V2 = - m2 g y2[t];
T = T1 + T2;
V = V1 + V2;

solution=smallOsc[T, 
        V, {θ1[t], θ2[t]}, {0, 0}, t]

Part b

eigenmatrix = Array[a, 2] //. solution[[2]]

thetaRule = {θ1[t], θ2[t]} -> 
            Transpose[eigenmatrix] . {c1 Cos[ω1 t + ϕ1], 
                c2 Cos[ω2 t + ϕ2]}  // Thread  // 
        Simplify

freqrule = {ω1, ω2} -> Sqrt[solution[[1]]]  // 
        Thread

values = {m1 -> 3, m2 -> 1, len -> 1, 
          g -> 1, ϕ1 -> 0, ϕ2 -> 0};

allRules = Join[thetaRule, values, freqrule] ;

line[c1_, c2_] = 
        Line[{{0, 0}, {x1[t], -y1[t]}, {x2[t], -y2[t]}}] //. 
          allRules;

Show[Graphics[Table[N[line[1, 0]], {t, 0, 4, 0.1}]]];

Show[Graphics[Table[N[line[0, 1]], {t, 0, 2.5, 0.2}]]];

Show[
        Graphics[Table[N[line[0.5, 0.5]], {t, 0, 6, 0.2}]]];

Part c

frames = 1 ; 
Do[Show[Graphics[N[line[0.5, 0.5]], PlotRange -> {{-2, 2}, {0, -2}}]], {t, 1, 25, 25/frames}];

Problem 5:  Understanding the User-defined procedure smallOsc[]

Solution

Clear["Global`*"]

smallOsc2[T_, V_, coord_List, equlib_List, 
        t_:  t] := 
Module[{}, 
q = 
          coord; 
p = D[coord, t]; 
        dim = Length[q]; 
equlibRule = 
          Thread[coord -> equlib]; 
Vij = 
          Outer[D[V, #1, #2] &, q, q ]; 
Tij = 
          Outer[D[T, #1, #2] &,  p, p ]; 
matrix = 
          Vij - ω2 Tij /. 
            equlibRule; 
ω2Sol = 
          Simplify[
            Solve[Det[matrix] == 
                0, ω2]]; 
freq = ω2 /. ω2Sol; 
        vector = Array[a, dim]; 
eqs = 
          Simplify[
            ExpandAll[
              Thread[0 == matrix . vector]]]; 
vSol = 
          Table[
Solve[
                Join[eqs /. ω2Sol[[i]], {vector . vector == 
                      1}]
, vector][[2]], {i, 1, 
              dim}]; 
{freq, vSol} ]

T=(1/2) m ((q1'[t])^2)+(1/2) m ((q2'[t])^2); 
V = 
        (1/2) k1 (((q2[t] - q1[t]))^2);

smallOsc2[T, V, {q1[t], q2[t]}, {0, 0}, t]

{q, p, dim} // ColumnForm

equlibRule

{Vij, Tij} // ColumnForm

ω2Sol

Solve[Join[eqs /. ω2Sol[[1]], {vector . vector == 1}], 
      vector]

{freq, vSol}

3.3 Oscillating Circuits

Overview of Circuits

Problem 1:  Series RC  Circuit

Remarks and Outline

Solution

Clear["Global`*"];

Part a

eq1 = V[t] == R q'[t] + q[t]/C

dsol = DSolve[eq1, q[t], t][[1]]

Part b

V[t_] = V0 Cos[ω t];
q[t_] = q[t] /. dsol//Simplify

c1Sol = Solve[q[0] == q0, C[1]][[1]]

q[t_]= q[t] /.c1Sol//FullSimplify

Part c

values1 = {q0 -> 1, V0 -> 1, C -> 1, 
          R -> 1, ω -> 6};

Plot[q[t] /. values1  // Evaluate,
{t, 0, 6},
AxesLabel->{ "time", "charge"},
PlotLabel ->  "RC Loop"];

values2 = {q0 -> 1, V0 -> 1, 
          R -> 1, ω -> 6};

\(\(Plot3D[
        q[t] /. values2\  // Evaluate, 
{t, 0, 
          6}, 
{C, 0.01, 
          0.25}, 
AxesLabel -> {"Time", "Capacitance", "Charge"}, 
Ticks -> {Automatic, { .01,  .1,  .2\ }, 
            Automatic}, 
FaceGrids -> {{0, \(-1\), 
              0}}, 
ViewPoint -> { .3, \(-1\),  .3}, 
PlotPoints -> 35];\)\)

Problem 2:   Series RL Loop

Remarks and Outline

Solution

Clear["Global`*"];

Part a

eq1 = 
      V[t] == 
        L  q''[t] + R  q'[t]

\(\(\(V[t_] = V0\ Cos[ω\ t];\)\(
\)
    \)\)
initial = {q[0] == q0, q'[0] == i0};

qsol = 
      q[t] /. (DSolve[{eq1, initial} // Flatten, q[t], t] // 
              Flatten //  Simplify)

isol = D[qsol, t] // FullSimplify

Part b

values = {R -> 0.1, L -> 2, 
          V0 -> 1, ω -> 1, i0 -> 0.1, 
          q0 -> 0};

\(\(Plot[{qsol, isol} /. values\  // Evaluate, 
{t, 0, 
          20}, 
PlotStyle -> {
{Hue[ .3], \ Thickness[0.01]}, {Hue[ .9], Dashing[{0.03}], 
              Thickness[0.02]}}, 
Frame -> 
          True, \ 
FrameLabel -> {\ "Time"\ , \ " "\ \ }, 
GridLines -> Automatic];\)\)

Part c

<<Statistics`NormalDistribution`

fourierdata = 
        Table[qsol RandomArray[NormalDistribution[1, 0.1], 
                    1][[1]]  /. values // N, {t, 0, 100 π, 
            100 π/300}] ;

ListPlot[
Abs[
          Fourier[fourierdata]], 
PlotJoined -> 
          True, 
PlotStyle -> 
          Hue[Random], 
PlotRange -> {0, 
            5}];

Problem 3 :   RLC  Loop

Remarks and Outline

Solution

Clear["Global`*"];

Part a

eq1 = L q''[t] + R q'[t] + q[t]/C == 0;

eq2 = LaplaceTransform[eq1, t, s] // Simplify

sol0 = Solve[eq2, 
          LaplaceTransform[q[t], t, s]][[1]]

initial = {q[0] -> q0, q'[0] -> i0};

eq3 = 
      InverseLaplaceTransform[LaplaceTransform[q[t], t, s] /. sol0 /. 
            initial, s, t] // FullSimplify

freqrule=C (R^2) - 4 L -> -4 C (L^2) ω(0^2);

charge  =   eq3 //. freqrule // PowerExpand  // 
          ExpToTrig //  FullSimplify

current  =  D[charge , t] // FullSimplify

Part b

values = {R -> 1, C -> 1, L -> 4, q0 -> 1, 
          i0 -> 0};

ωrule=Solve[(Sqrt[C (R^2) - 4 L]) == 
                2 I (Sqrt[C]) L ω0,ω0][[1]]

text = {Text[
            StyleForm["Charge" , FontFamily -> "Times-Bold", 
              FontSize -> 15]  , {5, 0.8}] , 
Text[
            StyleForm["Current" , FontFamily -> "Times-Bold", 
              FontSize -> 15 ], {8, 0.3}]};

\(\(Plot[\({charge, current} /. ωrule\) /. values // 
          Evaluate, 
{t,  .001, 
          20}, 
PlotStyle -> {
{Hue[
                0.3], Dashing[{0.01}], Thickness[0.010]}, {Hue[0.9], 
              Dashing[{0.02}], 
              Thickness[0.008]}}, 
Frame -> 
          True, 
FrameLabel -> {"Time", ""}, 
GridLines -> 
          Automatic, 
Epilog -> text];\)\)

values2 = {R -> 1, C -> 1, q0 -> 1, 
          i0 -> 1};

\(\(p1 = 
        Plot3D[\(charge /. ωrule\) /. values2\  // 
            Evaluate, 
{t, 0.01, 
            30}, 
{L, 1, 
            4}, 
ViewPoint -> {1, \(+3\), 
              1}, 
FaceGrids -> {{0, 0, 
                1}, {0, \(-1\), 0}}, 
          AxesLabel -> {"Time", "Inductance", "Charge"}, 
PlotPoints -> 30];\)\)

3.4 Exercises

Exercise 1: Blocks connect with a Spring on a Table

Exercise 2: Falling Blocks connect with a Spring

Exercise 3: Parabolic coordinates and Harmonic Oscillator

Exercise 4: Forced Oscillator

Exercise 5: Pendulums and Spring

Exercise 6: Five Particles connected with Springs

Exercise 7: Vibrating Atoms on an Equilateral Triangle

Exercise 8: Vibrating Charged Spheres

Exercise 9:  Quadratic Damping

Exercise 10: Motion of a Charged Particle in a Harmonic Potential