About this worksheet

Equation of motion. Define the system

Newton's laws immediately yield the following equation of motion for a pendulum
with a point mass m at the end of a massless stiff rod of length L.  measures the
angle at the suspension point between the rod and the vertical.

>	EQM:=m*L*diff(theta(t),t,t)=-m*g*sin(theta(t)): EQM;

For the sake of the following investigations choose suitable units such that g=L.
To write the 2nd order differential equation as a system of first order DEs intro-
duce the state variables     and .
The followiong forces MAPLE to inertly pretty-print the equations......

>	SYS:=array([[diff(y[1](t),t)=y[2]],             [diff(y[2](t),t)=-sin(y[1](t))]]): print(SYS);

For later calculations conveniently introduce the two functions -- i.e. if changes
are desired (e.g. diffferent lengths, or including damping) changes need to be made
only here:

>	f1:=(y1,y2)->y2; f2:=(y1,y2)->-sin(y1);

The phase-plane portrait

First take a look at the phase-plane portrait on both a small and a large scale.

>	with(plots);

Warning, the name changecoords has been redefined

>	setoptions3d(scaling=constrained);

>	sL1:=fieldplot([f1,f2],-2..2,-2..2,      scaling=constrained,grid=[13,13],      color=blue,thickness=2,      labels=[`theta`,`omega`],      title=`Simple pendulum: Small angles`): sL1;

Note that the solution curves are NOT circles, but are slightly distorted. Take
a look at a few overlaid circles:

>	display([sL1,seq(plot([r*cos(t),r*sin(t),t=0..2*Pi],                      thickness=1,color=magenta),                r=[1/2,1,3/2,2,5/2])],         view=[-2..2,-2..2]);

A view on a larger scale even more clearly shows the differences:

>	sL2:=fieldplot([f1,f2],-2..13,-5..5,      scaling=constrained,grid=[15,15],      color=blue,thickness=2,      labels=[`theta`,`omega`],      title=`Simple pendulum: Large angles`): sL2;

Numerical integration. "Euler's method performs poorly."

Define the basic parameters: initial conditions, length of time interval, step-size::

>	t0:=0; tfin:=7; N:=20; dt:=(tfin-t0)/N; y10:=1; y20:=0;

Use Euler's method to numerically integrate the system of differential equations.
The following implementation first creates empty tables to reserve contiguous
memory space to later hold the calculated values.

>	TT:=array([t0,seq(t0+k*dt,k=1..N)]): Y1:=array([y10,seq(0,k=1..N)]): Y2:=array([y20,seq(0,k=1..N)]):

By default MAPLE will use exact expressions, like sin(23/100) -- the result will
be huge -- instead we want decimal approximations. In order to insure that
round-off errors play NO significant rule we work with unusually high accuracies.

>	for k from 1 to N do     Y1[k+1]:=evalf(Y1[k]+dt*f1(Y1[k],Y2[k]),25):     Y2[k+1]:=evalf(Y2[k]+dt*f2(Y1[k],Y2[k]),25):     od:

In order to plot the points it is most elegant to "zip" together the arrays of points.
For practical purposes it is convenient to convert the internal data format to lists
and to round to a few digit -- after all the calcultions are done, and for plotting
one does not need such high accuracy::

>	TT:=evalf(convert(TT,list),4); Y1:=evalf(convert(Y1,list),4); Y2:=evalf(convert(Y2,list),4);

>	display([    plot(zip((a,b)->[a,b],TT,Y1),color=blue),    plot(zip((a,b)->[a,b],TT,Y2),color=red)]);

While this looks quite nice, a different plot is more disturbing:

>	curve:=plot(zip((a,b)->[a,b],Y1,Y2),      scaling=constrained,color=magenta,thickness=2,      labels=[`theta`,`omega`],      tickmarks=[4,4]): curve;

>	display([curve,sL1]);

Try it with different initial data, different time intervals, different step-sizs (change
the value of the number N of steps used).
Instead of scrolling up and re-executing the above commands it is opportune to
bundle them in a short procedure:

>	op(1,3..4);

>

pendeuler:=proc(F,y10,y20,trange,N)    local f1, f2, t0, tfin, dt, k,          Y1, Y2, TT, curve, pp;    f1:=F[1];    f2:=F[2];    t0:=op(1,trange);    tfin:=op(2,trange);    dt:=(tfin-t0)/N; TT:=array([t0,seq(t0+k*dt,k=1..N)]): Y1:=array([y10,seq(0,k=1..N)]): Y2:=array([y20,seq(0,k=1..N)]): for k from 1 to N do     Y1[k+1]:=evalf(Y1[k]+dt*f1(Y1[k],Y2[k]),25):     Y2[k+1]:=evalf(Y2[k]+dt*f2(Y1[k],Y2[k]),25):     od: TT:=evalf(convert(TT,list),4); Y1:=evalf(convert(Y1,list),4); Y2:=evalf(convert(Y2,list),4); curve:=plot(zip((a,b)->[a,b],Y1,Y2),              color=magenta,thickness=2):   pp:=fieldplot(F,           min(op(Y1))..max(op(Y1)),           min(op(Y2))..max(op(Y2)),           grid=[13,13],           color=blue,           thickness=1):            display([curve,pp,      textplot([0.6*max(op(Y1)),1.09*max(op(Y2)),                cat("N=",N,",   dt=",convert(evalf(dt,3),string))])],      labels=[`theta`,`omega`],      tickmarks=[3,3],      scaling=constrained,      title=`Euler's method performs very poorly`); end:

>	pendeuler([f1,f2],1,0,0..7,70);

Then try a larger number of steps, and an intial condition closer to the rest-point::

>	pendeuler([f1,f2],1/2,0,0..7,350);

Even with larger numbers of steps it still spirals!

Think about what this means for the total energy in the system -- e.g.
consider the points where the purported solution curve crosses the
coordinate axes.
Euler's method provides a very crude approximation -- but this simple
example shows that there is NO HOPE for obtaining closed solution
curves -- the solution curves obtained by Euler's method for this system
WILL ALWAYS SPIRAL OUTWARDS.  WHY????

The explanation deserves a nicely written paragraph...

Continuing for a larger time interval with the same step size dt:
pendeuler([f1,f2],1/2,0,0..40,2000);

>	pendeuler([f1,f2],1,0,0..22,110);

Morale:
Euler's method gives an intuitive understanding about what numerical algorithms
for solving differential equations do -- but in the end one definitively wants more
accurate solution algorithms.
Practically any software package has a whole battery of such algorithms for diffe-
rent purposes. In MAPLE load the DEtools package, and browse the help-pages.

>

Using professional solution algorithms

>	with(DEtools);

>	?phaseportrait

>	?dsolve[classical]

>

phaseportrait([D(theta)(t)=omega(t),                D(omega)(t)=-sin(theta(t))],               [theta(t),omega(t)],                t=0..7,               [[theta(0)=1/2,omega(0)=0]],                stepsize=.1,                color=blue,                linecolour=sin(t),                method=classical[foreuler],                scaling=constrained,                title=`Classical forward Euler's method with dt=0.1`);

>

phaseportrait([D(theta)(t)=omega(t),                D(omega)(t)=-sin(theta(t))],               [theta(t),omega(t)],                t=0..7,               [[theta(0)=1/2,omega(0)=0]],                stepsize=.1,                color=blue,                linecolour=sin(t),                method=rkf45,                scaling=constrained,                title=`Default method RKF45 with dt=0.1`);

>

phaseportrait([D(theta)(t)=omega(t),                D(omega)(t)=-sin(theta(t))],               [theta(t),omega(t)],                t=0..100,               [[theta(0)=1/2,omega(0)=0]],                stepsize=.1,                color=blue,                linecolour=sin(t),                method=rkf45,                scaling=constrained,                title=`Default method RKF45 over [0,100] with dt=0.1`);

There are many alternative tools provided in the DEtools package. E.g. note
that phaseportrait draws all arrows at fixed length, whereas fieldplot leaves
their relative sizes as defined.
There are several other different tools for plotting solutions of differential equns,
and dsolve itself allows for many different options. It takes quite a while to become
familiar with the versatility of these tools.

>	with(DEtools);

Animations of the simple pendulum

>	restart; with(plots): with(plottools): with(DEtools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Warning, the name translate has been redefined

>

animcurve:=proc(y10,y20,tfin,N) local sol, mydata, NN,       y1min, y1max, y2min, y2max; sol:=dsolve({D(theta)(t)=omega(t),              D(omega)(t)=-sin(theta(t)),              theta(0)=y10,omega(0)=y20},              [theta(t),omega(t)],type=numeric,method=rkf45,              value=array([seq(k*tfin/N,k=0..N)])); mydata:=sol[2,1]: NN:=max(op(map(q->q[1],[indices(mydata)]))): y1min:=min(seq(mydata[k,2],k=1..NN)): y1max:=max(seq(mydata[k,2],k=1..NN)): y2max:=max(seq(mydata[k,3],k=1..NN)): y2min:=min(-y2max,seq(mydata[k,3],k=1..NN)): y2max:=max(y2max,-y2min): display([seq(    display([fieldplot([y2,-sin(y1)],                        y1=y1min..y1max,y2=y2min..y2max,                        grid=[13,13],thickness=1,color=blue),             plot([seq([mydata[k,2],mydata[k,3]],k=1..nf)],                  color=red,thickness=2)             ]),             nf=1..NN)],insequence=true,          labels=[`theta`,`omega`],      title=`Animation of solution curve via RKF45`,      scaling=constrained); end:

>

animpend:=proc(y10,y20,tfin,N) local NN, sol, mydata; sol:=dsolve({D(theta)(t)=omega(t),              D(omega)(t)=-sin(theta(t)),              theta(0)=y10,omega(0)=y20},              [theta(t),omega(t)],type=numeric,method=rkf45,              value=array([seq(k*tfin/N,k=0..N)])); mydata:=sol[2,1]: NN:=max(op(map(q->q[1],[indices(mydata)]))): display([seq(    display([plot([[0,0],[0.9*sin(mydata[k,2]),-0.9*cos(mydata[k,2])]],                   thickness=2,color=black),             disk([sin(mydata[k,2]),-cos(mydata[k,2])],0.1,color=red),             arrow([sin(mydata[k,2]),-cos(mydata[k,2])],                   [sin(mydata[k,2])+(0.0001+mydata[k,3]/3*cos(mydata[k,2])),                   -cos(mydata[k,2])+(0.0001+mydata[k,3]/3*sin(mydata[k,2]))],                   0.02,0.05,0.2,color=green)             ]),             k=1..NN)],insequence=true,      axes=NONE,      title=`Animation of simple pendulum via RKF45`,      scaling=constrained); end:

>	animcurve(Pi/6,0,6.4,18); animpend(Pi/6,0,6.3,18);

Warning, the 'value' option is deprecated, use 'output' instead

Warning, the 'value' option is deprecated, use 'output' instead

>	animcurve(0,1.998,18,50); animpend(0,1.998,17.5,50);

Warning, the 'value' option is deprecated, use 'output' instead

Warning, the 'value' option is deprecated, use 'output' instead

>	animcurve(0,2.002,12,50);

Warning, the 'value' option is deprecated, use 'output' instead

>	animpend(0,2.002,8.75,50);

Warning, the 'value' option is deprecated, use 'output' instead

>

>