Math 331.1  •  19 October 2005

van der Pol equation

Copyright © 2005 Murray Eisenberg.  All rights reserved.

Try to solve exactly

DSolve[
{
x''[
              t]  -  (1  -  x[t]^2) x'[t]  +  
            x[t]  ==  0,  
x[0]  ==  
          1,  x'[0]  ==  
          1
},  
x[t],  t]

Put into system form: the vector field

Clear[F]
F[t_,  {x_,  
          y_}]  :=  {y,  -x  +  (1  -  (x^2)) y}

Euler's Method for Systems

The basic step:

eulerStep[F_,  t_,  Y_,  Δt_]  :=  
      Y  +  Δt F[t,  Y]

Iterating the basic step and collecting the results:

euler[F_,  t0_,  Y0_,  Δt_,  
        b_]  := 
Module[{t,  Y,  n,  
          k}, 
n  =  
          Round[((b  -  t0)/Δ)t]; 
t[
            0]  =  t0;  
        Y[0]  =  
          Y0; 
Do [
t[k]  =  
            t0 +  k Δt; 
Y[k]  =  
            eulerStep[F,  t[k  -  1],  
              Y[k  -  1],  Δt], 
{k,  
            1,  n}
]; 
Table[{t[k],  
            Y[k]},  {k,  0,  n}]
]

Pretty-printing the result of function euler:

padLast[v_]  :=  {v[[1]],  
        Map[NumberForm[#,  {8,  6},  
              NumberPadding  ->  {" ",  ""}] &,  
          v[[2]]]}

prettyEuler[F_,  t0_,  Y0_,   Δt_,  b_]  :=  
      TableForm[
        padLast  /@  euler[F,  t0,  Y0,  Δt,  b],  
        TableHeadings  ->  {None,  (TraditionalForm  /@  {HoldForm[(Subscript[t, i])],  HoldForm[(Subscript[Y, i])]})}]

Fancy plotting in phase plane:

eulerTrajectory[F_,  t0_,  Y0_,   Δt_,  
        b_]  :=  
Module[{pts}, 
pts  =  Last  /@  euler[F,  t0,  Y0,  Δt,  b]; 
        Show[{ListPlot[pts,  PlotStyle  ->  PointSize[0.015],  
              PlotJoined  ->  False,  AspectRatio  ->  1,  
              DisplayFunction  ->  
                Identity],  
ListPlot[pts,  
              PlotStyle  ->  {Red,  Thickness[0.01]},  
              PlotJoined  ->  True,  AspectRatio  ->  1,  
              DisplayFunction  ->  
                Identity], 
ListPlot[pts,  
              PlotStyle  ->  {Blue,  PointSize[0.015]},  
              PlotJoined  ->  False,  AspectRatio  ->  1,  
              DisplayFunction  ->  Identity]
},  
          DisplayFunction  ->  $DisplayFunction];
]

Use Euler's Method for van der Pol equation

Step size Δt  =  0.2 for 0  ≤  t  ≤  8, initial point ((Subscript[x, 0]),  (Subscript[y, 0]))  =  (1,  1)

Iterating manually:

eulerStep[F,  0,  {1,  1},  0.2]

eulerStep[F,  0.2,  %,  0.2]

eulerStep[F,  0.4,  %,   0.2]

Iterating automatically:

euler[F,  0,  {1,  1},  0.2,  8]

prettyEuler[F,  0,  {1,  1},  0.2,  8]

Plotting the trajectory:

pts  =  Last  /@  euler[F,  0,  {1,  1},  0.2,  8]

ListPlot[pts,  PlotJoined  ->  True,  
        AspectRatio  ->  1];

eulerTrajectory[F,  0,  {1,  1},  0.2,  8]

Step size Δt  =  0.1 for 0  ≤  t  ≤  8, initial point ((Subscript[x, 0]),  (Subscript[y, 0]))  =  (1,  1)

eulerTrajectory[F,  0,  {1,  1},  0.1,  8]

Step size Δt  =  0.1 for 0  ≤  t  ≤  16, initial point ((Subscript[x, 0]),  (Subscript[y, 0]))  =  (1,  1)

eulerTrajectory[F,  0,  {1,  1},  0.1,  16]

Step size Δt  =  0.1 for 0  ≤  t  ≤  32, new ((Subscript[x, 0]),  (Subscript[y, 0]))  =  (0.5,  0.25)

eulerTrajectory[F,  0,  {0.5,  0.25},  0.1,  32]

Step size Δt  =  0.1 for 0  ≤  t  ≤  32, new ((Subscript[x, 0]),  (Subscript[y, 0]))  =  (0.5,  4)

eulerTrajectory[F,  0,  {0.5,  4},  0.1,  32]

A conclusion

There seems to be a stable limit cycle for the van der Pol oscillator.

M. L. Cartwright discovered the period of that limit cycle.