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

Template

This document will illustrate what I plan to achieve in this semester.  Basic overview of what functions I plan to have as well as the ones I already have will be shown.

SetPendParam[force, freq, phase, damp]

Description

This function sets the three derivative functions used by Runge.  This function is built specifically to solve the pendulum DE's.

Initialization

This function takes input of four parameters.  Each of these parameters can be identified with one of the pieces of the DE's.  The first three pieces all have to do with the driving force.  You can input any numerical value for the amplitude, frequency, and phase of the driving function.  The last parameter is the strength of the dampening coefficient.

Example

SetPendParam[0.9, 2/3, 0.0, 0.5]

Output

Requirements and Limitations

This function is required to use the Runge function later.  If the frequency parameter is set to zero, time does not evolve, and thus we have no equations of motion.  The phase factor must lie in the range from 0 < phase < 2π.

Function

SetPendParam[driveforce_, drivefreq_ /; drivefreq > 0, 
          phase_ /; 0 ≤ phase < 2  π, dampen_] := 
        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]};
];

SetTimeParam[initialt, finalt, steps]

Description

This function sets the value for the time data in the calculations.  This function will take three values as input, and output a list of three values.  This list can then be used as input in later functions and the time interval does not need to be set for every calculation.

Initialization

The parameters here are very self-explanitory, first there is the initial time and the final time, and then the total number of steps you want to calculate.

Example

SetTimeParam[0, 10, 1000]

Output

Requirements and Limitations

N/A

Function

N/A

Runge[initialt, finalt, intialx, initialy, steps]

Description

This function integrates three derivative functions and outputs the original function as a set of data points in the range 
{initialt ≤ t ≤ finalt} in the form of a table of length steps, and each step is a 3 - vector.

Initialization

This function takes five inputs; initial and final times, initial x and y, and total number of steps.  The z coordinate has an automatic initial value equal to that of the initial time, this is done because in my set of equations, the third function is representative of time.  An example calling of this function is shown below:

Example

Runge[0, 10, 1, 1, 100]

Requirements and Limitations

As I previously mentioned the z coordinate is set as a default to the time coordinate.  This however must match up with

Function

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];
];

ThetaPlot[{dforce,dfreq,dphase,dampe},{inititalt,finalt,steps},{initialθ,initialω}]

Description

This function outputs a listplot of the solutions aquired from the calling of the function Runge.  This function takes two of the three lists output by Runge and sets them in a new list, also taking the Mod of theta over a range of 2π.   This table is then output using the built in ListPlot command, giving theta on the y-axis, and time on the x-axis.

Initialization

There are many inputs into this function, though they are broken into sets.  The first set is exactly what was called in SetPendParam.  These are the parameters for the DE's which are going to be solved.  The second set is exactly what is output by SetTimeParam, these are the initial and final times, as well as the number of total steps.  The last set consists of the initial theta and omega values.

Example

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

Requirements and Limitations

This function requires that the functions Runge and SetPendParam have been defined.  ThetaPlot explicitly calls both of these functions in it's own code.  As this is true, all the limitations of the previous two functions still stand true here.

Function

ThetaPlot[{driveforc_, drivefre_, phas_, dampe_}, {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[{dforce,dfreq,dphase,dampe},{inititalt,finalt,steps},{initialθ,initialω}]

Description

This function outputs a listplot of the solutions aquired from the calling of the function Runge.  This function takes two of the three lists output by Runge and sets them in a new list, also taking the Mod of omega over a range of 2π.  This table is then output using the built in ListPlot command, giving omega on the y-axis, and time on the x-axis.

Initialization

There are many inputs into this function, though they are broken into sets.  The first set is exactly what was called in SetPendParam.  These are the parameters for the DE's which are going to be solved.  The second set is exactly what is output by SetTimeParam, these are the initial and final times, as well as the number of total steps.  The last set consists of the initial theta and omega values.

Example

OmegaPlot[{0.9, 2/3, 0, 0.5}, {0, 10, 1000}, {1, 1}]

Requirements and Limitations

This function requires that the functions Runge and SetPendParam have been defined.  ThetaPlot explicitly calls both of these functions in it's own code.  As this is true, all the limitations of the previous two functions still stand true here.

Function

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];];

PhasePlot[{dforce,dfreq,dphase,dampe},{inititalt,finalt,steps},{initialθ,initialω}]

Description

This function outputs a listplot of the solutions aquired from the calling of the function Runge.  This function takes two of the three lists output by Runge and sets them in a new list, also taking the Mod of theta and omega over a range of 2π.  This table is then output using the built in ListPlot command, giving omega on the x-axis, and theta on the y-axis.

Initialization

There are many inputs into this function, though they are broken into sets.  The first set is exactly what was called in SetPendParam.  These are the parameters for the DE's which are going to be solved.  The second set is exactly what is output by SetTimeParam, these are the initial and final times, as well as the number of total steps.  The last set consists of the initial theta and omega values.

Example

PhasePlot[{0.9, 2/3, 0, 0.5}, {0, 10, 1000}, {1, 1}]

Requirements and Limitations

This function requires that the functions Runge and SetPendParam have been defined.  ThetaPlot explicitly calls both of these functions in it's own code.  As this is true, all the limitations of the previous two functions still stand true here.

Function

PhasePlot[{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[{Mod[solu[[i]][[1]] - π, 2  π], 
                Mod[solu[[i]][[2]] - π, 2  π]}, {i, 1, 
                stepp}]; 
ListPlot[
            newsolu2];];

PoinPlot[{dforce,dfreq,dphase,dampe},{inititalt,finalt,steps},{initialθ,initialω}]

Description

Initialization

Example

Requirements and Limitations

Function

CalcLyapunov[]

Description

Initialization

Example

Requirements and Limitations

Function

IsChaotic[]

Description

Initialization

Example

Requirements and Limitations

Function

Bifurcation[]

Description

Initialization

Example

Requirements and Limitations

Function

MapParamSpace[]

Description

Initialization

Example

Requirements and Limitations

Function

Description

Initialization

Example

Requirements and Limitations

Function

Functions Not Here Yet

Bifurcation

I already have a bifurcation function which uses the logistic equation to produce a nice map of bifurcation.  I believe that I can take the same type of idea and replace the function equation to solve this set of DE's and produce another bifurcation diagram.

CalcLyapunov

Solving DE's

This will use SetPendParam, SetTimeParam, and Runge.  Output of equations of motion into a stored table.  This is and can be done already.

Comparing Solutions to similar DE's

Once the DE's can be solved, we can take identical parameter sets, and arbitrarily close initial conditions and compare the set of two solutions.  We can then make a table of the differences in the solution as a function of time, from this the LE can be calculated.  To do this we will need to somehow fit an exponential curve to out table.  Not sure how to do that yet.

IsChaotic

Solving DE's

This will use SetPendParam, SetTimeParam, and Runge.  Output of equations of motion into a stored table.  This is and can be done already.

Is the solution Chaotic?

This will essentially be just a test on the value of the LE, if positive this function will note the parameter set used, and mark off the corresponding value in another table, this table corresponding to each parameter set value over as fine a gradiant as is given.

MapParamSpace

Solving DE's

This will use SetPendParam, SetTimeParam, and Runge.  Output of equations of motion into a stored table.  This is and can be done already.

Is the solution Chaotic?

This will require calculating the Lyapunov exponents (LE)  for some range of  3 dimensional parameter space.  Then calculate the average LE over a large average of solutions of the same parameter set, with different initial conditions.  Then if the average LE is positive for some parameter set, we make a note of where we are in parameter space and mark that point as chaotic, and mark that point in a different way if the LE is negative or zero

There is an issue with this still, as the notion of transience still makes me a little uneasy when dealing with the long term behavior of trajectories.  The fact that the solutions can only get so far from each other, the beginning of the divergence is going to be the most important piece of calculating the LE.  The transience however is only affecting the solution only in the first short amount of time.  So these two things mask their individual results together into one result that seems impossible to peel them apart.

To do this effectivley, I will need to calculate a very large number of solutions such that taking an average will hopfully cancel the transience behavior and leave me with an actual average LE that does match up to the actual reality of the solution.

Plotting the Parameter Space

Once we have data which describes the parameter space, ie whether or not solutions are chaotic, we need to plot that out somehow.  Either a 2d density or contour map, or some kind of 3d imaging.  Not sure what will be easiest or most visually appealing