DiscreteMath DiffEqs

The purpose of this notebook is to show how to use the Mathematica packages DiscreteMath'DiffEqs', which is available via
	 ftp://www.mat.univie.ac.at/pub/teschl/book-jac/DiffEqs.m
and needs to be  stored in the AddOns/Applications/DiscreteMath subfolder of your Mathematica folder. On multi-user systems, you can install an add-on either in the central Mathematica directory (provided you have access to it) or else in your individual user's Mathematica directory (usually ~/.Mathematica/3.0/). (You can also put them  somewhere else, as long as you change the $PATH variable in Mathematica such that the files can be found by Mathematica or evaluate the packages before using this  notebook.)You might want to note that the code is optimized for ease of use and not for speed.
The latest version of this notebook is avaiable from
	 ftp://www.mat.univie.ac.at/pub/teschl/book-jac/DiffEqs.nb

©1999-2013 by Gerald Teschl  <mailto:Gerald.Teschl@univie.ac.at>, <http://www.mat.univie.ac.at/~gerald/>

First order difference equations

We first load the package DiscreteMath'DiffEqs'.

Needs["DiscreteMath`DiffEqs`"]

The package DiscreteMath'DiffEqs' provides some basic commands for dealing with differences and difference  equations.

({{S[f[n],{n,j}], shifts the sequence f[n] by j places. }, {FDiff[f[n],n], (=f[n+1]-f[n] )  forward difference operator}, {BDiff[f[n],n], (= f[n-1]-f[n]) backward difference operator}, {SymbSum[f[j],{j,m,n}],  sum for symbolic manipulations}, {SymbProduct[f[j],{j,m,n}],  product for symbolic manipulations}, {SplitSum[expr], splits sums into smaller parts}, {SplitProduct[expr], splits products into smaller parts}, {SplitAll[expr], splits sums  and products into smaller parts}, {TestDiff[f[n]==g[n],{n,m}], simple test  for equality based on the fact that  f=g if equality holds at one point m and FDiff[f[n]-g[n],n]=0.}})

The main purpose of the built-in Sum  command is to search for a closed form. It is, however, not suitable for  manipulating sums which cannot be brought to a closed form. For example,

Simplify[ (∑_(j=1)^(n+1)) f[j] -  (∑_(j=1)^n)f[j] ]

Moreover,  the built-in Sum is zero if the upper limit is smaller than the lower limit. This can lead  to, at first sight, surprising results:

(Subscript[S, 1])[n_]:=SymbSum[j,{j,1,n}];
(Subscript[S, 2])[n_]:=Sum[j,{j,1,n}];
(Subscript[S, 3])[n_]=Sum[j,{j,1,n}];
{{(Subscript[S, 1])[n],(Subscript[S, 1])[-3]},{(Subscript[S, 2])[n],(Subscript[S, 2])[-3]},{(Subscript[S, 3])[n],(Subscript[S, 3])[-3]}} //MatrixForm

(note that Mathematica assumes n to be positive). Observe that SymbSum does not look for a closed form; which can be obtained by switching to  the built-in Sum:

(Subscript[S, 1])[n] /. SymbSum -> Sum

Observe also  the definition for situations where the upper limit is smaller then the lower limit:

{(Underoverscript[(∑^⋆), (j=0), 1])f[j],(∑_(j=0)^1)f[j],(Underoverscript[(∑^⋆), (j=0), (-3)])f[j],(∑_(j=0)^(-3))f[j] }

Similarly for symbolic products

{(Underoverscript[(∏^⋆), (j=0), 1])f[j],(∏_(j=0)^1)f[j],(Underoverscript[(∏^⋆), (j=0), (-3)])f[j],(∏_(j=0)^(-3))f[j] }

The main purpose of the commands SplitSum and SplitProduct is to simplify expressions involving sums and products (as pointed out  earlier, Mathematica's built-in capabilities in this respect are limited). The command SplitSum will try to break sums into parts and pull out constant factors such that  the rest can be done with built-in commands. In particular, it will split off  extra terms

SplitSum[(Underoverscript[(∑^⋆), (j=0), (n+1)])f[j]]

expand sums and pull out constants

SplitSum[(Underoverscript[(∑^⋆), (j=0), n])(h f[j]+ g[j])]

The command SplitProduct will perform similar operations with products.

SplitProduct[ (Underoverscript[(∏^⋆), (j=0), (n-1)])(g/(f[j]))]

Finally, the command SplitAll also handles nested expressions

SplitAll[(Underoverscript[(∑^⋆), (m=0), n])	(Underoverscript[(∑^⋆), (j=0), m])	(f[j]+(Underoverscript[(∏^⋆), (i=0), j])x g[i])]

We can shift a sequence  f(n) by j places using S[f[n],{n,j}]. For example,

S[f[g[(n^2)]],{n,1}]

or

S[f[n,m],{n,3}]

Similarly, we can compute forward, backward differences using FDiff[f[n],n], BDiff[f[n],n], respectively. For example

{FDiff[f[n],n],BDiff[f[n],n]}

A simple test for equality of sequences is based on the fact that f(n)=0 if f(n) is constant, that is, f(n+1)=f(n) and f((Subscript[n, 0]))=0 at one fixed point.

TestDiff[f[n]==0,{n,(Subscript[n, 0])}]

As an first application we can verify Abels formula (summation by parts)

(Underoverscript[(∑^⋆), (j=m), n])g[j]FDiff[f[j],j]==( g[n]f[n+1] -g[m-1] f[m]+(Underoverscript[(∑^⋆), (j=m), n])BDiff[g[j],j]f[j]) // TestDiff[#,{n,m}]&

Note also the following product rules

FDiff[f[n]g[n],n] == f[n] FDiff[g[n],n] + g[n+1] FDiff[f[n],n] //Simplify

BDiff[f[n]g[n],n] == f[n] BDiff[g[n],n] + g[n-1] BDiff[f[n],n] //Simplify

Next,  the product

P[n_]:=(Underoverscript[(∏^⋆), (j=0), (n-1)])f[j]

satisfies the difference equation P(n+1)=f(n) P(n) plus the initial condition P(0)=1

{P[0],SplitProduct[P[n+1] - f[n] P[n]]}

Moreover, the sum

F[n_]:=(Underoverscript[(∑^⋆), (j=0), (n-1)])f[j]

satisfies the difference equation ∂F(n) = F(n+1)-F(n) = f(n) plus the initial condition F(0)=0

{F[0],SplitSum[FDiff[F[n],n] - f[n]]}

Similarly, the general solution of the difference equation F(n+1) = f(n)F(n) + g(n) with initial condition F(0)=(Subscript[F, 0])is given by

F[n_]:=P[n] ((Subscript[F, 0]) + (Underoverscript[(∑^⋆), (i=0), (n-1)])((g[i])/(P[i+1])))

In fact, we can easily verify this claim:

{F[0],SplitAll[F[n+1] -( f[n] F[n] + g[n])]}

Clear[F,P,Subscript];

The Frobenius method for ODE's

We first load the package DiscreteMath'DiffEqs'.

Needs["DiscreteMath`DiffEqs`"]

In addition the the commands from the previous section, there are two more commands which are usefull when dealing with power series:

({{ShiftPowers[expr,j], shift the index of all sums and products by  j places}, {ShiftIndex[expr,z], shift the index of all sums and products according to powers of z}})

For example, let us solve Bessel's equation

BesselEquation[u_]:=z^2 (∂_(z,z))u+ z (∂_z) u+ ((z^2)-(ν^2) )u;

using the Frobenius method. By the Frobenius theory, there is a solution of the form

u[z_]= (Underoverscript[(∑^⋆), (j=0), ∞])(z^(j+ν)) h[j];

For convennience later on we set negative coefficients equal to 0

h[j_]:=0 /; (j<0);

This results in

BesselEquation[u[z]]

after plugging it into Bessel's equation. Now let us expand all sums and arrange the powers in all resulting sums:

ShiftPowers[ExpandSum[(z^(-ν)) BesselEquation[u[z]] ],z]

Splitting of the extra terms in all sums and simplifying the expression gives the recurrence relation for the coefficients.

SplitSum[%,0] //FactorSum //Simplify

Now you could try to solve it:

RSolve[h[-2+j]+j (j+2 ν) h[j]==0,h[j],j]

In addition, one can also andle some nonlinear equations. Here is a simple Riccati equation:

RiccatiEqn[u_]:=(∂_z) u  + (u^2)-z;

u[z_]:=SymbSum[h[j] (z^j),{j,0,∞}];h[j_]:=0 /; (j<0);

Here is the equation after plugging in our ansatz:

RiccatiEqn[u[z]]

In order to be able to compare coefficients we use the Cauchy product formula

CauchyProduct[%]

The package chooses some unique summation indices to avoid conflicts with existing ones. You can choose more readable ones if you want.

% /. {n1->j,m1->k}

Now we procced as before

FactorSum[SplitSum[ShiftPowers[ExpandSum[%],z],2]] //Simplify

Clear[u,h]