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>
Return to computing page for the second course APMA0330
Return to Mathematica tutorial page for the first course
Return to Mathematica tutorial page for the second course
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
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:
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:
�� 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:
If the generating function Sum[y[n]*z^n , {n,0,infinity}] is known, then (to see the first 10 coefficients)
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 = b; b = c, {n, 20}];
lst
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:
Example. Consider the function
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:
c[5] x^27 + c[6] x^32 + c[7] x^37 + c[8] x^42 + c[9] x^47 +O[x]^48}
coefeq = LogicalExpand[de];
coeffs = Solve[coefeq, Table[c[i], {0, 1, 17}]];
y = y /. coeffs
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)
Then we set n to be 13, the number of terms in the power series:
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
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
Equate the coefficients to 0. Recall that the operator && means And.
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.
To get the solution series, set a[0] equals to zero:
3. Series Solutions
Example: solve the IVP:
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}]]
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:
exactsol=DSolve[{(4-x^2)y'[x]+y[x]==0,y[0]==1},y[x],x]
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]}}