PHY 410-505 Home Home  |  Course Outline  |  Lectures  |  Homework  |  Files

Lecture 23: October 27


Poincaré Maps and Poincaré Sections

We have seen that the Poincaré Map, defined on page 157, provides a much simpler picture of the period-doubling and chaotic motion of the forced damped pendulum, for example, that a plot of the trajectory of the system in phase space. This map is essentially a stroboscopic picture of the phase-space trajectory with strobe period equal to the period of the driving force. Another example is the kicked rotor, in which the Poincaré Map at the kicking frequency is equivalent to the Standard Map.

There are many interesting systems, however, in which there is no convenient strobe frequency at which to examine a trajectory in phase space. An example is the double pendulum discussed on pages 166-169. This is a time-independent Hamiltonian system with a 4-dimensional phase space. It exhibits chaotic motion for some values of the energy E and initial conditions. Poincaré introduced an extremely useful concept called the Poincaré Section to examine the behavior of systems such as the double pendulum for which there is no obvious natural or driving period. A Poincaré section is simply a slice of the phase space. Every time the trajectory crosses this slice in a particular direction, we plot the phase point projected on some convenient two-dimensional plane. It turns out that the section points of periodic or regular trajectories lie on one-dimensional curves, but section points of chaotic trajectories spread out in the two-dimensional plane.

The Double Pendulum

The following applet is based generated by the program DoublePendulum.java which follows the suggestions on pages 168-169.

This applet animates one or two double pendulums in the left panel, and plots the Poincaré section in the right panel. Section points are generated when the second mass swings through the vertically downward direction in the counter-clockwise direction. At this instant, the angular position and generalized momentum of the first mass are recorded on the section plot.

Here is the code for this applet:

// DoublePendulum.java

import comphys.*;
import comphys.graphics.*;

public class DoublePendulum extends Animation {

    double t;
    double dt = 0.003;

    class Pendulum {
        double L = 1;
        double m = 1;
        double g = 9.8;

        String name;
        double E;
        double[] Y0 = new double[4];
        double[] Y = new double[4];
        double[] YSection = new double[4];
        double q2Old;           // for figuring section crossing
        int count;              // section point number

        Pendulum (double E, double q1, double q2, String name) {
            this.E = E;
            Y0[0] = q1;
            Y0[1] = q2;
            reset();
            this.name = name;
        }

        void reset () {
            double q1 = Y0[0];
            double q2 = Y0[1];
            // set p1 = 0 and attempt to solve for p2
            double p1 = Y0[2] = 0;
            double p2 = E - m * g * L * (3 - 2 * Math.cos(q1) - Math.cos(q2));
            double s = Math.sin(q1 - q2);
            p2 *= m * L * L * (1 + s * s);
            if (p2 < 0) {
                System.out.println(name + ": E = " + Easy.format(E, 6) +
                                   " too small, setting q1 = q2 = 0");
                Y0[0] = Y0[1] = Y0[2] = 0;
                Y0[3] = Math.sqrt(E);
            } else {
                Y0[3] = Math.sqrt(p2);
            }
            for (int i = 0; i < 4; i++)
                YSection[i] = Y[i] = Y0[i];
            q2Old = Y[1];
            count = 0;
        }

        boolean findSection;

        void NewtonsEquations (double[] Y, double[] dYdt) {
            double q1 = Y[0];
            double q2 = Y[1];
            double p1 = Y[2];
            double p2 = Y[3];
            
            double s = Math.sin(q1 - q2);
            double c = Math.cos(q1 - q2);
            double d = 1 + s * s;
            double a = (p1 * p1 + 2 * p2 * p2) * c - p1 * p2 * (2 + c * c);
            a *= s / (m * L * L * d * d);
        
            dYdt[0] = (p1 - p2 * c) / (m * L * L * d);
            dYdt[1] = (2 * p2 - p1 * c) / (m * L * L * d);
            dYdt[2] = a - 2 * m * g * L * Math.sin(q1);
            dYdt[3] = - a - m * g * L * Math.sin(q2);

            if (findSection) {
                double denom = dYdt[1];
                for (int i = 0; i < 4; i++)
                    dYdt[i] /= denom;
            }
        }

        double energy () {
            double q1 = Y[0];
            double q2 = Y[1];
            double p1 = Y[2];
            double p2 = Y[3];
        
            double s = Math.sin(q1 - q2);
            double c = Math.cos(q1 - q2);
            double c1 = Math.cos(q1);
            double c2 = Math.cos(q2);
                
            return 1 / (2*m*L*L) * (p1*p1 + 2*p2*p2 - 2*p1*p2*c)
                / (1 + s*s) + m*g*L * (3 - 2*c1 - c2);
        }

        double[] k1 = new double[4];
        double[] k2 = new double[4];
        double[] k3 = new double[4];
        double[] k4 = new double[4];
        double[] YTemp = new double[4];
    
        void RungeKutta4 (double dt) {
            NewtonsEquations(Y, k1);
            for (int i = 0; i < 4; i++)
                YTemp[i] = Y[i] + dt * k1[i] / 2;
            NewtonsEquations(YTemp, k2);
            for (int i = 0; i < 4; i++)
                YTemp[i] = Y[i] + dt * k2[i] / 2;
            NewtonsEquations(YTemp, k3);
            for (int i = 0; i < 4; i++)
                YTemp[i] = Y[i] + dt * k3[i];
            NewtonsEquations(YTemp, k4);
            for (int i = 0; i < 4; i++)
                Y[i] += (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) * dt / 6;
        }

        void compactify () {
            // bring q1, q2 in range -pi, pi
            while (Y[0] > Math.PI)
                Y[0] -= 2 * Math.PI;
            while (Y[0] < -Math.PI)
                Y[0] += 2 * Math.PI;
            while (Y[1] > Math.PI)
                Y[1] -= 2 * Math.PI;
            while (Y[1] < -Math.PI)
                Y[1] += 2 * Math.PI;
        }

        void saveQ2 () {
            q2Old = Y[1];
        }

        double[] YSave = new double[4];
        boolean newIntersection;
        
        boolean intersection () {
            if ( !(Y[3] > 0 && q2Old * Y[1] < 0) )
                return newIntersection = false;
            findSection = true;
            for (int i = 0; i < 4; i++)
                YSave[i] = Y[i];
            RungeKutta4(-Y[1]);
            for (int i = 0; i < 4; i++) {
                YSection[i] = Y[i];
                Y[i] = YSave[i];
            }
            findSection = false;
            ++count;
            return newIntersection = true;
        }
    }

    Pendulum pendulum1 = new Pendulum(15, 0, 0, "pendulum 1");
    Pendulum pendulum2 = new Pendulum(15, 1.1, 0, "pendulum 2");
    boolean twoPendula;

    class RealSpace extends Plot {
        RealSpace () {
            setSize(300, 300);
            setBackground("white");
        }

        void drawPendulum (double q1, double q2, String color) {
            double x1 = Math.sin(q1);
            double x2 = x1 + Math.sin(q2);
            double y1 = - Math.cos(q1);
            double y2 = y1 - Math.cos(q2);
            setColor("green");
            plotLine(0, 0, x1, y1);
            plotLine(x1, y1, x2, y2);
            setColor(color); 
            double r = 0.1;
            floodCircle(x1 - r, x1 + r, y1 - r, y1 + r);
            floodCircle(x2 - r, x2 + r, y2 - r, y2 + r);
        }

        public void paint () {
            clear();
            setWindow(-2.15, 2.15, -2.15, 2.15);
            drawPendulum(pendulum1.Y[0], pendulum1.Y[1], "red");
            if (twoPendula)
                drawPendulum(pendulum2.Y[0], pendulum2.Y[1], "blue");
            double r = 0.05;
            setColor("black");
            floodCircle(-r, r, -r , r);

            plotOverlay();
            setColor("black");
            plotStringCenter("t = " + Easy.format(t, 6), 0, 1.9);
            setColor("red");
            plotString("E = " + Easy.format(pendulum1.energy(), 10),
                       -1.9, -2.05);
            if (twoPendula) {
                setColor("blue");
                plotStringLeft("E = " + Easy.format(pendulum2.energy(), 10),
                               1.9, -2.05);
            }
        }
    }

    double qMin = -Math.PI / 2;
    double qMax = Math.PI / 2;
    double pMin = -6;
    double pMax = 8;

    class Section extends Plot {
        Section () {
            setSize(300, 300);
            setBackground("white");
        }

        boolean needToClear = true;
        
        public void paint () {
            if (needToClear) {
                clear();
                setColor("gray");
                drawAxes(qMin, qMax, pMin, pMax);
                needToClear = false;
            }
            if (pendulum1.newIntersection) {
                setColor("red");
                plotPoint(pendulum1.YSection[0], pendulum1.YSection[2]);
            }
            if (pendulum2.newIntersection) {
                setColor("blue");
                plotPoint(pendulum2.YSection[0], pendulum2.YSection[2]);
            }

            plotOverlay();
            double dq = (qMax - qMin) * 0.02;
            double dp = (pMax - pMin) * 0.02;
            setColor("magenta");
            plotStringCenter("Poincare Section q2 = 0, p2 > 0",
                             0, pMax + dq);
            setColor("red");
            double q = pendulum1.YSection[0];
            double p = pendulum1.YSection[2];
            plotBox(q - dq, q + dq, p - dp, p + dp);
            plotString("Crossing No. " + pendulum1.count,
                       qMin, pMin - 3 * dp);

            if (twoPendula) {
                setColor("blue");
                q = pendulum2.YSection[0];
                p = pendulum2.YSection[2];
                plotBox(q - dq, q + dq, p - dp, p + dp);
                plotStringLeft("Crossing No. " + pendulum2.count,
                               qMax, pMin - 3 * dp);
            }
        }
    }

    RealSpace realSpace;
    Section section;
    int skip = 10;
    Reader E1Reader, E2Reader, dtReader, skipReader;
    Reader p1q1Reader, p1q2Reader, p2q1Reader, p2q2Reader;
    Reader twoReader;
    Reader qMinReader, qMaxReader, pMinReader, pMaxReader;

    public void init () {
        add(realSpace = new RealSpace());
        add(section = new Section());
        add(E1Reader = new Reader("Red E = ", pendulum1.E));
        add(p1q1Reader = new Reader("Red q1 = ", pendulum1.Y0[0]));
        add(p1q2Reader = new Reader("Red q2 = ", pendulum1.Y0[1]));
        add(twoReader = new Reader("Blue pendulum?", twoPendula));
        add(E2Reader = new Reader("Blue E = ", pendulum2.E));
        add(p2q1Reader = new Reader("Blue q1 = ", pendulum2.Y0[0]));
        add(p2q2Reader = new Reader("Blue q2 = ", pendulum2.Y0[1]));
        add(dtReader = new Reader("dt =", dt));
        add(skipReader = new Reader("skip steps", skip));
        add(qMinReader = new Reader("qMin = ", qMin, 4));
        add(qMaxReader = new Reader("qMax = ", qMax, 4));
        add(pMinReader = new Reader("pMin = ", pMin, 4));
        add(pMaxReader = new Reader("pMax = ", pMax, 4));
        addControlPanel();
    }

    boolean updateSection;
    
    void takeStep () {
        t += dt;
        updateSection = false;
        boolean yes = false;
        pendulum1.saveQ2();
        pendulum1.RungeKutta4(dt);
        pendulum1.compactify();
        if (pendulum1.intersection())
            updateSection = true;
        if (!twoPendula)
            return;
        pendulum2.saveQ2();
        pendulum2.RungeKutta4(dt);
        pendulum2.compactify();
        if (pendulum2.intersection())
            updateSection = true;
    }

    public void step () {
        int s = skip;
        do {
            takeStep();
            if (updateSection)
                break;
        } while (--s > 0);
        if (updateSection)
            section.repaint();
        realSpace.repaint();
    }

    public void reset () {
        pendulum1.E = E1Reader.readDouble();
        pendulum1.Y0[0] = p1q1Reader.readDouble();
        pendulum1.Y0[1] = p1q2Reader.readDouble();
        twoPendula = twoReader.readBoolean();
        pendulum2.E = E2Reader.readDouble();
        pendulum2.Y0[0] = p2q1Reader.readDouble();
        pendulum2.Y0[1] = p2q2Reader.readDouble();
        dt = dtReader.readDouble();
        skip = skipReader.readInt();
        qMin = qMinReader.readDouble();
        qMax = qMaxReader.readDouble();
        pMin = pMinReader.readDouble();
        pMax = pMaxReader.readDouble();
        t = 0;
        pendulum1.reset();
        pendulum2.reset();
        realSpace.repaint();
        section.needToClear = true;
        section.repaint();
    }

    public static void main (String[] args) {
        DoublePendulum doublePendulum = new DoublePendulum();
        doublePendulum.frame("Double Pendulum", 630, 550);
    }
}
This program uses the fourth order Runge-Kutta algorithm with time step dt which can be set by the user. The program monitors section crossings, and then uses a simple change of variable from t to q2 to estimate the section point precisely.


UB Physics Home Questions or comments: phygons@acsu.buffalo.edu