Numerical Computation of the Equations of Motion and Much More For The Damped Driven Pendulum

Introduction

The main function which drives this whole package is in the Runge function.  This function will integrate up to a third order ordinary differential equation.  This function takes the input of initial and final times, initial values for the derivative functions as well as the total number of steps.  The other integral function which works with Runge, is SetPendParam.  This function sets the three derivative functions which are later evaluated in the Runge function.  Therefore Runge will not be able to calculate without SetPendParam being called prior.

SetPendParam and Runge called together, then give output in table form.  This is all well and good, but there are better ways to recieve the output.  In the case of the chaotic pendulum, the two variables being solved for are θ and ω, where those represent the angle and angular momentum of the pendulum as it proceeds dynamically through the third output variable t, time.  There are two other functions, ThetaPlot and OmegaPlot, which take the table form solutions of Runge, and output them in ListPlot form.

This finds a solution for θ with t in the range 0 to 10, using an initial condition for θ and ω at t==1 for a set of parameters leading to periodic behavior.

<< Graphics`Pend`

ThetaPlot[{0.9, 2/3, 0.5, 0}, {0, 100, 1000}, {1, 1}]

This finds a solution for θ with t in the range 0 to 10, using an initial condition for θ and ω at t==1 for a set of parameters leading to chaotic behavior.

ThetaPlot[{1.5, 2/3, 0.5, 0}, {0, 100, 1000}, {1, 1}]

Testing

Here I will run some test cases on different differential equations which have analytic solutions.  From here I show that both my ODE solver is able to calculate correctly, efficiently, with as many cases as I can think to run.  To start, I use a very basic differential equation, (dx/dt)= -x ,  this equation has the solution x(t) = Exp[-x].  With this information I plot out the solution to this DE as well as the analytic solution.  Then I plot the difference between the two solutions, which should be scaled by the time factor I input.

Runge[initt_, finalt_, initx_, inity_, 
          steps_] := 
Module[{it = initt, ft = finalt, j, 
            m = steps}, 
h = (ft - it)/
              m; 
T = 
            Table[a + (j - 1) h, {j, 1, m + 1}]; 
R = 
            Table[{initx, inity, initt}, {j, 1, 
                m + 1}]; 
For[j = 1, j < m + 1, j++, 
            k1 = h F[T[[j]], R[[j]]]; 
k2 = 
              h F[T[[j]] + h/2, 
                  R[[j]] + k1/2]; 
k3 = 
              h F[T[[j]] + h/2, 
                  R[[j]] + k2/2]; 
k4 = 
              h F[T[[j]] + h, 
                  R[[j]] + k3]; 
R[[j + 1]] = 
              R[[j]] + 
                1/6 (k1 + 2 k2 + 2 k3 + 
                      k4)]; 
Return[
            R];
];
SetPendParam[] := 
Module[{q = driveforce, 
            p = phase, b = dampen, 
            w = drivefreq}, 
f1[t_, x_, y_, 
              z_] =  -x; 
f2[t_, x_, y_, 
              z_] =  -y; 
f3[t_, x_, y_, z_] =  
            1; 
F[t_, {x_, y_, z_}] = {f1[t, x, y, z], 
              f2[t, x, y, z], 
              f3[t, x, y, z]};
];

SetPendParam[]

testcase  =  Runge[0, 4, 1, 1, 500];

codesol  =  
        Table[{testcase[[i]][[3]], testcase[[i]][[1]]}, {i, 1, 500}];
analyticsol  =  
        Table[{testcase[[i]][[3]], 
            Exp[-testcase[[i]][[3]]]}, {i, 1, 
            500}];
diffsol = 
        Table[{testcase[[i]][[3]], 
            Abs[((Exp[-testcase[[i]][[3]]] - testcase[[i]][[1]])/Exp)[-testcase[[i]][[3]]]]}, {i, 1, 
            500}];

ListPlot[codesol];

ListPlot[analyticsol];

ListPlot[diffsol];

When the derivative values are set to zero, the function output returns a set of points which all agree with the set initial condition.  This was to be expected as if there is no rate of change over any step, that function will remain whereever it was originally.

SetPendParam[] := 
Module[{q = driveforce, 
            p = phase, b = dampen, 
            w = drivefreq}, 
f1[t_, x_, y_, z_] =  
            0; 
f2[t_, x_, y_, z_] =  
            0; 
f3[t_, x_, y_, z_] =  
            0; 
F[t_, {x_, y_, z_}] = {f1[t, x, y, z], 
              f2[t, x, y, z], 
              f3[t, x, y, z]};
];

SetPendParam[]

testcase  =  Runge[0, 4, 1, 1, 500];

codesol  =  
        Table[{testcase[[i]][[3]], testcase[[i]][[1]]}, {i, 1, 500}];

ListPlot[codesol];

If the input parameter for the driving frequency is set to zero, the output is all compiled onto a line in which t = 0, there is a visible range of points over the whole range of 0 ≤ θ ≤ 2π.  This happens because the third derivative function is essentially the rate at which time progresses in the system, so the rate of change of time is equal to zero, which just doesn't work.

SetPendParam[driveforce_, drivefreq_, dampen_, phase_] := 
        Module[{q = driveforce, p = phase, b = dampen, 
            w = drivefreq}, 
f1[t_, x_, y_, z_] = 
            y; 
f2[t_, x_, y_, z_] = 
            q * Cos[z + p] - Sin[x] - b * y; 
f3[t_, 
              x_, y_, z_] = 
            w; 
F[t_, {x_, y_, z_}] = {f1[t, x, y, z], 
              f2[t, x, y, z], f3[t, x, y, z]};];

SetPendParam[0.9, 0, 0.5, 0]

testcase  =  Runge[0, 40, 1, 1, 500];

codesol  =  
        Table[{testcase[[i]][[3]], 
            Mod[testcase[[i]][[1]] - π, 2  π]}, {i, 1, 
            500}];

ListPlot[codesol];

I now realize that I can set the limits on what values are allowed in the input.  So if a value is input that will not work, the function doesn't evaluate.  You must change the input and reevaluate the cell to make it work.

SetPendParam[driveforce_, drivefreq_ /; drivefreq > 0, dampen_, 
          phase_ /; 0 ≤ phase < 2  π] := 
        Module[{q = driveforce, p = phase, b = dampen, 
            w = drivefreq}, 
f1[t_, x_, y_, z_] = 
            y; 
f2[t_, x_, y_, z_] = 
            q * Cos[z + p] - Sin[x] - b * y; 
f3[t_, 
              x_, y_, z_] = 
            w; 
F[t_, {x_, y_, z_}] = {f1[t, x, y, z], 
              f2[t, x, y, z], 
              f3[t, x, y, z]};
];

SetPendParam[ .9, 0.1, 0.5, 2  π]

SetPendParam[ .9, 0.1, 0.5, 2  π - 0.0001]

SetPendParam[ .9, 0, 0.5, 0]

SetPendParam[ .9, 0.1, 0.5, 0]

The first and third commands were not accepted, though the altered versions which fall within the acceptable inputs, were evaluated and the derivative functions can be used.

This little bit I wrote out to show how the error scales with the time steps to ensure that the code was working.  It's probably not the most efficient way to do it, but it shows how the error grows for known analytic solutions for both (dx/dt)= -x and (dy/dt)= y.  The error given by the theory of this algorithm claims to drop in error at a rate of O [(h^4)], which this seems to agree with.

SetPendParam[] := 
        Module[{q = driveforce, p = phase, b = dampen, 
            w = drivefreq}, 
f1[t_, x_, y_, 
              z_] = -x; 
f2[t_, x_, y_, z_] = 
            y; 
f3[t_, x_, y_, z_] = 
            1; 
F[t_, {x_, y_, z_}] = {f1[t, x, y, z], 
              f2[t, x, y, z], 
              f3[t, x, y, z]};
];

SetPendParam[]

\(\(CheckError[
          input_] := 
Module[{i, j, k, l, 
            fin\  = \ input}, 
Do[
            sol\  = \ Runge[0, 3, 1, 1, 10\^i]; 
new\  = \ 
              Table[{\(sol[\([10\^i]\)]\)[\([3]\)], 
                  N[Abs[\(Exp[\(-\(sol[\([10\^i]\)]\)[\([3]\)]\)] - \(sol[\([10\^i]\)]\)[\([1]\)]\)\/Exp[\(-\(sol[\([10\^i]\)]\)[\([3]\)]\)]]]}]; 
Print["The percent error after 10^", 
              i, " steps is ", 
              N[100*new[\([2]\)]], "%, with the value of stepsize equal to ", 10\^i\/3, " steps per unit time."]
, {i, 1, 
              fin}
]; 
Do[
            sol\  = \ Runge[0, 3, 1, 1, 10\^j]; 
new2\  = \ 
              Table[{\(sol[\([j]\)]\)[\([3]\)], 
                  N[Abs[\(Exp[\(-\(sol[\([j]\)]\)[\([3]\)]\)] - \(sol[\([j]\)]\)[\([1]\)]\)\/Exp[\(-\(sol[\([j]\)]\)[\([3]\)]\)]]]}, {j, 1, 
                  10\^j + 1}]; 
ListPlot[new2, 
              PlotLabel -> \ \ "{error, steps/unit t} ="\ \ TraditionalForm[{N[\(Exp[\(-\(sol[\([10\^j]\)]\)[\([3]\)]\)] - \(sol[\([10\^j]\)]\)[\([1]\)]\)\/Exp[\(-\(sol[\([10\^j]\)]\)[\([3]\)]\)]], 
                      10\^j\/3}]];
, {j, 1, 
              fin}
]; 
Do[
            sol\  = \ Runge[0, 3, 1, 1, 10\^k]; 
new3\  = \ 
              Table[{\(sol[\([10\^k]\)]\)[\([3]\)], 
                  N[Abs[\(Exp[\(sol[\([10\^k]\)]\)[\([3]\)]] - \(sol[\([10\^k]\)]\)[\([2]\)]\)\/\(sol[\([10\^k]\)]\)[\([2]\)]]]}]; 
Print["The percent error after 10^", 
              k, " steps is ", 
              N[100*new3[\([2]\)]], "%, with the value of stepsize equal to ", 10\^k\/3, " steps per unit time."]
, {k, 1, 
              fin}
]; 
Do[
            sol\  = \ Runge[0, 3, 1, 1, 10\^l]; 
new4\  = \ 
              Table[{\(sol[\([l]\)]\)[\([3]\)], 
                  N[Abs[\(Exp[\(sol[\([l]\)]\)[\([3]\)]] - \(sol[\([l]\)]\)[\([2]\)]\)\/\(sol[\([l]\)]\)[\([2]\)]]]}, {l, 1, 
                  10\^l + 1}]; 
ListPlot[new4, 
              
              PlotLabel -> \ \ "{error, steps/unit t} ="\ \ TraditionalForm[{N[\(Exp[\(sol[\([10\^l]\)]\)[\([3]\)]] - \(sol[\([10\^l]\)]\)[\([2]\)]\)\/Exp[\(sol[\([10\^l]\)]\)[\([3]\)]]], 
                      10\^l\/3}]];

, {l, 
              1, fin}];
];\)\)

CheckError[3]

The Package

(* : Package Version : 1.0*)
(* : Copyright : Copyright 2007, Jesse Welch*)
(* : Name : Graphics`Pendulum`*)
(* : Author : Jesse Welch*)


Unprotect["Graphics`Pendulum`*"];
ClearAll["Graphics`Pendulum`*"];
BeginPackage["Graphics`Pendulum`"]

Runge::usage = "Runge[initialtime, finaltime, initialx, initialy, numberofsteps] 
 Uses a fourth order runge kutta algorithm to integrate third order ODE's.  The output is in the form of a table of
 3 x steps values.  The first of the three values of the set is the x function, the second value is the y function, 
 and the third value is the time value, for each step";

SetPendParam::usage = "SetPendParam[drivforce,dampen,drivefreq] 
 This sets the values of the parameters of your derivative functions for the ODE solver.  This specifically puts in
 the derivative functions for the chaotic pendulum,though you could make your own function to set the functions and
 parameters of different derivative functions to solve a different problem";

ThetaPlot::usage = "ThetaPlot[{drivforce, drivfreq, drivphase, dampen}, {initt, fint, steps}, {initth, initw}]
 This gives graphical output of theta vs T.";

OmegaPlot::usage = "OmegaPlot[{drivforce, drivfreq, drivphase, dampen}, {initt, fint, steps}, {initth, initw}]
 This gives graphical output of omega vs T.";

Begin["`Private`"]

Runge[initt_, finalt_, initx_, inity_, steps_] := 
Module[{it = initt, ft = finalt, j, m = steps}, 
   h = (ft - it)/m;
   T = Table[a + (j - 1) h, {j, 1, m + 1}];
   R = Table[{initx, inity, initt}, {j, 1, m + 1}];
   For[j = 1, j < m + 1, j++, 
      k1 = h F[T[[j]], R[[j]]];
      k2 = h F[T[[j]] + h/2, R[[j]] + k1/2];
      k3 = h F[T[[j]] + h/2, R[[j]] + k2/2];
      k4 = h F[T[[j]] + h, R[[j]] + k3];
      R[[j + 1]] = R[[j]] + 1/6 (k1 + 2 k2 + 2 k3 + k4)
   ];
   Return[R];
];

SetPendParam[driveforce_, drivefreq_, dampen_, phase_] :=
Module[{q = driveforce, p = phase, b = dampen, w = drivefreq},
   f1[t_, x_, y_, z_] = y;
   f2[t_, x_, y_, z_] = q Cos[z + p] - Sin[x] - b y;
   f3[t_, x_, y_, z_] = w;
   F[t_, {x_, y_, z_}] = {f1[t, x, y, z], f2[t, x, y, z], f3[t, x, y, z]};
];

ThetaPlot[{driveforc_, drivefre_, dampe_, phas_}, {int_, fint_, step_}, {intx_, inty_}] :=
Module[{ddforce = driveforc, ddrivfre = drivefre, ddampe = dampe, pphase = phas, intt = int,
fintt = fint, intxx = intx, intyy = inty, stepp = step,finishh = finish},
   SetPendParam[ddforce, ddrivfre, ddampe, pphase];
   solu = Runge[intt, fintt, intxx, intyy, stepp];
   newsolu2 = Table[{solu[[i]][[3]], Mod[solu[[i]][[1]] - π, 2π]}, {i, 1, stepp}];
   ListPlot[newsolu2];
];

OmegaPlot[{driveforc_, drivefre_, dampe_, phas_}, {int_, fint_, step_}, {intx_, inty_}] :=
Module[{ddforce = driveforc, ddrivfre = drivefre, ddampe = dampe, pphase = phas, intt = int,
fintt = fint, intxx = intx, intyy = inty, stepp = step,finishh = finish},
   SetPendParam[ddforce, ddrivfre, ddampe, pphase];
   solu = Runge[intt, fintt, intxx, intyy, stepp];
   newsolu2 = Table[{solu[[i]][[3]], Mod[solu[[i]][[2]] - π, 2π]}, {i, 1, stepp}];
   ListPlot[newsolu2];
];

End[]
EndPackage[]