Poincaré maps

for Mathematica versions 5.2 or earlier

By Gianluca Gorni, Department of Mathematics and Computer Science, University of Udine, Italy. Latest revision: June 16, 2009.
Thanks to Murray Eisenberg for his suggestions.

Description: The Poincaré map at time t associated with a (vector) differential equation x' (t) = f(x(t), t) is the map that associates to each starting point (Subscript[x, 0]) the state of the associated solution at the time t. The Poincaré map is important for example in the study of long-term behaviour of solutions, the existence of periodic solutions, and chaos. The package contained in this notebook file defines the following functions:

PoincareMap and PoincareMapPlot to help visualize the Poincaré map of a two-dimensional system;

PoincareMapFixedPointDensityPlot, PoincareMapFixedPointPlot3D and FindPoincareMapFixedPoint assist in finding fixed points of the Poincaré map.

The package is embedded in this notebook file. You load the package by clicking on the following button (whose effect is just to "Evaluate Initialization"), and then you can evaluate the tutorial examples, and make your own variations. Read on for instruction on how to make a more permanent installation.

Button[
      StyleForm["Click here to load package!", 
        FontColor -> RGBColor[1, 0, 0], FontSize -> 16, 
        FontWeight -> "Bold"], 
      ButtonFunction :> 
        FrontEndExecute[FrontEndToken["EvaluateInitialization"]], 
      ButtonEvaluator -> None, Active -> True]

PoincareMaps was formerly part of the ODEPlot package, now superseded by CurveGraphics.

I wrote the package for myself and for anybody interested. I make no guarantee of any kind. A version of the package, called PoincareMaps6, is available for Mathematica version 6 or later.

Instructions: load the package by clicking on the button below, then evaluate the examples. This notebook is not meant for absolute beginners of Mathematica. It assumes some familiarity with notebooks and with the syntax of basic standard Mathematica commands such as Plot and ParametricPlot.

If you wish to be able to load this package from other notebooks, you can generate a package file by resetting the notebook option  AutoGeneratedPackage to Automatic. The next time you save this notebook, Mathematica will generate a new package file called PoincareMaps.m inside the same directory as this notebook file. Then you should place this new PoincareMaps.m file somewhere Mathematica will find it, for example in the directory that you will get from evaluating the following command:

ToFileName[{$UserAddOnsDirectory, "Applications"}]

Mathematica package

This is the program. It is meant to be read only by Mathematica programmers.

Preliminaries and usage documentation

Off[General::"spell"]

Off[General::"spell1"]

\(\(BeginPackage["PoincareMaps`", "Utilities`FilterOptions`", "Graphics`Arrow`", "Graphics`FilledPlot`", "Graphics`PlotField`", "DrawGraphics`DrawingMaster`"];\)\)

Unprotect["`*"];
ClearAll["`*"];

\(\(PoincareMap::usage = 
"PoincareMap[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},options] calculates the Poincare map at time t1 starting at time t0 of the differential system x'=fx, y'=fy on on the region uMin<u<uMax, vMin<v<vMax, where u,v are parameters, and the default parameterization is cartesian (i.e., u=x, v=y). If you want (u,v) to be the polar parameters (r,theta) use the option GridParameterization->Polar. You can specify as GridParameterization any function that takes two variables and gives back a 2-element list, like GridParameterization->Function[{u,v},{u^2,v+u}]. PoincareMap will numerically solve the differential system in a rectangular grid of (u,v) points. The number of points in the two directions is governed by the option GridPoints (that behaves like PlotPoints in ParametricPlot3D), that defaults to 5 in each direction. The Poincaré map in the rest of the region is obtained by ListInterpolation.
Example: PoincareMap[{{x,y},{y,-2x}},{t,0,Pi},{{0,1},{-Pi,Pi}},GridParameterization->Polar].";\)\)

\(\(PoincareMapPlot::usage = 
"PoincareMapPlot[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},options] makes a plot of PoincareMap. The syntax for the non-optional arguments is the same as for PoincareMap. The style of the grid lines is governed by the option GridStyle. Normally you get two plots side-by-side: on the left there is the map at time t0 and on the right the map at time t1. With OnlyEndGrid->True you only get the map at time t1. With StackVertical->True you get the two plots stacked vertically instead of horizontally. With Superimposed->True you get the two plots superimposed.
Example: PoincareMapPlot[{{x,y},{y,-2x}},{t,0,Pi},{{0,1},{-Pi,Pi}},GridParameterization->Polar].";\)\)

\(\(PoincareMapFixedPointDensityPlot::usage = 
"PoincareMapFixedPointDensityPlot[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},epsilon,options] helps in visually searching for fixed points in the PoincareMap. It calculates first the PoincareMap phi, then the function d(x,y)=|phi(x,y)-(x,y)|^2, and makes a DensityPlot of Min[d(x,y),epsilon]. This way the regions where d(x,y) is less than epsilon will be highlighted in grey-black, while the rest of the plot will be white. The plot will be parameterized with (u,v), not with (x,y).
Example: PoincareMapFixedPointDensityPlot[{{x,y},{y,-(x/(1+x^2)+3x(x^2/(1+x^2)))}},{t,0,Pi},{{0,1},{-Pi,Pi}},1/100,GridParameterization->Polar].";\)\)

\(\(PoincareMapFixedPointPlot3D::usage = 
"PoincareMapFixedPointPlot3D[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},{u,u0},{v,v0},options] is similar to PoincareMapFixedPointDensityPlotm, but it makes a ParametricPlot3D instead of a DensityPlot. You will have to adjust PlotRange and BoxRatios to see the minima.
Example: PoincareMapFixedPointPlot3D[{{x,y},{y,-(x/(1+x^2)+3x(x^2/(1+x^2)))}},{t,0,Pi},{{0,1},{-Pi,Pi}},1/100,GridParameterization->Polar].";\)\)

\(\(FindPoincareMapFixedPoint::usage = 
"FindPoincareMapFixedPoint[{{x,y},{fx,fy}},{t,t0,t1},{{uMin,uMax},{vMin,vMax}},{u,u0},{v,v0},options] applies a FindMinimum to PoincareMap with starting point (u0,v0).
Example: FindPoincareMapFixedPoint[{{x,y},{y,-(x/(1+x^2)+3x(x^2/(1+x^2)))}},{t,0,5},{{.65,.85},{-.1,.1}},{r,.76},{θ,0},GridPoints->{40,11},GridParameterization->Polar].";\)\)

GridPoints::usage = 
"GridPoints is an option for \
PhasePlot, PoincareMap, PoincareMapPlot, PoincareMapFixedPointDensityPlot, \
PoincareMapFixedPointPlot3D and FindPoincareMapFixedPoint. It specifies how \
many divisions there will be in the u,v directions, just like PlotPoints does \
for ParametricPlot3D. The default is 5 in each direction.";

\(\(GridParameterization::usage = 
"GridParameterization is an option for PoincareMap, PoincareMapPlot, PoincareMapFixedPointDensityPlot, PoincareMapFixedPointPlot3D and FindPoincareMapFixedPoint. It specifies how the parameters (u,v) are mapped to (x,y). The default is the identity, that is, ({##}&). A predefined value is GridParameterization->Polar, which is shorthand for GridParameterization->Function[{r,theta},r{Cos[theta],Sin[theta]}]. In general as a value for GridParameterizationyou can give any object c that when in the pattern c[x,y], with numeric x,y, gives a list of two numbers {x,y}.";\)\)

Polar::usage = 
"Polar is a possible value for \
the option GridParameterization to be used for PoincareMap, PoincareMapPlot, \
PoincareMapFixedPointDensityPlot, PoincareMapFixedPointPlot3D and \
FindPoincareMapFixedPoint. It is shorthand for \
Function[{r,theta},r{Cos[theta],Sin[theta]}].";

GridStyle::usage = 
"GridStyle is an option for \
PoincareMapPlot and PoincareMapFixedPointDensityPlot. It specifies the style \
of the grid lines. The default is a hue from 0 to 1/2 in one direction and \
from 1/2 to 1 in the other direction. With GridStyle->None all lines will be \
black. To make your own style, just make a list {styleU,styleV} where styleU \
and styleV are each in the format required by ParametricPlot.";

OnlyEndGrid::usage = 
"OnlyEndGrid is an option \
for PoincareMapPlot. With the default OnlyEndGrid->False you will get two \
grids, one at the initial and one at the final time. With OnlyEndGrid->True \
you only get the grid at the final time.";

Superimposed::usage = 
"Superimposed is an option \
for PoincareMapPlot. With the default Superimposed->False you will get two \
grids, one at the initial and one at the final time in separate plots. With \
Superimposed->True you only get the two grid superimposed in one \
plot.";

StackVertical::usage = 
"StackVertical is an \
option for PoincareMapPlot. With the default StackVertical->False you will \
get two grids, one at the initial and one at the final time, side by side. \
With StackVertical->True you get the two grid one below the other.";

SamePlotRange::usage = 
"SamePlotRange is an \
option for PoincareMapPlot. With the default SamePlotRange->True you will get \
two grids, one at the initial and one at the final time, with the same \
PlotRange. With SamePlotRange->False you get Automatic settings for the \
PlotRange in each plot separately.";

tail

Begin["`Private`"];

Poincaré map and fixed point thereof

Default options and error messages

Options[PoincareMap] = {GridPoints -> 5, 
          GridParameterization -> ({##} &)};

Options[PoincareMapPlot] = {GridStyle -> Automatic, 
          OnlyEndGrid -> False, Superimposed -> False, 
          StackVertical -> False, 
          SamePlotRange -> True};

PoincareMap::"BadGridPoints" = "Bad GridPoints"; 
    PoincareMap::"BadParameterization" = "Bad \
GridParameterization";

This utility processes the grid data given by the user into a form more readily usable by the other functions.

\(\(baseGrid[{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
            GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
            GridVMin < gridVMax, 
          opts___?OptionQ] := 
Module[{hm, grdPrmtrzt, 
            ptsU, ptsV, grdPtsU, grdPtsV, 
            r, θ}, 
{hm, 
              grdPrmtrzt} = \({GridPoints, GridParameterization} /. {opts}\) /. 
              Options[PoincareMap]; 
{ptsU, ptsV} = 
            Which[MatchQ[hm, {_Integer?Positive, __Integer?Positive}], hm, 
              MatchQ[hm, _Integer?Positive], {hm, hm}, True, 
              Message[PoincareMap::"BadGridPoints"]; 
              Abort[]]; 
grdPtsU = 
            GridUMin + \((gridUMax - GridUMin)\) 
                Range[0, 
                    ptsU - 1]/\((ptsU - 1)\); 
grdPtsV = 
            GridVMin + \((gridVMax - GridVMin)\) 
                Range[0, 
                    ptsV - 1]/\((ptsV - 
                      1)\); 
Which[
grdPrmtrzt === Polar, 
            grdPrmtrzt = 
              Function[{r, θ}, r {Cos[θ], Sin[θ]}], 
            grdPrmtrzt === Identity, grdPrmtrzt = \(({##} &)\), 
            Not[And[Map[NumericQ, grdPrmtrzt[GridUMin, GridVMin]]]], 
            Message[PoincareMapPlot::"BadParameterization"; 
              Abort[]]]; 
{grdPtsU, grdPtsV, 
            grdPrmtrzt}];\)\)

This utility plots a parametric grid.

gridPlot[grdPtsU_List, grdPtsV_List, grdPrmtrzt_, opts___?OptionQ] := 
        Module[{grdstl, ptsU, ptsV, stlU, stlV, a, b, trajU, 
            trajV}, 
grdstl = GridStyle /. {opts} /. 
              Options[PoincareMapPlot]; 
{ptsU, ptsV} = 
            Map[Length, {grdPtsU, grdPtsV}]; 
stlU = 
            Which[grdstl === Automatic, 
              Map[{Hue[#]} &, Range[0, ptsV]/(2  ptsV)], 
              grdstl == None, Array[{} &, grdPtsV], True, 
              grdstl[[1]]]; 
          stlV = Which[grdstl === Automatic, 
              Map[{Hue[#]} &, -Range[1, ptsU]/(2  ptsU)], 
              grdstl == None, Array[{} &, grdPtsU], True, 
              grdstl[[2]]]; 
trajU = 
            ParametricPlot[
              Evaluate[Map[grdPrmtrzt[a, #] &, grdPtsV]], {a, Min[grdPtsU], 
                Max[grdPtsU]}, DisplayFunction -> Identity, 
              PlotStyle -> stlU, 
              Evaluate[
                FilterOptions[ParametricPlot, 
                  opts]]]; 
trajV = 
            ParametricPlot[
              Evaluate[Map[grdPrmtrzt[#, b] &, grdPtsU]], {b, Min[grdPtsV], 
                Max[grdPtsV]}, DisplayFunction -> Identity, 
              PlotStyle -> stlV, 
              Evaluate[
                FilterOptions[ParametricPlot, 
                  opts]]]; 
Show[{trajU, trajV}, 
            DisplayFunction -> Identity, FilterOptions[Graphics, opts], 
            FilterOptions[Graphics, 
              Sequence @@ Options[PoincareMapPlot]]]];

This function calculates the Poincaré map on a grid of points and interpolates between the values.

PoincareMap[data___] = internalPoincareMap[data];

Clear[internalPoincareMap];

internalPoincareMap[{{x_Symbol, y_Symbol}, {xprime_, 
              yprime_}}, {t_Symbol, t0_?NumericQ, 
            tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
              GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
              GridVMin < gridVMax}, 
          opts___?OptionQ] := internalPoincareMap[{{x, y}, {xprime, 
                yprime}}, {t, t0, 
              tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, opts] = 
          Module[{grdPtsU, grdPtsV, grdPrmtrzt, strtPtsGrid, endPt, x0, y0, 
              u, v, endPtsGrid, 
              pairOfInterp}, 
{grdPtsU, grdPtsV, 
                grdPrmtrzt} = 
              baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
                opts]; 
strtPtsGrid = 
              Apply[grdPrmtrzt, 
                Outer[List, grdPtsU, grdPtsV], {2}]; 
endPt[
                x0_, y0_] := {u[tMax], 
                  v[tMax]} /. NDSolve[{u'[t] == 
                          xprime, v'[t] == yprime, u[0] == x0, 
                        v[0] == y0} /. {x -> u[t], 
                        y -> v[t]}, {u, v}, {t, t0, tMax}, 
                    FilterOptions[NDSolve, 
                      opts]][[1]]; 
endPtsGrid = 
              Apply[endPt, 
                strtPtsGrid, {2}]; 
pairOfInterp = 
              Map[ListInterpolation[#, {{GridUMin, gridUMax}, {GridVMin, 
                        gridVMax}}] &, 
                Transpose[
                  endPtsGrid, {2, 3, 1}]]; 
Function[
              Evaluate[Through[pairOfInterp[#1, #2]]]]];

This function plots the Poincaré map.

PoincareMapPlot[{{x_Symbol, y_Symbol}, {xprime_, 
              yprime_}}, {t_Symbol, t0_?NumericQ, 
            tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
              GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
              GridVMin < gridVMax}, 
          opts___?OptionQ] := 
Module[{onlyEndQ, superImpQ, 
            stkVertQ, samePR, grdPtsU, grdPtsV, grdPrmtrzt, baseGridPlot, 
            poinc, poincPlot, 
            plRng}, 
{onlyEndQ, superImpQ, stkVertQ, 
              samePR} = {OnlyEndGrid, Superimposed, StackVertical, 
                  SamePlotRange} /. {opts} /. 
              Options[PoincareMapPlot]; 
{grdPtsU, grdPtsV, 
              grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
baseGridPlot := 
            gridPlot[grdPtsU, grdPtsV, grdPrmtrzt, 
              DisplayFunction -> Identity, 
              opts]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
poincPlot = 
            gridPlot[grdPtsU, grdPtsV, poinc, 
              DisplayFunction -> Identity, 
              opts]; 
plRng := 
            If[samePR == True, 
              If[$VersionNumber < 4, 
                PlotRange -> 
                  FullOptions[Show[{baseGridPlot, poincPlot}], PlotRange], 
                AbsoluteOptions[Show[{baseGridPlot, poincPlot}], PlotRange]], 
              Sequence @@ {}]; 
Show[
Which[
onlyEndQ === 
                True, poincPlot, 
superImpQ, {baseGridPlot, 
                poincPlot}, 
stkVertQ, 
              GraphicsArray[{{Show[baseGridPlot, plRng]}, {Show[poincPlot, 
                      plRng]}}], 
True, 
              GraphicsArray[{Show[baseGridPlot, plRng], 
                  Show[poincPlot, plRng]}]], 
FilterOptions[
              Graphics, opts], 
            DisplayFunction -> $DisplayFunction]];

Let φ(x, y) be the Poincaré map at some time for some system. Consider the function ψ(x, 
          y) := ∥φ(x, y) - (x, 
              y) (∥ ^2). The following function makes a DensityPlot of Min[ψ(x, y), ϵ] for a give tolerance ϵ. This way the point of the plot that are gray to black are candidates for being a fixed point for φ.

PoincareMapFixedPointDensityPlot[{{x_Symbol, y_Symbol}, {xprime_, 
                yprime_}}, {t_Symbol, t0_?NumericQ, 
              tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
                GridUMin < gridUMax, {GridVMin_?NumericQ, 
                  gridVMax_?NumericQ} /; GridVMin < gridVMax}, 
            Optional[tolerance_, (1/100)], opts___?OptionQ] /; 
          tolerance > 0 := 
Module[{grdPtsU, grdPtsV, 
            grdPrmtrzt, poinc, diff, u, 
            v}, 
{grdPtsU, grdPtsV, grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
rectGrid = 
            Last[baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
                GridParameterization -> Identity, 
                opts]]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
diff[u_, v_] := (# . # &)[
              poinc[u, v] - 
                grdPrmtrzt[u, v]]; 
Show[{DensityPlot[
                Min[diff[u, v], tolerance], {u, GridUMin, gridUMax}, {v, 
                  GridVMin, gridVMax}, DisplayFunction -> Identity, 
                Evaluate[FilterOptions[DensityPlot, opts]], 
                Mesh -> False], 
              gridPlot[grdPtsU, grdPtsV, rectGrid, 
                DisplayFunction -> Identity, PlotPoints -> 2]}, 
            FilterOptions[Graphics, opts], 
            DisplayFunction -> $DisplayFunction, 
            PlotRange -> All]];

The following function makes a 3D plot of the function ψ  instead of a DensityPlot as above.

PoincareMapFixedPointPlot3D[{{x_Symbol, y_Symbol}, {xprime_, 
              yprime_}}, {t_Symbol, t0_?NumericQ, 
            tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
              GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
              GridVMin < gridVMax}, 
          opts___?OptionQ] := 
Module[{grdPtsU, grdPtsV, 
            grdPrmtrzt, poinc, diff, u, 
            v}, 
{grdPtsU, grdPtsV, grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
diff[u_, v_] := (# . # &)[
              poinc[u, v] - grdPrmtrzt[u, v]]; 
          ParametricPlot3D[
            Join[grdPrmtrzt[u, v], {diff[u, v]}], {u, GridUMin, 
              gridUMax}, {v, GridVMin, gridVMax}, 
            Evaluate[FilterOptions[ParametricPlot3D, opts]]]];

The following function tries to find a fixed point of the Poincaré map by minimizing the function ψ defined above near a point (u0,v0) given by the user.

FirstSmallest[x_, 
          f_] := x[[Position[#, Min[#], 1, 1][[1, 1]]]] &[
          f /@ x];

FindPoincareMapFixedPoint[{{x_Symbol, y_Symbol}, {xprime_, 
              yprime_}}, {t_Symbol, t0_?NumericQ, 
            tMax_?NumericQ}, {{GridUMin_?NumericQ, gridUMax_?NumericQ} /; 
              GridUMin < gridUMax, {GridVMin_?NumericQ, gridVMax_?NumericQ} /; 
              GridVMin < gridVMax}, {uFixed_Symbol, vFixed_Symbol}, 
          opts___?OptionQ] := 
Module[{grdPtsU, grdPtsV, 
            grdPrmtrzt, poinc, diff, a, 
            b}, 
{grdPtsU, grdPtsV, grdPrmtrzt} = 
            baseGrid[{GridUMin, gridUMax}, {GridVMin, gridVMax}, 
              opts]; 
poinc = 
            PoincareMap[{{x, y}, {xprime, yprime}}, {t, t0, 
                tMax}, {{GridUMin, gridUMax}, {GridVMin, gridVMax}}, 
              opts]; 
diff[{a_, 
                b_}] := (poinc[a, b][[1]] - grdPrmtrzt[a, 
                        b][[1]])^2 + (poinc[a, 
                        b][[2]] - grdPrmtrzt[a, 
                        b][[2]])^2; 
strtPts = 
            Flatten[Outer[List, grdPtsU, grdPtsV], 
              1]; 
correctedDiff[{a_, b_}] := 
            diff[{a, 
                  b}] + cutoff[{GridUMin, gridUMax}, {GridVMin, 
                    gridVMax}][a, b]; 
fxdPt = 
            FirstSmallest[strtPts, diff]; 
{diff[fxdPt], 
            Thread[{uFixed, vFixed} -> fxdPt]}];

Tail end

End[];

Protect["`*"];

EndPackage[];

On[General::"spell1"]

On[General::"spell"]

Examples

Evaluate the input cells below to see the package in action and learn its usage.

Simple examples of Poincaré map plots

PoincareMapPlot takes a system of two differential equations of the form {({{x' = a(t, x, y)}, {y' = b(t, x, y)}})and draws the initial and final position of a user-specified grid under the Poincaré map at a given time. The system is specified with the syntax {{x, y}, {a, b}}, where x, y are variable names and a, b the right-hand sides:

field = {{x, y}, {1/(1 + x^2), -y + x}};

The grid is given as {{xmin, xmax}, {ymin, ymax}}.

grid = {{0, 1}, {0, 1}};

Here is the graph of the Poincaré map at time 1/2 with all default settings. On the left is the original grid at time 0, and on the right where the grid has drifted at time 1/2

PoincareMapPlot[field, {t, 0, 1/2}, grid];

With the option GridPoints you can control the number of lines in the grid. The default is 5 in each direction:

PoincareMapPlot[field,{t, 0, 1/2},grid,GridPoints -> {3, 5}];

The grid does not need to be cartesian. You can assign a polar grid, for example:

PoincareMapPlot[field,{t, 0, 1/2},{{0, 1}, {0, 2  π}},GridParameterization -> Polar,AspectRatio -> Automatic,GridPoints -> {5, 10}];

You can have an arbitrarily parameterized grid:

PoincareMapPlot[field,{t, 0, 1/2},{{0, 1}, {0, 1}},GridParameterization -> 
              Function[{a, b}, {a + b, a - b}],AspectRatio -> Automatic,GridPoints -> {5, 10}];

You can stack the two plots vertically:

PoincareMapPlot[field,{t, 0, 1/2},{{0, 1}, {0, 1}},GridParameterization -> 
            Function[{a, b}, {a + b, a - b}],AspectRatio -> Automatic,GridPoints -> {5, 10},StackVertical -> True];

You can also decide to draw the final grid only:

PoincareMapPlot[field,{t, 0, 1/2},{{0, 1}, {0, 1}},GridParameterization -> 
            Function[{a, b}, {a + b, a - b}],AspectRatio -> Automatic,GridPoints -> {5, 10},OnlyEndGrid -> True,AxesOrigin -> {0, 0}];

You can also have the initial and final grid superimposed:

PoincareMapPlot[field,{t, 0, 1/2},{{0, 1}, {0, 1}},GridParameterization -> 
            Function[{a, b}, {a + b, a - b}],AspectRatio -> Automatic,GridPoints -> {5, 10},Superimposed -> True];

By default the initial and final grid are drawn with the same PlotRange. You can override this behaviour and let each plot go on its own with the option SamePlotRange->False.

PoincareMapPlot[field,{t, 0, 1},{{0, 1}, {0, 1}},GridParameterization -> 
            Function[{a, b}, {a + b, a - b}],GridPoints -> {5, 10},SamePlotRange -> False];

Poincaré map of x'' (t) = -f(x(t)), and search for periodic solutions of a given period.

Take the following function of one variable

f = x/(1 + x^2) + 3  x (x^2/(1 + x^2));

Here in red is the graph of f. Near x = 0 f(x) ≈ x, while f(x) ≈ 3  x for large x.

Plot[{x, f, 3  x}, {x, 0, 4}, 
        PlotStyle -> {{}, {Thickness[ .01], RGBColor[1, 0, 0]}, {}}, 
        AspectRatio -> 1];

Solutions to x'' (t) = -f(x(t)) that start with small values of (x^2) + (((x'))^2) will be periodic of period approximately  2  π ≈ 6.28, and those with (x^2) + (((x'))^2)large will have period approximately  ((2 π)/(Sqrt[3])) ≈ 3.63. We will try to locate a solution with an intermediate period, namely 5.

The program wants the second-order differential equation x'' (t) = -f(x(t)) to be framed as a system of first-order equation of the general form {({{x' = a(t, x, y)}, {y' = b(t, x, y)}}), and the data to be arranged as {{x, y}, {a(t, x, y), b(t, x, y)}}:

field = {{x, y}, {y, -f}};

Let us take a polar grid with 0 ≤ r ≤ 2, \(TraditionalForm\`\(-π\) ≤ θ ≤ π\):

grid = {{0, 2}, {-π, π}};

Here is the graph of the Poincaré map at time 5:

PoincareMapPlot[field, {t, 0, 5}, grid, GridPoints -> {10, 41}, 
        GridParameterization -> Polar, 
        AspectRatio -> Automatic];

The plot on the left is the starting grid, while on the right you see the position of the grid lines after the map is applied. If for any reason you only want the plot on the right just add the option OnlyEndGrid->True:

PoincareMapPlot[field, {t, 0, 5}, grid, GridPoints -> {10, 41}, 
        GridParameterization -> Polar, AspectRatio -> Automatic, 
        OnlyEndGrid -> True];

You can also have the two grids on top of each other with the option Superimposed->True. You had better do that with less dense grids than the present one. You will se an example later.

If  ϕ(x, y) the Poincaré map at time 5, consider the real function d(x, 
          y) := ∥ϕ(x, y) - (x, 
              y) (∥ ^2). A very small value for d means a solution pretty much periodic with period 5. Here is a 3D plot of  d(x, y):

PoincareMapFixedPointPlot3D[field, {t, 0, 5}, grid, 
        GridPoints -> {10, 41}, PlotPoints -> {20, 41}, 
        GridParameterization -> Polar, 
        AxesLabel -> {"x", "y", "z"}, 
        PlotRange -> {0, 1}];

With a density plot of Min[d(x, y), ϵ] the regions where d(x, y)is smaller than ϵ will stand out as grey or black against a white background. In the following picture ϵ = 1/100. The parameters on the two axes are not x, y but r, θ:

PoincareMapFixedPointDensityPlot[field, {t, 0, 5}, grid, 1/100, 
        GridPoints -> {20, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

You will see two distinct grey-black stripes: one around  r = 0, which corresponds to the trivial periodic solution x(t) ≡ y(t) ≡ 0, whilst the other seems to be the graph of a function from θ to r. Let us zoom around the gray point r =  .75, θ = 0:

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 5}, {{ .65,  .85}, {- .2,  .2}}, 1/1000, 
        GridPoints -> {20, 21}, PlotPoints -> {20, 20}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

Let us close in on the fixed point by minimizing the function  d near  r =  .756, θ = 0 with the built-in minimization routines:

{squareDistance, initialData} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 5}, {{ .65,  .85}, {- .1,  .1}}, {r0, θ0}, 
        GridPoints -> {40, 11}, 
        GridParameterization -> Polar]

Here is a plot of the solution corresponding the fixed point found above, seen as a function t ↦ x(t):

solutionX = 
        NDSolve[{x''[t] == (-f /. x -> x[t]), 
            x[0] == r0*Cos[θ0] /. 
              initialData, x'[0] == 0}, 
          x, {t, 0, 5}];
Plot[x[t] /. solutionX, {t, 0, 5}];

A check to see how much the values at t = 0 and t = 5 are close to each other:

{x[5] - x[0], x'[5] - x'[0]} /. solutionX

A graph of our periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

solutionXY = 
        NDSolve[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == r0*Cos[θ0] /. initialData, 
            y[0] == 0}, {x, y}, {t, 0, 
            5}];
phaseSpacePlot = 
        ParametricPlot[{x[t], y[t]} /. solutionXY, {t, 0, 5}, 
          AspectRatio -> Automatic, 
          PlotStyle -> {Thickness[ .007], 
              Dashing[{ .02,  .04}]}];

The following picture superimposes the 5-periodic solution to the Poincare map plot:

Show[{PoincareMapPlot[field, {t, 0, 5}, grid, 
            GridPoints -> {10, 41}, GridParameterization -> Polar, 
            OnlyEndGrid -> True, AspectRatio -> Automatic, 
            DisplayFunction -> Identity], phaseSpacePlot}, 
        DisplayFunction -> $DisplayFunction];

A zoom of Poincaré map around the starting point of the 5-periodic solution, with the initial and the final grids superimposed:

Show[{phaseSpacePlot, 
          PoincareMapPlot[
            field, {t, 0, 
              5}, {0.756650394233193` + {- .1,  .1}, {- .05,  .05}}, 
            GridPoints -> {7, 5}, GridParameterization -> Polar, 
            Superimposed -> True, AspectRatio -> Automatic, 
            DisplayFunction -> Identity]}, 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> {{ .62,  .86}, {- .282,  .22}}, 
        Frame -> True, Axes -> False];

Animated Poincaré maps

Plots of Poincaré maps at a sequence of times can be shown as a movie. When all the frames have been produced, double-click the first cell in the sequence to start the animation

f = x/(1 + x^2) + 3  x (x^2/(1 + x^2));

field = {{x, y}, {y, -f}};

grid = {{1/4, 2}, {-1, 1}};

filmFrame[instant_] := 
        PoincareMapPlot[field, {t, 0, instant}, grid, 
          OnlyEndGrid -> True, GridPoints -> {5, 5}, 
          GridParameterization -> Polar, AspectRatio -> Automatic, 
          PlotRange -> {{-3, 3}, {-3, 3}}, 
          AxesOrigin -> {0, 0}];

Do[filmFrame[τ], {τ, 0, 5, 1/5}]

Poincaré map of x'' (t) = -f(t, x(t)), and search for periodic solutions of period 2  π.

Take a function f(t, x) that is 2  π-periodic with respect to t and that is asymptotic to x/2 as x -> 0 and to 3  x as x -> +∞:

f = 
        x/(2 + x^2) + 3  x (x^2/(1 + x^2)) + 
          x^2*Cos[t]/(3  x^4 + 2);

Solutions to x'' (t) = -f(t, x(t)) that start with small values of (x^2) + (((x'))^2) will be periodic of period approximately  (Sqrt[8]) π ≈ 8.9, and those with (x^2) + (((x'))^2)large will have period approximately  ((2 π)/(Sqrt[3])) ≈ 3.63. We will try to locate solutions which are periodic with the same t-period as f(t, x), that is, 2  π-periodic solutions. Any such solution will correspond to a fixed point for the Poincaré map at time 2  π.

field = {{x, y}, {y, -f}};

grid = {{0, 2}, {-π, π}};

Here is a plot of Poincaré map:

PoincareMapPlot[field, {t, 0, 2  π}, grid, 
        GridPoints -> {20, 41}, GridParameterization -> Polar, 
        AspectRatio -> Automatic];

If  ϕ(x, y) the Poincaré map at time  2  π, consider the real function d(x, 
          y) := ∥ϕ(x, y) - (x, 
              y) (∥ ^2). A very small value for d means a solution pretty much periodic with period 2  π. With a density plot of Min[d(x, y), ϵ] the regions where d(x, y)is smaller than ϵ will stand out as grey or black against a white background. In the following picture ϵ = 1/1000. The parameters on the two axes are not x, y but r, θ:

PoincareMapFixedPointDensityPlot[field, {t, 0, 2  π}, grid, 
        1/1000, GridPoints -> {40, 81}, PlotPoints -> {100, 100}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

We make out one grey stripe plus two grey spots. The stripe around r = 0 corresponds to the trivial null solution. The two grey spots are approximately at the points ( .45, 0) and ( .7, ±π).

A zoom around (.45,0):

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 2  π}, {{ .42,  .48}, {- .5,  .5}}, 1/10000, 
        GridPoints -> {40, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

Minimizing d(x, y) with numerical methods we close in on the fixed point:

{squareDistance1, initialData1} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 
          2  π}, {{ .42,  .48}, {- .5,  .5}}, {r0, θ0}, 
        GridPoints -> {40, 11}, 
        GridParameterization -> Polar]

Here is a plot of the solution corresponding the fixed point found above, seen as a function t ↦ x(t):

solutionX1 = 
        NDSolve[{x''[t] == (-f /. x -> x[t]), 
            x[0] == r0*Cos[θ0] /. 
              initialData1, x'[0] == 0}, 
          x, {t, 0, 2  π}];
Plot[x[t] /. solutionX1, {t, 0, 2  π}, 
        Ticks -> {PiScale[0, 2  π], Automatic}];

A check to see how much the values at t = 0 and t = 2  π are close to each other:

{x[2  π] - x[0], x'[2  π] - x'[0]} /. 
      solutionX1

A graph of our periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

solutionXY1 = 
        NDSolve[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == r0*Cos[θ0] /. initialData1, 
            y[0] == 0}, {x, y}, {t, 0, 
            2  π}];
ParametricPlot[{x[t], y[t]} /. solutionXY1, {t, 0, 2  π}, 
        AspectRatio -> Automatic];

A zoom of Poincaré map around the starting point of this 2  π-periodic solution, with the initial and the final grids superimposed:

PoincareMapPlot[
        field, {t, 0, 
          2  π}, {{0.44404 -  .1, 0.44404 +  .1}, {- .1,  .1}}, 
        GridPoints -> {5, 5}, GridParameterization -> Polar, 
        AspectRatio -> Automatic, Superimposed -> True, 
        Frame -> True, Axes -> False];

A zoom of the density plot around the second grey spot, the one near ( .7, ±π):

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 2  π}, {{ .6,  .8}, {π -  .1, π +  .1}}, 
        1/10000, GridPoints -> {40, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> Polar, 
        FrameLabel -> {"r", "θ"}];

Minimize d(x, y) in that spot:

{squareDistance2, initialData2} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 2  π}, {{ .6,  .8}, {3, 3.2}}, {r0, θ0}, 
        GridPoints -> {40, 11}, 
        GridParameterization -> Polar]

Here is a plot of the solution corresponding the second fixed point, seen as a function t ↦ x(t):

solutionX2 = 
        NDSolve[{x''[t] == (-f /. x -> x[t]), 
            x[0] == r0*Cos[θ0] /. 
              initialData2, x'[0] == 0}, 
          x, {t, 0, 2  π}];
Plot[x[t] /. solutionX2, {t, 0, 2  π}];

A check to see how much the values at t = 0 and t = 2  π are close to each other:

{x[5] - x[0], x'[5] - x'[0]} /. solutionX2

A graph of our second periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

solutionXY2 = 
        NDSolve[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == r0*Cos[θ0] /. initialData2, 
            y[0] == 0}, {x, y}, {t, 0, 
            2  π}];
ParametricPlot[{x[t], y[t]} /. solutionXY2, {t, 0, 2  π}, 
        AspectRatio -> Automatic];

A zoom of Poincaré map around the starting point of the second 2  π-periodic solution, with the initial and the final grids superimposed:

PoincareMapPlot[
        field, {t, 0, 
          2  π}, {{0.70665 -  .1, 
            0.70665 +  .1}, {π -  .1, π +  .1}}, 
        GridPoints -> {5, 5}, GridParameterization -> Polar, 
        AspectRatio -> Automatic, Superimposed -> True, 
        Frame -> True, Axes -> False];

Poincaré map of x'' (t) = -f(x(t)) for which there is an invariant annulus

Consider this function of x:

f = 2  x^3;

The solutions of x'' = -f(x) starting with small (x^2) + (((x'))^2) will have large period, while the ones starting with large (x^2) + (((x'))^2) will have small period. We set out to find periodic solutions with period 4.

The orbits of this system, seen as {({{x' = y}, {y' = -f(x)}}) lie on the level sets of the function (x, y) ↦ (x^4) +(y^2), which look somewhat like ellipses. What follows is a pseudo-polar parameterization of the plane where  θ is the usual argument but r = (Sqrt[(x^4) + (y^2)]). The advantage of this parameterization is that the sets r = constant coincide with the level sets of (x, y) ↦ (x^4) +(y^2). The main parameterization has a singularity along the y-axis, so that we need a separate formula for that case.

myParameterization[
          r_?NumericQ, θ_?NumericQ /; 
            Abs[Cos[θ]] < $MachineEpsilon] = {0, 
          r*Sign[Sin[θ]]};
myParameterization[
          r_?NumericQ, θ_?NumericQ /; 
            Abs[Cos[θ]] ≥ $MachineEpsilon] = 
        Sign[Cos[θ]] (Sqrt[-(1/2) Tan[θ]^2 + 
                (1/2) (Sqrt[4 (r^2) + Tan[θ]^4])]) {1, 
            Tan[θ]};

Since the system is autonomous, if there is a solution with period 4, all points (x(t), y(t)) along its trajectory are starting points of a solution with period 4. Hence the fixed points of the Poincaré map are not isolated points, but a whole line. Actually, it will be a whole level set of the function (x, y) ↦ (x^4) +(y^2).

The data structure for the differential system:

field = {{x, y}, {y, -f}};

Let us take a pseudo-polar grid with 1 ≤ r ≤ 2, \(TraditionalForm\`\(-π\) ≤ θ ≤ π\):

grid = {{1, 2}, {-π, π}};

A plot of the Poincaré map at time 4:

PoincareMapPlot[field, {t, 0, 4}, grid, GridPoints -> {7, 25}, 
        GridParameterization -> myParameterization, 
        AspectRatio -> Automatic];

If  ϕ(x, y) the Poincaré map at time 4, consider the real function d(x, 
          y) := ∥ϕ(x, y) - (x, 
              y) (∥ ^2). A very small value for d means a solution pretty much periodic with period 4. Here is a 3D plot of  d(x, y):

PoincareMapFixedPointPlot3D[field, {t, 0, 4}, grid, 
        GridPoints -> {10, 41}, PlotPoints -> {20, 41}, 
        GridParameterization -> myParameterization, 
        AxesLabel -> {"x", "y", "z"}, 
        PlotRange -> {0, 1}];

With a density plot of Min[d(x, y), ϵ] the regions where d(x, y)is smaller than ϵ will stand out as grey or black against a white background. In the following picture ϵ = 1/100. The parameters on the two axes are not x, y but r, θ:

PoincareMapFixedPointDensityPlot[field, {t, 0, 4}, grid, 1/100, 
        GridPoints -> {20, 41}, PlotPoints -> {50, 50}, 
        GridParameterization -> myParameterization, 
        FrameLabel -> {"r", "θ"}];

We make out a grey-black stripe that looks like the graph of a set  r = constant. Let us zoom in on the point r = 1.72, θ = 0:

PoincareMapFixedPointDensityPlot[
        field, {t, 0, 4}, {{1.6, 1.8}, {- .2,  .2}}, 1/1000, 
        GridPoints -> {20, 21}, PlotPoints -> {20, 20}, 
        GridParameterization -> myParameterization, 
        FrameLabel -> {"r", "θ"}];

Let us minimize d on a grid of points in the rectangle 1.6 ≤ r ≤ 
        1.8,  - .2 ≤ θ ≤  .2:

{squareDistance, initialData} = 
      FindPoincareMapFixedPoint[
        field, {t, 0, 4}, {{1.6, 1.8}, {- .2,  .2}}, {r0, θ0}, 
        GridPoints -> {40, 21}, 
        GridParameterization -> myParameterization]

Here is a plot of the solution corresponding the fixed point found above, seen as a function t ↦ x(t):

solutionX = 
        NDSolve[{x''[t] == (-f /. x -> x[t]), 
            x[0] == (myParameterization[r0, θ0] /. 
                    initialData)[[1]], x'[
                0] == (myParameterization[r0, θ0] /. 
                    initialData)[[2]]}, 
          x, {t, 0, 4}];
Plot[x[t] /. solutionX, {t, 0, 4}];

A check to see how much the values at t = 0 and t = 4 are close to each other:

{x[4] - x[0], x'[4] - x'[0]} /. solutionX

A graph of our 4-periodic solution in the phase space, i.e., seen as a function t ↦ (x(t), y(t)):

solutionXY = 
        NDSolve[{x'[t] == 
              y[t], y'[t] == (-f /. x -> x[t]), 
            x[0] == (myParameterization[r0, θ0] /. 
                    initialData)[[1]], 
            y[0] == (myParameterization[r0, θ0] /. 
                    initialData)[[2]]}, {x, y}, {t, 0, 
            4}];
ParametricPlot[{x[t], y[t]} /. solutionXY, {t, 0, 4}, 
        AspectRatio -> Automatic];