Demo 5.2

Series Solutions near an ordinary point.

We consider the equation  P(x)y'' + Q(x)y' +R(x)y = 0.   We seek a series solution near (Subscript[x, 0]).  We will begin by assuming that  (Subscript[x, 0])  =  0.   We first define a generic power series with terms up to order n.  
The last term stands for terms with powers of x higher than n.

s[n_] := (∑_(i = 0)^n) (Subscript[a, i]) (x^i) + O[x]^(n + 1)

s[5]

Notice that taking the derivative causes the order (that is the highest exponent of the polynomial  given explicitely) to drop.

Define for abbreviation the derivatives of first and second order as follows.

\(D1[F_, x_] := D[F, {x, 1}]; 
D2[F_, x_] := D[F, {x, 2}]\)

D1[s[5], x]

And this is the second derivative.

D2[ s[5], x]

The following is the first few terms of  s + s''.

s[5] + D2[ s[5], x]

Example 1.               y'' + xy' + y = 0

The coefficients of the equation are

{P[x_], Q[x_], R[x_]} = {1, x, 1}

We want to substitute s[n] into the differential equation.  We use the sign := because we can't evaluate until we give a value to n.

\(eqn[n_] := \ 
      P[x]*D2[s[n], x]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  + 
          Q[x]*D1[s[n], x] + R[x]*s[n] == 0\)

eqn[4]

Note that, since  D2[s[n], x]  is of order  n-2, eqn[n] shows the terms up to power n-2 of x.

We want to set the coefficients equal to 0, and get relations among the (Subscript[a, i]). The following command LogicalExpand does this. 
ColumnForm writes the output as a column of equations.

ColumnForm[LogicalExpand[eqn[4]]]

relations[n_] := ColumnForm[LogicalExpand[eqn[n]]]

relations[7]

Let's solve these equations for  (Subscript[a, 5]), (Subscript[a, 4]), (Subscript[a, 3]), and  (Subscript[a, 2])  in terms of  (Subscript[a, 0]) and  (Subscript[a, 1]). We don't need the relations , but  can do this directly with the output of "eqn".

Solve[eqn[5], {(Subscript[a, 5]), (Subscript[a, 4]), (Subscript[a, 3]), (Subscript[a, 2])}]

In general we want to solve eqn[n] for  (Subscript[a, n]),  (Subscript[a, n - 1]),  ... ,  (Subscript[a, 2]) in terms of  (Subscript[a, 0]) and  (Subscript[a, 1]).
Let's store the result in a list.

coeffs[n_] := 
      Simplify[Solve[eqn[n], Table[(Subscript[a, n - i]), {i, 0, n - 2}]]]

coeffs[7]

We will next substitute these values of the coefficients in terms of (Subscript[a, 0]) and (Subscript[a, 1]) into s[n] to get an nth order approximation to the general solution.

y[x_, n_] := s[n] /. First[coeffs[n]]

y[x, 7]

This would be better if we could collect the coefficients of (Subscript[a, 0]) and (Subscript[a, 1]).  In fact there is a command which does this for us, and, if there is a term O[x]^n, it keeps only the polynomial terms of order less than n. This is accomplished by the command Collect.

Here is the final form of y[x,n]

y[x_, n_] := Collect[s[n] /. First[coeffs[n]], {(Subscript[a, 0]), (Subscript[a, 1])}]

y[x, 7]

Recall that (Subscript[a, 0])=y[0], and (Subscript[a, 1])=y'[0].  Then, for example, the approximate solution for the initial value problem
                y'' + xy' + y  =  0,    y[0]=2,   y'[0]=3 
is given by

yy[x_, n_] := y[x, n] /. {(Subscript[a, 0]) -> 2, (Subscript[a, 1]) -> 3}

yy[x, 5]

Finally we find nth order approximations to the usual fundamental set of solutions.

y1[x_, n_] := y[x, n] /. {(Subscript[a, 0]) -> 1, (Subscript[a, 1]) -> 0}

y2[x_, n_] := y[x, n] /. {(Subscript[a, 0]) -> 0, (Subscript[a, 1]) -> 1}

y1[x, 7]

y2[x, 7]

We can now use all these functions as we wish.  For example,

Plot[Evaluate[{y2[x, 4], y2[x, 10], y2[x, 30]}], {x, -5, 5}, 
      AspectRatio -> 1]

What if (Subscript[x, 0])  ≠  0?

As we did in class, we set t = x - (Subscript[x, 0]), and x = t+(Subscript[x, 0]), and solve about the point (Subscript[t, 0])=0.  Afterwards we replace again t by x-(Subscript[x, 0]).

Now the procedure is as follows:

A palette to  enter the data:

({{Equation Data}, {\(Clear[P, Q, R, s, a, eqn, coeffs, x, y, t]; 

            	{P[x_], Q[x_], R[x_], x\_0} = {\[SelectionPlaceholder], 
                \[SelectionPlaceholder], \[SelectionPlaceholder], 
                \[SelectionPlaceholder]}\)}})

Example 2.    y'' - xy = 0  near the point 1.

Thus, we have to replace x by t + 1.  But this can also be implemented in the commands.

\(Clear[P, Q, R, s, a, eqn, coeffs, x, y, t]; 
    
{P[x_], Q[x_], R[x_], x\_0} = {1, 0, \(-x\), 1}\)

s[n_] := (∑_(i = 0)^n) (Subscript[a, i]) (t^i) + O[t]^(n + 1)

\(eqn[n_] := 
	\ \ \ \ 
      P[t + x\_0]*D2[s[n], t] + Q[t + x\_0]*D1[s[n], t] + 
          R[t + x\_0]\ s[n] == 0\)

coeffs[n_] := 
      Simplify[Solve[eqn[n], Table[(Subscript[a, n - i]), {i, 0, n - 2}]]]

y[x_, n_] := 
      Collect[s[n] /. First[coeffs[n]], {(Subscript[a, 0 ]), (Subscript[a, 1])}] /. {
          t -> x - (Subscript[x, 0])}

y1[x_, n_] := y[x, n] /. {(Subscript[a, 0]) -> 1, (Subscript[a, 1]) -> 0}

y2[x_, n_] := y[x, n] /. {(Subscript[a, 0]) -> 0, (Subscript[a, 1]) -> 1}

This is the general procedure for getting the expressions for the solution as a power series.

Notice the little vertical line in the cell brackets for the last seven cells.  We have used the "Cell Properties" command under the Cell menu to make these into "Initialization Cells."  You can evaluate all of them at once by choosing "Evaluate Initialization " from the Kernel/Evaluation menu.

Once we have these cells evaluated, we can use the resulting data as we like.

y[x, 7]

y1[x, 5]

The next plot shows some of the lower term approximations to the solution which has y(1) = 1 and y'(1) = 0.

Plot[
      Evaluate[Table[y1[x, n], {n, 3, 45, 3}], {x, -8, 2},  
        AspectRatio  -> 1]]

Once we have finished the discussion of one equation, we clear, enter new coefficients and a new starting point, and then we can calculate again for the new data...

Problem 5.2/18

(1-x) y'' + x y' - y  =  0, 
        y(0) = -3, y'(0) = 2  near the origin

In this problem we are asked to find the first five nonzero terms in the solution and to plot the four-term and the five-term approximations to the solution in one graph. From the plot we shall estimate the interval in which the four-term approximation is reasonably accurate.

\(Clear[P, Q, R, s, a, eqn, coeffs, x, y, t]; 

    	{P[x_], Q[x_], R[x_], x\_0} = {1 - x, x, \(-1\), 0}\)

After entering the data we simply choose Evaluate Initialization from  the Kernel/Evaluation menu and then we can ask for all expressions -- if we like. 

Although you may not have seen anything happen, the seven evaluation cells were evaluated.

Let's first ask for the general 5th order approximation.

y[x, 5]

The approximation which has 5 nonzero terms is therefore

y[x, 4]

We can also ask for the underlying set of equations or coefficients.

eqn[5]

coeffs[5]

This is the 7th order approximation to the solution which has y(0) = 0 and y'(0) = 1.

y2[x, 7]

Here are the required approximations to the solution of the initial value problem.

y[x, 4] /. {(Subscript[a, 0]) -> -3, (Subscript[a, 1]) -> 2}

y[x, 3] /. {(Subscript[a, 0]) -> -3, (Subscript[a, 1]) -> 2}

We plot both, the 5 term and the 4 term, approximation in one graph.

Plot[
      Evaluate[Table[y[x, n], {n, 3, 4, 1}] /. {(Subscript[a, 0]) -> -3, 
            (Subscript[a, 1]) -> 2}], {x, -2, 2},  AspectRatio  ->  1]

The five-term approximation is the lower one. 
The four-term approximation is reasonably accurate (that is as good as the five-term appr.) for roughly |x| < 0.8.

The above method works even with non-polynomial, but analytic coefficients.

Problem 5.3/11                  y'' + (sin x) y  = 0  about the origin

\(Clear[P, Q, R, s, a, eqn, coeffs, x, y, t]; 

    	{P[x_], Q[x_], R[x_], x\_0} = {1, 0, Sin[x], 0}\)

y[x, 10]

eqn[5]

Mathematica  uses the Taylor series for the sine function and calculates as many coefficients  as desired.

The next plot shows the first approximations to the solution which goes through the origin and has slope 1 at that point., i. e. y(0) = 0 and y'(0) = 1.

Plot[
      Evaluate[Table[y2[x, n], {n, 1, 10, 1}], {x, -3, 3},  
        AspectRatio  -> 1]]