Phase portraits, Poincare maps and Stable Manifolds for the Van der Pol equation.

First write a function to take an initial condition (q0, v0) and produce an orbit.  The orbit will be a vector-valued function (q(t), v(t))

VdPOrbit[{q0_,v0_},{a_,b_},eps_]:={q[t],v[t]}/.NDSolve[{q'[t]==v[t],v'[t]==-q[t]-eps*(q[t]^2-1)v[t],q[0]==q0,v[0]==v0},{q[t],v[t]},{t,a,b}][[1]]

Here is a typical orbit with eps = 0.5 which seems to converge to a periodic orbit.

orb1=VdPOrbit[{1,0},{0,20},0.5];
ParametricPlot[orb1,{t,0,20}]

We can use the Manipulate command to produce an interactive orbit calculator.

MakeVdPPlot[{q0_,v0_},{a_,b_},eps_]:= With[{orb=VdPOrbit[{q0,v0},{a,b},eps]},
ParametricPlot[orb,{t,a,b},PlotRange->{{-4,4},{-4,4}}]]

Manipulate[MakeVdPPlot[p,{0.001,tmax},eps],{{p,{0,0}},Locator},{{tmax,5},0,20},{{eps,0.5},0,2}]

Evidently there is an attracting periodic orbit.  Now we fix eps = 0.5 and try to compute the stable bundle and stable foliation of the attracting orbit.  First we find the attractor numerically and compute the period.  To find the attractor we first construct the Poincaré  return map where the section is the positive q axis.  This is conveniently done using an event detector in the NDSolve command.  The event will be that the v coordinate reaches 0 and changes from positive to negative.  We also need a return time function.  Then the Poincare map of the q axis is the first component of the return point.

VdPOrbitUntilReturn[{q0_,v0_},eps_]:={q[t],v[t]}/.NDSolve[{q'[t]==v[t],v'[t]==-q[t]-eps*(q[t]^2-1)v[t],q[0]==q0,v[0]==v0},{q[t],v[t]},{t,0,100},Method->{ "EventLocator","Direction"->1,"Event"->(v[t]<0)}][[1]]

ReturnTime[{q0_,v0_},eps_]:= With[{orb=VdPOrbitUntilReturn[{q0,v0},eps]},
orb[[1]][[0]]["Domain"[]][[1,2]]]

PoincareMap[q0_?NumericQ,eps_]:=Module[{q1},
NDSolve[{q'[t]==v[t],v'[t]==-q[t]-eps*(q[t]^2-1)v[t],q[0]==q0,v[0]==0},{q,v},{t,0,100},Method->{ "EventLocator","Direction"->1,"Event"->(v[t]<0),"EventAction":>{q1=q[t];"StopIntegration" }}];
q1]

orb=VdPOrbitUntilReturn[{3,0},0.5];
ParametricPlot[orb,{t,0,ReturnTime[{3,0},0.5]}]

Here is a plot of the Poincare map as the initial value of q runs from 1 to 4

Plot[{q0,PoincareMap[q0,0.5]},{q0,0,4},PlotStyle->Thick]

We can solve numerically for the fixed point {q0,0}.  Using this we can get the period of the periodic orbit.

fp = {q0,0}/.FindRoot[PoincareMap[q0,0.5]==q0,{q0,1}]
perorb=VdPOrbitUntilReturn[fp,0.5];
period = ReturnTime[fp,0.5]
perorbplot= ParametricPlot[perorb,{t,0,period},PlotStyle->{Thick,Black}]

To find the stable bundle we need to solve the variational equations.  Here is a function which solves the variational equation starting from a given initial point with a given initial vector.

VarOrbitAlongVdPOrbit[{q0_,v0_},{dq0_,dv0_},{a_,b_},eps_]:={q[t],v[t],dq[t],dv[t]}/.NDSolve[{q'[t]==v[t],v'[t]==-q[t]-eps*(q[t]^2-1)v[t],
dq'[t]==dv[t],dv'[t]==-dq[t]-eps*(q[t]^2-1)dv[t]-2*eps*q[t]v[t]dq[t],
q[0]==q0,v[0]==v0,dq[0]==dq0,dv[0]==dv0},{q[t],v[t],dq[t],dv[t]},{t,a,b}][[1]]

To find the columns of the monodromy matrix, solve the variational equations for one time period starting with initial vectors {1,0} and {0,1}.

VarOrbitAlongVdPOrbit[fp,{1,0},{0,period},0.5]/.t->period;
{a,c} =  %[[{3,4}]];
VarOrbitAlongVdPOrbit[fp,{0,1},{0,period},0.5]/.t->period;
{b,d} =  %[[{3,4}]];
PhiT = {{a,b},{c,d}};
PhiT//MatrixForm

Next we find the eigenvalues and single out the stable eigenvector.

Eigenvalues[PhiT]
Eigenvectors[PhiT]
es0 = %[[2]]

To get the stable bundle, just solve the variational equation using es0 as initial vector.  Then we can plot a line segment along E^s at each point of the orbit to visualize the stable bundle.  Or plot a parametric surface based on these.

es= VarOrbitAlongVdPOrbit[fp,es0,{0,2*period},0.5];

ParametricPlot[es[[{1,2}]]+s*es[[{3,4}]]/EuclideanDistance[es[[{3,4}]],{0,0}],{t,0, period},{s,-0.6,0.6},MeshFunctions->{#3&},Mesh->40,MeshStyle->Thick,PlotRange->{{-3.5,3.5},{-3.5,3.5}},PlotPoints->30];
Show[%,perorbplot]

The stable bundle gives an approximation to the stable foliation (the stable manifold foliated into the stable manifolds of the individual points of the orbit).  We can get a better approximation by integrating the backward vectorfield starting at points of the stable subspace very close to the orbit.  To approximate the stable curve at x, first flow x forward for time deltat.  Then choose n points on the stable subspace of phi(deltat, x) very close to this new point.   This approximates a segment on the stable manifold of the new point.  Following this backwards by time deltat, the segment stretches out into an appoximation of the stable manifold of the original point.

BackwardVdPOrbit[{q0_,v0_},{a_,b_},eps_]:={q[t],v[t]}/.NDSolve[{q'[t]==-v[t],v'[t]==q[t]+eps*(q[t]^2-1)v[t],q[0]==q0,v[0]==v0},{q[t],v[t]},{t,a,b}][[1]]

BackwardTimeTMap[{q0_,v0_},elapsedtime_,eps_]:=({q[t],v[t]}/.NDSolve[{q'[t]==-v[t],v'[t]==q[t]+eps*(q[t]^2-1)v[t],q[0]==q0,v[0]==v0},{q[t],v[t]},{t,0,elapsedtime}][[1]])/.t->elapsedtime

GetStableSubspacePoints[t0_,n_]:= With[{est0= es/.t->t0},
Table[est0[[{1,2}]]+0.1*h^3*est0[[{3,4}]]/EuclideanDistance[est0[[{3,4}]],{0,0}],{h,-0.7,1.3,2/n}]]

FlowCurvePointsBackward[points_,elapsedtime_,eps_]:=Map[BackwardTimeTMap[#,elapsedtime,eps]&,points]

ThrowAwayPoints[list_,maxsize_]:=Select[list,Function[q,q[[1]]^2+q[[2]]^2≤maxsize^2]]

StableManifold[t0_,n_,deltat_]:=FlowCurvePointsBackward[GetStableSubspacePoints[t0+deltat,n],deltat,0.5]

Here is a plot of the stable manifolds of 20 points on the orbit.  There are some difficulties with points going to infinity under backward integration, so it took some tweaking to get the picture !

ListPlot[Table[ThrowAwayPoints[StableManifold[t0,300,6],6],{t0,0,period-0.001,period/20}],Joined->True,PlotRange->{{-4,4},{-4,4}},PlotStyle->Blue];
Show[perorbplot,%,PlotRange->{{-4,4},{-4,4}}]

The stable manifold of the periodic orbit is the whole plane except (0,0).  This is foliated into the curves of the stable foliation. These curves are sometimes called isochrons.  All of the points on one of the curves converge to one point on the moving orbit, so they are somehow synchronized in time with one another.