A Nonlinear Oscillator,
Numerical Solution

Problem 3, on P.S. #3, 2002.  Numerical check of perturbation 
solution, problem not solved here.

Finds a numerical solution to:
y'' + y +((ϵy)^3)=0 
with y'(0)=0 and y(0)=1.
The period in the numerical solution is then compared to the prediction made from the perturbation expansion.  The math for the
perturbation expansion is not shown here.

Clear[eps]

The following line was changed to text.  You may want to change the Format>Style to "Input" to execute the statement, and confirm that Mathematica cannot find an analytical solution:

DSolve[{y''[t]+y[t]+eps*y[t]^3==0,y[0]==1,y'[0]==0},y[t],t]

So we resort to a numerical solution for a specified value of eps:

nonosc[eps_]:=NDSolve[{y''[t]+y[t]+eps*y[t]^3==0,y[0]==1,y'[0]==0},{y},{t,0,10}]

tryit=nonosc[1.5]

Plot the analytical solution  for eps=0 in black and the numerical
solution for eps=0.5  in red:

Plot[{Cos[t],y[t]/.tryit},{t,0,10},PlotStyle->{RGBColor[0,0,0],RGBColor[1,0,0]}]

Here we make a more convenient interface to access the interpolating function that contains the numerical solution.  Notice the replacement rule is inside a nested list, which is accessed by the [[1,1]]:

yfun[t_]=y[t]/.tryit[[1,1]];

Let's give it a try:

yfun[2]

Let's find  t where y(t)=0 in the numerical solution:

FindRoot[yfun[t],{t,{0,1}}]

Now make a module to find the first t where y(t)=0 for a
specified value of ϵ.

FirstZero[eps_]:=Module[{t,tryit,yfun,rt},tryit=nonosc[eps];
yfun[t_]=y[t]/.tryit[[1,1]];rt=FindRoot[yfun[t]==0,{t,1,2}];
t/.rt[[1]]
]

Try it:

FirstZero[-.1]

FirstZero[.1]

Make a table of ratio of the period of oscillation to that with ϵ=0.

Table[FirstZero[x]*2/Pi,{x,-.5,.5,.1}]

PertResult1[eps_]:=1/(1+3*eps/8)

Table[PertResult1[x],{x,-.5,.5,.1}]

Compare the numerical result, which we assume to be exact, (green) with the perturbation result (red).

Plot[{FirstZero[x]*2/Pi,PertResult1[x]},{x,-.5,.5},PlotStyle->{RGBColor[0,1,0],RGBColor[1,0,0]}]

Now compare with the perturbation result out to (ϵ^2):

PertResult2[eps_]:=1/(1+3*eps/8 -21*eps^2/256)

Plot[{FirstZero[x]*2/Pi,PertResult2[x]},{x,-.5,.5},PlotStyle->{RGBColor[0,1,0],RGBColor[1,0,0]}]