# Series and Recurrences

This tutorial contains many Mathematica scripts. You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to <Vladimir_Dobrushkin@brown.edu>

# Part V. Series and Recurrences

Recurrences
Sturm-Liouville problem
Heat transfer equations
Wave equations
Laplace equation
Applications

# 5.1. Recurrences

## I. Case Sensitivity

To solve the linear constant coefficient recurrence (also called the difference equation) of the second order

y_{n+2} =6y_{n+1} -y_n               subject to y_0 = 0,   y_1 = 1

you may try some options:

LinearRecurrence[{6, -1}, {0, 1}, 10]

Here the first pair represents the coefficients of the difference equation of order 2; next pair tells Mathematica the first two initial values of the recurrence. The final number (10) is the number of terms the user wants to see as output. Another option:

Transpose[
�� NestList[Flatten[{Rest[#], ListCorrelate[{-1, 6}, #]}] &, {0, 1},
���� 10]][[1]]

If you know that its solution is represented through the Chebyshev polynomial of the second kind at the point x=3, then type:

a[n_] := ChebyshevU[n - 1, 3]

If the generating function Sum[y[n]*z^n , {n,0,infinity}] is known, then (to see the first 10 coefficients)

CoefficientList[Series[z/(1-6 z + z^2), {z, 0, 10}],z]

Consider another recurrence that generates the Pell numbers:

y_{n+2} = 2 y_{n+1} + y_n ,         y_0 =0, y_1 =1,     n=0,1,2,.... .

Then these numbers can be generated as follows:

a = 1; b = 0; c = 0; lst = {b}; Do[c = a + b + c; AppendTo[lst, c];
�� a = b; b = c, {n, 20}];
lst
CoefficientList[Series[x/(1 - 2*x - x^2), {x, 0, 20}], x] (* 20 numbers *)
Expand[Table[((1 + Sqrt[2])^n - (1 - Sqrt[2])^n)/(2 Sqrt[2]), {n, 0, 20}]]

Recurrences, although a very tedious computation method by hand, is very simple to do in Mathematica. The best way to learn how to do recurrences in Mathematica are by examples, and a perfect example for this topic is the Fibonacci sequence.

# 5.2. Generating Functions

Recurrences, although a very tedious computation method by hand, is very simple to do in Mathematica. The best way to learn how to do recurrences in Mathematica are by examples, and a perfect example for this topic is the Fibonacci sequence.

The process of shifting (or slipping) indices can be obtained by the following Mathematica code:

SlipIndices[Op_[expr_, {k_, a_, b_}], d_] := Op[(expr /. k -> k + d), {k, a - d, b - d}]
As an example, consider shifting indices in the sum:
SlipIndices[Sum[k*(k - 1)*c[k]*x^(k - 2), {k, 2, Infinity}], 2]
Out[2]= $Sum], {k 0, \[Infinity]} (1 + k)*(2 + k)*x^k * c[2 + k] Example. Consider the function \[ f(x) = 2\,\sin (x) -3\,\arcsin (x) .$
f[x_] := 2*Sin[x] + 2*ArcSin[x]
We expand this function into the Taylor series about 0:
ser = Series[f[x], {x, 0, 14}]
$4x + \frac{x^5}{6} -\frac{4\, x^7}{45} + \frac{5513\, x^9}{90720} + \frac{2537\, x^{11}}{56700} - \frac{4156001\,x^{13}}{119750400} + O\left( x^{15} \right)$
The 0 indicates that the series is centered about 0 (such series are often called Maclaurin series) and the 14 specifies the highest power sought. The Series command gives its output in a special form, with a big-Oh error term. If only the partial sum of the series is wanted, then use Normal:
Normal[ser]
There is a coincidental cancellation of the third-degree terms when sine and arcsine are added, and the function is close to being linear. The contribution of the fifth- and higher-order terms, whose coefficients are rather small on the interval [-1, 1], is very small.

# 5.3. Series Solutions for the first Order Equations

## I. How to

Consider the initial value problem for the Abel equation

y' = y^3 +x,     y(0) =0

We find its power series expansion using the following Mathematica commands:

y = Sum[c[i]*x^{2 + 5*i}, {i, 0, 9}] + O[x]^(48)
Out[1]= { c[0] x^2 + c[1] x^7 + c[2] x^12 + c[3] x^17 + c[4] x^22 +
c[5] x^27 + c[6] x^32 + c[7] x^37 + c[8] x^42 + c[9] x^47 +O[x]^48}
de = D[y, {x, 1}] - y^3 - x == 0;
coefeq = LogicalExpand[de];
coeffs = Solve[coefeq, Table[c[i], {0, 1, 17}]];
y = y /. coeffs
Out[5]= {{ x^2 /2 + x^7 /56 + x^12 / 896 + 33 x^17 / 426296 + 1475
x^22 / 262721536 + 879 x^27 / 2101772288 + }}

This problem can be solved in slightly different way. First, we define the differential operator (non-linear because the given equation is non-linear)

L[x_, y_] = y'[x] - x - (y[x])^3
Out[1]= -x - y[x]^3 + Derivative[1][y][x]

Then we set n to be 13, the number of terms in the power series:

n = 13
Out[2]= 13

The next cell says to create a sum of n terms and effectively turn it into a series saying that the terms beyond n are indefinite. O[x]^(n+1) indicates that we know nothing about terms of order n+1 and beyond. So we define the (finite) series

y[x_] = Sum[a[i]*x^i, {i, 0, n}] + O[x]^(n + 1)
Out[3]= a[0] + a[1] x + a[2] x^2 + a[3] x^3 + a[4] x^4 + a[5] x^5 +
a[6] x^6 + a[7] x^7 + a[8] x^8 + a[9] x^9 + a[10] x^10 + a[11] x^11 +
a[12] x^12 + a[13] x^13 + O[x]^14

Substitute this expression into the differential operator

LHS = Collect[L[x, y], x] + O[x]^n

Equate the coefficients to 0. Recall that the operator && means And.

sys = LogicalExpand[LHS == 0]

Find the coefficients a[i] in the poer series of y[x] using Reverse[Table[a[i],{i,0,n}]]. The latter makes the form of the
coefficients agree with what you would find by hand: the next term is expressed through the previous terms.

coeff = Solve[sys, Reverse[Table[a[j], {j, 0, n}]]]
Solve::svars: Equations may not give solutions for all "solve" variables. >>

To get the solution series, set a[0] equals to zero:

s[x_] = y[x] /. coeff[[1]] /. a[0] -> 0
Out[7]= x^2/2 + x^7/56 + x^12/896 +O[x]^14

# 3. Series Solutions

Example: solve the IVP:

$(4-x^2 )y'(x)+y(x)=0, \qquad y(0)=1.$
Clear[a,n]
a[n_]:=a[n]=((n-2) a[n-2] -a[n-1])/4/n
a[0] = a0
a[1] = -a0/4
TableForm[Table[{n,a[n]},{n,0,5}]]
Out[21]/TableForm=
0 a0
1 -a0/4
2 a0/32
3 -3 a0/128
4 11 a0/2048
5 -31 a0/8192

When we apply the initial condition, we set a0=1 in the general solution.

Now we solve the initial value problem exactly:

Clear[x,y]
exactsol=DSolve[{(4-x^2)y'[x]+y[x]==0,y[0]==1},y[x],x]
Out[63]= {{y[x]-> (2-x)^(1/4)/(2+x)^(1/4) }}
Note that the exact solution is |2-x|^(1/4)/|2+x|^(1/4) with absolute values instead of parenthethis.

formula=Simplify[(2-x)^(1/4)/(2+x)^(1/4)]
Out[64]= (2 - x)^(1/4)/(2 + x)^(1/4)
var=Table[D[y[x],i]/.x->0,{i,1,11}]
????sols=Solve[sysofeqs,vars]
sersol=Series[y[x],{x,0,11}/.sols[[1]]]

Legendre's equation
genform =
Solve[(n + 2)*(n + 1)*a[n + 2] + (-n*(n - 1) - 2 n + k*(k + 1))*
a[n] == 0, a[n + 2]]
Out[8]= {{a[2 + n] -> -(((k + k^2 - n - n^2) a[n])/((1 + n) (2 + n)))}}
Factor[genform[[1, 1, 2]]]
Out[9]= ((-k + n) (1 + k + n) a[n])/((1 + n) (2 + n))
We obtain a formula for a[n] by replacing each occurrence of n in a[n+2] by n-2.
genform[[1, 1, 2]] /. n -> n - 2
Out[10]= -(((2 + k + k^2 - (-2 + n)^2 - n) a[-2 + n])/((-1 + n) n))
a[n_] := -(((2 + k + k^2 - (-2 + n)^2 - n) a[-2 + n])/((-1 + n) n))
DSolve[(1 - x^2)*y''[x] - 2 x y'[x] + k*(k + 1)*y[x] == 0, y[x], x]
Out[12]= {{y[x] -> C[1] LegendreP[k, x] + C[2] LegendreQ[k, x]}}

# 5.4. Series Solutions for the Second Order Equations

## I. How to determine a polynomial approximation to solution

Set k equal to desired x_0 value (maintain a positive radius of convergence, x-x_0 must be greater than 0)
Set n equal to the highest power term desired in the power series
Set yInitial equal to the value of y when x equals 0

a[1] represents the coefficient in front of the x^1 term
a[2] represents the coefficient in front of the x^2 term
etc.

Note: In order to solve for instances when k (x_0) does not equal 0, I suggest you reduce the n term to something below 5 in order to prevent a significantly long run time

In[1]:= Clear[yInitial];
Clear[f];
Clear[k];

k=0;
n=8;
yInitial=2;

f=yInitial+Sum[a[i] (x-k)^i,{i,n}]+O[x]^(n+1);
D[f,x]-f^2==x^2;
LogicalExpand[%];
Solve[%]

Out[6]= {{a[8]->73029/140,a[7]->81989/315,a[6]->1948/15,a[5]->324/5,a[4]->97/3,a[3]->49/3,a[2]->8,a[1]->4}}

Airy functions are solutions to the Airy differential equation
y''[x] - x y[x] == 0
There are two linearly independent solutions, called by mathematica as
AiryAi[x] and AiryBi[x].

Series[AiryAi[x], {x, 0, 9}]:
N[%]
Out[2]= 0.355028 -0.258819 x + 0.0591713 x^3 - 0.0215683 x^4 +
0.00197238 x^6 - 0.000513531 x^7 + 0.0000273941 x^9 - 5.7059 x
10^^{-6} + )(x)^10

Consider Legendre's equation
(1-x^2)y'' -2xy' + n(n+1) y =0
when n=4. Legendre's equation comes up in many physical situations
involving spherical symmetry. First we set up a series to work with:

y = Sum[c[i] x^i, {i, 0, 6}] + O[x]^7
Out[1]= c[0] + c[1] x + c[2] x^2 + c[3] x^3 + c[4] x^4 + c[5] x^5 + c[6] x^6 + O[x]^7

To see that our input is correct we type

OutputForm[y]//Normal

We are going to insert this series into the differential equation

de = (1 -x^2) D[y, {x,2}] - 2 x D[y,x] + 20 y == 0
Out[2]= (20 c[0] + 2 c[2]) + (18 c[1] + 6 c[3]) x + (14 c[2] + 12 c[4])
x^2 + (8 c[3] + 20 c[5]) x^3 + 30 c[6] x^4 + O[x]^5 == 0
Then we get the coefficients
coeffeqns = LogicalExpand[de]
Out[3]= 20 c[0] + 2 c[2] == 0 && 18 c[1] + 6 c[3] == 0 &&
14 c[2] + 12 c[4] == 0 && 8 c[3] + 20 c[5] == 0 && 30 c[6] == 0
and express through first two c[0] and c[1]:
solvedcoeffs = Solve[ coeffeqns, Table[ c[i], {i,2,12}]]

Solve::svars: Equations may not give solutions for all "solve" variables. >>

Out[4]= {{c[2] -> -10 c[0], c[3] -> -3 c[1], c[4] -> (35 c[0])/3, c[5] -> (6 c[1])/5, c[6] -> 0}}
Put these back into the series expansion.
y = y /. solvedcoeffs
Out[5]= { c[0] + c[1] x -10 c[0] x^2 - 3 c[1] x^3 + (35/3) c[0] x^4 + (6/5) c[1] x^5 + O[x]^7 }
Extract out the two linearly independent solutions
Coefficient[ y, c[0]]
Out[6]= {1 - 10 x^2 + (35 x^4)/3}
Coefficient[ y, c[1]]
Out[7]= {x - 3 x^3 + (6 x^5)/5}
To check our answer, we use Mathematica again because it knows Legendre:
LegendreP[4,x]
Out[8]= 1/8 (3 - 30 x^2 + 35 x^4)
This is different from the polynomial obtained earlier, but only by a constant factor. The constant factor is used to set the normailization for the Legendre polynomials.

# 1.5.5. Series Solutions

Example: solve the IVP:
(4-x^2 )y'[x]+y[x]==0, y[0]==1

Clear[a,n]
a[n_]:=a[n]=((n-2) a[n-2] -a[n-1])/4/n
a[0] = a0
a[1] = -a0/4
TableForm[Table[{n,a[n]},{n,0,5}]]
Out[21]/TableForm=
0 a0
1 -a0/4
2 a0/32
3 -3 a0/128
4 11 a0/2048
5 -31 a0/8192

When we apply the initial condition, we set a0=1 in the general solution.

Now we solve the initial value problem exactly:

Clear[x,y]
exactsol=DSolve[{(4-x^2)y'[x]+y[x]==0,y[0]==1},y[x],x]
Out[63]= {{y[x]-> (2-x)^(1/4)/(2+x)^(1/4) }}
Note that the exact solution is |2-x|^(1/4)/|2+x|^(1/4) with absolute values instead of parenthethis.

formula=Simplify[(2-x)^(1/4)/(2+x)^(1/4)]

Out[64]= (2 - x)^(1/4)/(2 + x)^(1/4)
var=Table[D[y[x],i]/.x->0,{i,1,11}]
????sols=Solve[sysofeqs,vars]
sersol=Series[y[x],{x,0,11}/.sols[[1]]]

Legendre's equation

genform =
Solve[(n + 2)*(n + 1)*a[n + 2] + (-n*(n - 1) - 2 n + k*(k + 1))*
a[n] == 0, a[n + 2]]
Out[8]= {{a[2 + n] -> -(((k + k^2 - n - n^2) a[n])/((1 + n) (2 + n)))}}
Factor[genform[[1, 1, 2]]]
Out[9]= ((-k + n) (1 + k + n) a[n])/((1 + n) (2 + n))
We obtain a formula for a[n] by replacing each occurrence of n in a[n+2] by n-2.
genform[[1, 1, 2]] /. n -> n - 2
Out[10]= -(((2 + k + k^2 - (-2 + n)^2 - n) a[-2 + n])/((-1 + n) n))
a[n_] := -(((2 + k + k^2 - (-2 + n)^2 - n) a[-2 + n])/((-1 + n) n))
DSolve[(1 - x^2)*y''[x] - 2 x y'[x] + k*(k + 1)*y[x] == 0, y[x], x]
Out[12]= {{y[x] -> C[1] LegendreP[k, x] + C[2] LegendreQ[k, x]}}

# 1.5.5a. Bessel functions

phi5=Normal[Series[BesselJ[1,x],{x,0,5}]]
Out[1]= x/2-x3/16+x5/384
psi5=FullSimplify[phi5*Integrate[phi5^-2 *(1/x),x]]
Out[2]=
(1/4608)x (192 - 24 x^2 + x^4) ((576 (-8 + x^2))/( x^2 (192 - 24 x^2 + x^4)) +
12 Log[x] + (-3 - 3 I Sqrt[3]) Log[-12 - 4 I Sqrt[3] + x^2] +
3 I (I + Sqrt[3]) Log[-12 + 4 I Sqrt[3] + x^2])

Plot[{psi5, BesselY[1, x]}, {x, 0.1, 6},
PlotLegends -> {psi, BesselY}, PlotStyle -> {Blue, Black}]
Plot[{phi5, BesselJ[1, x]}, {x, 0.1, 6},
PlotLegends -> {phi, BesselJ}, PlotStyle -> {Blue, Black}]

# 5.6. Motivated Examples from Music

## I. How to

As an example, note A above middle C is the note on which most tunings of instruments is based: 440 Hz. In radians, we could describe this note A with

sin(440*(2*\pi t)) = sin (880*pi *t)

Musical tones are restricted to a limited set of frequencies we call
notes. How we judge combinations of notes is partly subjective, but
most people would agree on some simple basics. For example, notes C
and G blend nicely together and produce a stable harmony. Change note
G to F# and now the combination of C and F# is unstable, producing
tension. Adjacent note combinations such as C and C# simply clash and
are unpleasant to the ear. Can these effects be seen mathematically?
That is the purpose for graphing note combinations here.

For purposes of this project musical notes are assigned relative
frequencies with the lowest note, C, having a frequency of 1 cycle per
second, or 1Hz. (Hz is a "Hertz.") While these are not true
frequencies, the math is easier and the comparative results the same
as if true frequencies were used. (Actual middle C has a frequency
close to 262 Hz.) Subsequent notes are tuned by the well-tempered
system of tuning in use since 1720. The ratio of frequencies of
consecutive chromatic notes is the twelfth root of 2, approximately
1.05946. This means that multiplying the frequency of a note by
1.05946 will yield the frequency of the next note above it on the
chromatic (all-inclusive) scale.

For a chart of relative frequencies over two octaves (with lowest C