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.

In[1]:=

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:

In[2]:=

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

In[3]:=

tryit = nonosc[1.5]

Out[3]=

{{y -> InterpolatingFunction[{{0.`, 10.`}}, <>]}}

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

In[4]:=

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

[Graphics:HTMLFiles/index_7.gif]

Out[4]=

-Graphics -

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]]:

In[5]:=

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

Let's give it a try:

In[6]:=

yfun[2]

Out[6]=

-0.9652002501620248`

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

In[7]:=

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

Out[7]=

{t -> 1.0839591623123357`}

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

In[8]:=

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:

In[9]:=

FirstZero[-.1]

Out[9]=

1.6334535024077244`

In[10]:=

FirstZero[.1]

Out[10]=

1.5151623510238674`

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

In[11]:=

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

Out[11]=

{1.274611126722762`, 1.1998531515165958`, 1.1379310613803093`, 1.0853548439050222`, 1.0398887968758344`, 0.9999988917456456`, 0.9645823110087438`, 0.9328352499027122`, 0.9041442188648122`, 0.8780344030199516`, 0.8541306198179611`}

In[12]:=

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

In[13]:=

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

Out[13]=

{1.2307692307692308`, 1.1764705882352942`, 1.1267605633802817`, 1.081081081081081`, 1.0389610389610389`, 1.`, 0.9638554216867469`, 0.9302325581395349`, 0.898876404494382`, 0.8695652173913044`, 0.8421052631578947`}

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

In[14]:=

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

[Graphics:HTMLFiles/index_25.gif]

Out[14]=

-Graphics -

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

In[15]:=

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

In[16]:=

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

[Graphics:HTMLFiles/index_30.gif]

Out[16]=

-Graphics -

Converted by Mathematica  (March 3, 2003)