(* :Title: Jacobi Difference Equations *)

(* :Author: Gerald Teschl 
			http://www.mat.univie.ac.at/~gerald/
			Gerald.Teschl@univie.ac.at  *)

(* :Summary:
	This package defines Jacobi operators for the
	use with "DiscreteMath`DiffEqs`" *)

(* :Context: DiscreteMath`JacDEqs` *)

(* :Package Version: 1.0.2 *)

(* :Copyright: Copyright 1998, G. Teschl  *)

(* :History: 
	Version 1.0 - 1.0.2 by G. Teschl *)

(* :Keywords: Difference Equations, Jacobi equation *)

(* :Source: None. *)

(* :Warning: None. *)

(* :Mathematica Version: 3.0 *)

(* :Limitation: *)

(* :Discussion: This package defines the Jacobi
	difference eexpression JacobiDE plus somme additional stuff. *)

BeginPackage["DiscreteMath`JacDEqs`","DiscreteMath`DiffEqs`"]

JacobiDE::usage =
	"JacobiDE[f[n],n] computes the Jacobi difference
	expression of f[n].
	JacobiDE[k,f[n],n] applies the Jacobi difference k times."
SolutionJacobi::usage =
	"SolutionJacobi[u] tells Mathematica that u[z,n] solves
	JacobiDE u[z,n] = z u[z,n],
	SolutionJacobi[u,x,y,m] tells Mathematica that u[z,n] solves
	JacobiDE u[z,n] = z u[z,n] and satisfies the initial
	conditions u[z,m]=x, u[z,m+1]=y,
	SolutionJacobi[u,x,y] tells Mathematica that u[z,n,m] solves
	JacobiDE u[z,n,m] = z u[z,n,m] and satisfies the initial
	conditions u[z,m,m]=x, u[z,m+1,m]=y,
	SolutionJacobi[f,g] tells Mathematica that f[n] solves
	JacobiDE f[n] = g[n]."
a::usage =
	"a[n] coefficient of the Jacobi equation."
b::usage =
	"b[n] coefficient of the Jacobi equation."
c::usage =
	"c[z,n,m] sequence satisfying JacobiDE c[z,n,m] =
	z c[z,n,m] and the initial
	conditions c[z,m,m]=1, c[z,m+1,m]=0."
s::usage =
	"s[z,n,m] sequence satisfying JacobiDE s[z,n,m] =
	z s[z,n,m] and the initial
	conditions s[z,m,m]==0, s[z,m+1,m]=1."

Begin["`Private`"]

(* Jacobi operator *)

JacobiDE[f_,n_]:= a[n] S[f,{n,1}] + S[a[n] f,{n,-1}] + b[n] f;
JacobiDE[k_?IntegerQ,f_,n_]:= Nest[JacobiDE[#,n]&,f,k];
	
(* u satisfies homogenous Jacobi equation *)

	
SolutionJacobi[u_]:=
	(
	u[z_,i_]:=(-a[i-2] u[z,i-2] - b[i-1] u[z,i-1] +
		z u[z,i-1])/a[i-1] /; i>1;
	u[z_,i_]:=(-a[i+1] u[z,i+2] - b[i+1] u[z,i+1] +
		z u[z,i+1])/a[i] /; i<0;
	u[z_,n_+i_]:=(-a[n+i-2] u[z,n+i-2] - b[n+i-1] u[z,n+i-1] +
		z u[z,n+i-1])/a[n+i-1] /; i>1;
	u[z_,n_+i_]:=(-a[n+i+1] u[z,n+i+2] - b[n+i+1] u[z,n+i+1] +
		z u[z,n+i+1])/a[n+i] /; i<0;
	)


SolutionJacobi[u_,x_,y_,m_]:=
	(
	u[z_,n_]:=(-a[n-2] u[z,n-2] -
		b[n-1] u[z,n-1] + z u[z,n-1])/a[n-1] /; n-m>1;
	u[z_,n_]:=(-a[n+1] u[z,n+2] -
		b[n+1] u[z,n+1] + z u[z,n+1])/a[n] /; n-m<0;
	u[z_,n_]:=x /; (n==m);
	u[z_,n_]:=y /; (n==m+1);
	u[z_,n_+i_]:=(-a[n+i-2] u[z,n+i-2] - b[n+i-1]u[z,n+i-1] +
		z u[z,n+i-1])/a[n+i-1] /; i>1;
	u[z_,n_+i_]:=(-a[n+i+1] u[z,n+i+2] - b[n+i+1]u[z,n+i+1] +
        	z u[z,n+i+1])/a[n+i] /; i<0;
	)

SolutionJacobi[u_,x_,y_]:=
	(
	u[z_,n_,m_]:=(-a[n-2] u[z,n-2,m] -
		b[n-1] u[z,n-1,m] + z u[z,n-1,m])/a[n-1] /; n-m>1;
	u[z_,n_,m_]:=(-a[n+1] u[z,n+2,m] -
		b[n+1] u[z,n+1,m] + z u[z,n+1,m])/a[n] /; n-m<0;
	u[z_,n_,m_]:=x /; (n==m);
	u[z_,n_,m_]:=y /; (n==m+1);
	u[z_,n_+i_,m_]:=(-a[n+i-2] u[z,n+i-2,m] - b[n+i-1]u[z,n+i-1,m] +
		z u[z,n+i-1,m])/a[n+i-1] /; i>1;
	u[z_,n_+i_,m_]:=(-a[n+i+1] u[z,n+i+2,m] - b[n+i+1]u[z,n+i+1,m] +
		z u[z,n+i+1,m])/a[n+i] /; i<0;
	)

(* f satisfies inhomogenous Jacobi equation *)

SolutionJacobi[f_,g_]:=
	(
	f[n_+i_]:=(-a[n+i-2] f[n+i-2] - b[n+i-1] f[n+i-1] +
		g[n+i-1])/a[n+i-1] /; i>1;
	f[n_+i_]:=(-a[n+i+1] f[n+i+2] - b[n+i+1] f[n+i+1] +
	 	g[n+i+1])/a[n+i] /; i<0
	)	

(* basic solutions c[z,n,m] and s[z,n,m] *)

SolutionJacobi[s,0,1]; SolutionJacobi[c,1,0];
	
End[]

Protect[{JacobiDE,SolutionJacobi,a,b,c,s}]

EndPackage[]