Regular Power Series Solution

Solution of
y''+cos(x) y=0
for y(x) with y(0)=1 and y'(0)=0.
An exact solution can be found in this case, and can
be compared with the series solution.

In[1]:=

Off[General :: spell1]

Let's find a solution valid out to x^10.  The first step in finding a series solution is to find an expansion for Cos[x]:

In[2]:=

n = 10 ;

In[3]:=

cosser = Cos[x] + O[x]^(n + 1)

Out[3]=

1 - x^2/2 + x^4/24 - x^6/720 + x^8/40320 - x^10/3628800 + O[x]^11

In[4]:=

y[x_] = Sum[a[k] x^k, {k, 0, n, 2}] + O[x]^(n + 1)

Out[4]=

a[0] + a[2] x^2 + a[4] x^4 + a[6] x^6 + a[8] x^8 + a[10] x^10 + O[x]^11

In[5]:=

deqn = y ''[x] + cosser * y[x] == 0

Out[5]=

(a[0] + 2 a[2]) + (-a[0]/2 + a[2] + 12 a[4]) x^2 + (a[0]/24 - a[2]/2 + a[4] + 30 a[6]) x^4 + (-a[0]/720 + a[2]/24 - a[4]/2 + a[6] + 56 a[8]) x^6 + (a[0]/40320 - a[2]/720 + a[4]/24 - a[6]/2 + a[8] + 90 a[10]) x^8 + O[x]^9 == 0

Here are the unknowns in the expansion of y(x) that we must solve for:

In[6]:=

coeffs = Table[a[k], {k, 2, n, 2}]

Out[6]=

{a[2], a[4], a[6], a[8], a[10]}

Here are the algebraic equations for those coefficients:

In[7]:=

eqns = LogicalExpand[deqn]

Out[7]=

a[0] + 2 a[2] == 0 && -a[0]/2 + a[2] + 12 a[4] == 0 && a[0]/24 - a[2]/2 + a[4] + 30 a[6] == 0 && -a[0]/720 + a[2]/24 - a[4]/2 + a[6] + 56 a[8] == 0 && a[0]/40320 - a[2]/720 + a[4]/24 - a[6]/2 + a[8] + 90 a[10] == 0

We actually don't need the TableForm for the math, it is just nice to look at:

In[8]:=

TableForm[ReplaceAll[eqns, And -> List]]

Out[8]//TableForm=
a[0] + 2 a[2] == 0
-a[0]/2 + a[2] + 12 a[4] == 0
a[0]/24 - a[2]/2 + a[4] + 30 a[6] == 0
-a[0]/720 + a[2]/24 - a[4]/2 + a[6] + 56 a[8] == 0
a[0]/40320 - a[2]/720 + a[4]/24 - a[6]/2 + a[8] + 90 a[10] == 0

In[9]:=

solset = Solve[eqns, coeffs]

Out[9]=

{{a[2] -> -a[0]/2, a[4] -> a[0]/12, a[6] -> -a[0]/80, a[8] -> (11 a[0])/8064, a[10] -> -(17 a[0])/129600}}

In[10]:=

ys[x_] = y[x] /. solset[[1]]

Out[10]=

a[0] - 1/2 a[0] x^2 + 1/12 a[0] x^4 - 1/80 a[0] x^6 + (11 a[0] x^8)/8064 - (17 a[0] x^10)/129600 + O[x]^11

In[11]:=

yss[x_] = Normal[ys[x]]

Out[11]=

a[0] - 1/2 x^2 a[0] + 1/12 x^4 a[0] - 1/80 x^6 a[0] + (11 x^8 a[0])/8064 - (17 x^10 a[0])/129600

So here is the series solution that satisfies the specified initial conditions:

In[12]:=

yss1[x_] = yss[x] /. a[0] -> 1

Out[12]=

1 - x^2/2 + x^4/12 - x^6/80 + (11 x^8)/8064 - (17 x^10)/129600

Now find exact solution:

In[13]:=

Clear[y]

In[14]:=

exact = DSolve[{y ''[x] + Cos[x] y[x] == 0, y '[0] == 0, y[0] == 1}, y[x], x]

Out[14]=

{{y[x] -> MathieuC[0, -2, x/2]/MathieuC[0, -2, 0]}}

In[15]:=

ey[x_] = y[x] /. exact[[1]]

Out[15]=

MathieuC[0, -2, x/2]/MathieuC[0, -2, 0]

Plot approximate (red) and exact (green) solutions:

In[16]:=

Plot[{ey[x], yss1[x]}, {x, 0, 10}, PlotRange -> {-10, 10}, PlotStyle -> {RGBColor[0, 1, 0], RGBColor[1, 0, 0]}]

[Graphics:HTMLFiles/index_35.gif]

Out[16]=

-Graphics -

Converted by Mathematica  (March 3, 2003)