Jacobi Difference Equations
with
Mathematica

The purpose of this notebook is to show how Mathematica can be used to make some calculations with difference equations, in  particular, Jacobi difference equations. All files are supposed to work with Mathematica versions 3.x and 4.x.
It requires the packages DiscreteMath'DiffEqs' and DiscreteMath'JacOp'  which are available via
	 http://www.mat.univie.ac.at/~gerald/ftp/book-jac/DiffEqs.m
	 http://www.mat.univie.ac.at/~gerald/ftp/book-jac/JacDEqs.m
and need 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 available from
	 http://www.mat.univie.ac.at/~gerald/ftp/book-jac/JacDEqs.nb

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

The package DiffEqs and 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}, {Casoratian[(Subscript[f, 1]),..,(Subscript[f, k]),n], Casoratian of a list of sequences}, {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 package JacDEqs and Jacobi difference equations

Now we come to Jacobi difference equations. The necessary commands are provided in the package  DiscreteMath'JacDEqs' which we load first.

Needs["DiscreteMath`JacDEqs`"]

This will also load the package DiscreteMath'DiffEqs' if neccessary. The following comands are provided in the package JacDEqs :

({{JacobiDE[f[n],n], computes the Jacobi difference expression, associated with the sequences  a[n], b[n], of f[n]. }, {JacobiDE[k,f[n],n], applies the Jacobi difference expression  k  times. }, {SolutionJacobi[u], tells Mathematica that u[z,n] solves  JacobiDE[u[z,n],n] = z u[z,n].}, {SolutionJacobi[u,x,y,m], tells Mathematica that u[z,n] solves  JacobiDE[u[z,n],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],n] = 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] satisfies	 JacobiDE[f[n],n] =g[n]}})

The Jacobi difference expression J applied to a sequence gives the following result

JacobiDE[f[n],n]

This can also be written as

JacobiDE[f[n],n] ==-BDiff[a[n]FDiff[f[n],n],n] +(a[n-1] +a[n] +b[n])f[n] // Simplify

or

JacobiDE[f[n],n] ==-FDiff[a[n-1]BDiff[f[n],n],n] +(a[n-1] +a[n] +b[n])f[n] // Simplify

The command

SolutionJacobi[u];

will tell Mathematica that u(z,n) satisfies the Jacobi equation J u(z) = z u(z). Let's see how this works:

JacobiDE[u[z,n],n] == z u[z,n]

We can also define solutions u(z,n,m) satisfying the initial conditions u(z,m,m)=(Subscript[u, 0]) and u(z,m+1,m)=(Subscript[u, 1]) by

Clear[u];
SolutionJacobi[u,(Subscript[u, 0]),(Subscript[u, 1])];

Indeed, we obtain

{JacobiDE[u[z,n,m],n] == z u[z,n,m], u[z,m,m], u[z,m+1,m]}

The solutions c(z,n,m) and s(z,n,m) are predefined by the package

{JacobiDE[c[z,n,m],n] == z c[z,n,m], c[z,m,m], c[z,m+1,m]}

{JacobiDE[s[z,n,m],n] == z s[z,n,m], s[z,m,m], s[z,m+1,m]}

Note that you can tell Mathematica that  u(n) satisfies (J -z)u = 0 using

Clear[u];
SolutionJacobi[u,(z u[#])&];

Indeed,

JacobiDE[u[n],n] == z u[n]

Clear[u];

Simple properties of Jacobi difference equations

Needs["DiscreteMath`JacDEqs`"]

The Wronskian of two solutions (Subscript[W, n])(u(z),v(z)) is defined as  a(n) times the Casoratian of two solutions (Subscript[C, n])(u(z),v(z))

W[u_,v_,n_]:= a[n] Casoratian[u,v,n]

Please note that the argument n tells Mathematica what to consider as index. To compute the value of the Wronskian use

W[u[n],v[n],n] /. n->0

and not

W[u[n],v[n],0]

Defining two solution u(z,n) and v(z,n)

Clear[u,v]; SolutionJacobi[u]; SolutionJacobi[v] ;

we can easily verify that the Wronskian is independent of n

Expand[FDiff[W[u[z,n],v[z,n],n],n]]

Similarly, we can verify Green's formula

(Underoverscript[(∑^⋆), (j=m), n])(f[j] JacobiDE[g[j],j]- g[j]JacobiDE[f[j],j])==W[f[n],g[n],n] -(W[f[n],g[n],n] /. n->m-1)// TestDiff[#,{n,m}]&

We can also define a second solution using the formula

Clear[v];
v[z_,n_]:= (Underoverscript[(∑^⋆), (j=0), (n-1)])((u[z,n])/(a[j]u[z,j]u[z,j+1]));

Indeed, the sequence v(z,n)defined as above satisfies the Jacobi equation

SplitSum[JacobiDE[v[z,n],n] - z v[z,n]] //Factor

The Wronskian of u(z,n) and v(z,n) is one:

SplitSum[W[u[z,n],v[z,n],n]] // Factor

Next recall that the equality u(n)=v(n) holds if equality holds at two consecutive points n=m,m+1and the  Jacobi difference equations applied to both sides gives the same  result, that is, J u(z) = J v(z). We can define a simple test using this fact:

TestJac[eqn_Equal,{n_,m_}]:= Simplify[(JacobiDE[#,n] -z #)& /@ eqn] &&(eqn /. n->m) &&(eqn /. n->m+1);

As an application we verify that any solution can be written as a linear  combination of the two fundamental solutions  c(z,n,m) and s(z,n,m) as follows

SolutionJacobi[u,(Subscript[u, 0]),(Subscript[u, 1])];
u[z,n,m]==( (Subscript[u, 0]) c[z,n,m] +(Subscript[u, 1]) s[z,n,m]) //TestJac[#,{n,m}]&

The solutions  c(z,n,m) and s(z,n,m) have some interesting properties.

s[z,n,m+1] == -((a[m+1])/(a[m]))c[z,n,m]//TestJac[#,{n,m}]&

s[z,n,m-1] == c[z,n,m]+((z-b[m])/(a[m]))s[z,n,m]//TestJac[#,{n,m}]&

c[z,n,m+1] == ((z-b[m+1])/(a[m]))c[z,n,m]+s[z,n,m]//TestJac[#,{n,m}]&

c[z,n,m-1] == -((a[m-1])/(a[m]))s[z,n,m]//TestJac[#,{n,m}]&

Unfortunately, this does not work for the following formula

c[z,m,n] == ((a[n])/(a[m]))s[z,n+1,m]//TestJac[#,{n,m}]&

and we can only check it for specific values of n:

Simplify[(c[z,m,n] == ((a[n])/(a[m]))s[z,n+1,m]) /. n->m+3]

Similarly for the following asymptotics (w=(1/z))

(((∏_(j=m+1)^(n-1))a[j])((w^(n-m-2))/(-a[m]))c[(1/w),n,m]-(1- w(∑_(j=m+2)^(n-1))b[j]+((O[w])^2)) )/. n->m+3

(((∏_(j=m+1)^(n-1))a[j])(w^(n-m-1))s[(1/w),n,m]-(1- w(∑_(j=m+1)^(n-1))b[j]+((O[w])^2)) )/. n->m+3

(((∏_(j=n)^(m-1))a[j])(w^(m-n))c[(1/w),n,m]-(1- w(∑_(j=n+1)^m)b[j]+((O[w])^2)) )/. n->m-3

(((∏_(j=n)^(m-1))a[j])((w^(m-n-1))/(-a[m]))s[(1/w),n,m]-(1- w(∑_(j=n+1)^(m-1))b[j]+((O[w])^2)) )/. n->m-3

Defining the Jacobi matrix

Delta[i_,j_]:=If[Simplify[i==j],1,0];
J[i_,j_]:=a[j] Delta[i-1,j]+a[i] Delta[i+1,j]+b[j]Delta[i,j];JacobiMatrix[m_,n_]:=Table[J[m+i,m+j],{i,n-m-1},{j,n-m-1}]

MatrixForm[JacobiMatrix[m,m+6]]

we can check

(s[z,#,m] ==(( Det[z IdentityMatrix[#-m-1]-JacobiMatrix[m,#]])/((∏_(j=m+1)^(#-1))a[j])) )& [m+4] //Simplify

or compute traces

Tr[M_]:=(∑_(i=1)^(Length[M]))M[[i,i]]; (* Not need for Mathematica 4.x *)

Expand[(Tr[JacobiMatrix[m,#].JacobiMatrix[m,#]]-(2 (∑_(j=m+1)^(#-2))((a[j])^2)+(∑_(j=m+1)^(#-1))((b[j])^2)))&[m+3]]

Let us next investigate the inhomogeneous Jacobi operator. The command

SolutionJacobi[f,g];

tells Mathematica that f(n) satisfies the inhomogenous Jacobi equation J f = g. In fact,

JacobiDE[f[n],n] == g[n]

Next, let us express f(n) in terms of s(z,n,m). We define the kernel

K[z_,n_,j_]:=s[z,n,j]/a[j]

which satisfies

Factor[{K[z,n,n-2],K[z,n,n-1],K[z,n,n],K[z,n,n+1],K[z,n,n+2]}]

and the sequence f(n)

Clear[f];
f[n_]:= (Underoverscript[(∑^⋆), (j=0), n])K[z,n,j]g[j];

SplitSum[JacobiDE[f[n],n]-z f[n]] // Factor

Finally, let us consider the case with a(n)=1/2, b(n)=0.

Unprotect[a,b];
a[n_]=(1/2); b[n_]=0;

Then we have

JacobiDE[(((z+(Sqrt[(z^2)-1])))^n),n]== z (((z+(Sqrt[(z^2)-1])))^n) //Simplify

JacobiDE[(((z-(Sqrt[(z^2)-1])))^n),n]== z (((z-(Sqrt[(z^2)-1])))^n) //Simplify

and

s[z,n,0]==(((((z+(Sqrt[(z^2)-1])))^n)-(((z-(Sqrt[(z^2)-1])))^n))/(2(Sqrt[(z^2)-1]))) // TestJac[#,{n,0}]&

c[z,n,0]== -s[z,n-1,0] // TestJac[#,{n,0}]&

Clear[a,b,u,v,f,g]; Protect[a,b];

Orthogonal Polynomials

Jacobi operators are closely related to orthogonal polynomials and the  moment problem. To see this we first introduce the matrices C(k), D(k), and P(z,k). We use Cm[k] for the matrix C(k) and Cd[k]for its determinant. Similar for the other three matrices.

m[0]=1;
Cd[0]=1;Cd[k_]:=Det[Cm[k]];
Cm[k_?IntegerQ]:=Table[m[i+j],{i,0,k-1},{j,0,k-1}];
Dd[0]=0;Dd[k_]:=Det[Dm[k]];
Dm[k_?IntegerQ]:=Table[If[j==k-1,m[i+j+1],m[i+j]],{i,0,k-1},{j,0,k-1}];
Pd[z_,0]=0;Pd[z_,k_]:=Det[Pm[z,k]];
Pm[z_,k_?IntegerQ]:=Table[If[i==k-1,(z^j),m[i+j]],{i,0,k-1},{j,0,k-1}];

Here is how the matrices look like.

Print[MatrixForm[Cm[3]] ,"   ", MatrixForm[Dm[3]],"   ",  MatrixForm[Pm[z,3]] ]

The determinants of C(k) and D(k) appear as the first two coefficients in the expansion of the determinant  of P(z,k) for large z.

k=3;
Simplify[Series[Evaluate[Pd[1/t,k]],{t,0,-1}]== Cd[k-1] (((1/t))^(k-1))- Dd[k-1] (((1/t))^(k-2))]
k=.;

Unprotect[a,b, s]; (* If JacDEqs is loaded *)

The polynomials s(z,k)

s[z_,k_]:=((Pd[z,k])/((Sqrt[Cd[k]]) (Sqrt[Cd[k-1]])));

are orthogonal with respect to the measure ρ whose moments are m(k)=∫(λ^k)dρ(λ). Moreover, defining

a[k_]:=(((Sqrt[Cd[k+1]]) (Sqrt[Cd[k-1]]))/(Cd[k])); a[0]:=0;
b[k_]:=((Dd[k])/(Cd[k]))-((Dd[k-1])/(Cd[k-1]));

we have that s(z,k) solves the corresponding Jacobi equation.

Simplify[a[k] s[z,k+1]+a[k-1] s[z,k-1]+(b[k]-z) s[z,k]  /. k->3]

Moreover, the polynomials s(z,k) are orthogonal with respect to the scalar product

ScalarProduct[P_,Q_,z_]:= (∑_(i=0)^(Exponent[P,z]))(∑_(j=0)^(Exponent[Q,z]))m[i+j]Coefficient[P,z,i]Coefficient[Q,z,j]/; PolynomialQ[Q,z]&&PolynomialQ[Q,z];

which can be checked as follows

Simplify[ScalarProduct[s[z,j],s[z,k],z] ==If[j==k,1,0]/. {k-> 3, j->3}]

Clear[a,b,m,Cm,Cd,Dm,Dd,Pm,Pd];s[z_,k_]=.; Protect[a,b, s];

Recursions

Now we want to look at  the matrix elements (Subscript[g, j])(n)=⟨(Subscript[δ, n]),(H^j) (Subscript[δ, n])⟩, (Subscript[h, j])(n)=2a(n)⟨(Subscript[δ, n+1]),(H^j) (Subscript[δ, n])⟩ which are important when dealing with high energy expansions for the  resolvent of our Jacobi operator H.

Needs["DiscreteMath`JacDEqs`"]

For computational purpose it is most convenient to define them using the following recursions.

g[0,n_]=1; h[0,n_]=0;
g[j_?IntegerQ,n_]:=g[j,n]=(1/2) (h[j-1,n]+h[j-1,n-1])+b[n] g[j-1,n];
h[j_?IntegerQ,n_]:=h[j,n]=2 ((a[n])^2) (∑_(l=0)^(j-1))g[j-l-1,n] g[l,n+1]-(1/2) (∑_(l=0)^(j-1))h[j-l-1,n] h[l,n];

These definitions can be used to compute arbitrary coefficients  as needed.

Simplify[g[2,n]]

Simplify[h[3,n]]

In addition, we have the third sequence ((Subscript[γ, j])^β)(n)=⟨((Subscript[δ, n+1])+β (Subscript[δ, n])),(H^j) ((Subscript[δ, n+1])+β (Subscript[δ, n]))⟩

γ[j_,n_]:=f[j,n+1]+(β/(a[n])) h[j,n]+(β^2) g[j,n];

But now we want to see that the recursions indeed produce the matrix elements from above. For this purpose we need a delta sequence.

Delta[n_]:=If[n==0,1,0];

Using this we have

g[j,n]==JacobiDE[j,Delta[n-m],m]  /. j->3 /. m->n // Simplify

and

h[j,n]==2 a[n]JacobiDE[j,Delta[n-m],m]  /. j->3 /. m->n+1 // Simplify

Hence the recursions produce indeed the desired matrix elements.
There is  another relation for (Subscript[g, j])and (Subscript[h, j])which can be verified for given j.

Simplify[h[j+1,n] - h[j+1,n-1]==2(((a[n])^2) g[j,n+1]-((a[n-1])^2) g[j,n-1])+b[n]( h[j,n]-h[j,n-1]) /. j->4]

Moreover, there is a recursion involving only (Subscript[g, j]).

Clear[g];
g[0,n_]=1;
g[1,n_]=b[n];
g[2,n_]=((a[-1+n])^2)+((a[n])^2)+((b[n])^2);
g[j_?IntegerQ,n_]:= -((a[n])^2) g[j-2,n+1]+((a[n-1])^2) g[j-2,n-1]-((b[n])^2) g[j-2,n]+2 b[n] g[j-1,n]-(1/2) (∑_(l=0)^(j-2))k[j-l-2,n] k[l,n]+2 ((a[n])^2) ((∑_(l=0)^(j-2))g[j-l-2,n+1] g[l,n]-2 b[n] (∑_(l=0)^(j-3))g[j-l-3,n] g[l,n+1]+((b[n])^2) (∑_(l=0)^(j-4))g[j-l-4,n] g[l,n+1]);
k[0,n_]=-b[n];
k[j_,n_]:=((a[n])^2) g[j-1,n+1]-((a[n-1])^2) g[j-1,n-1]+((b[n])^2) g[j-1,n]-2 b[n] g[j,n]+g[j+1,n];

Again, it can be checked, that it is equivalent to our other one for given  j.

g[j,n]==JacobiDE[j,Delta[n-m],m]  /. j->4 /. m->n // Simplify

In particular, these relations will be used in the sections on the Toda and Kac van Moerbeke lattice.

Some further recursions :

mp[0,n_]=1; mp[1,n_]=b[n+1];
mp[j_?IntegerQ,n_]:=mp[j,n]=b[n+1] mp[j-1,n]+((a[n+2])^2)(∑_(l=0)^(j-2))mp[j-l-2,n]mp[l,n+1];
mm[0,n_]=1; mm[1,n_]=b[n-1];
mm[j_?IntegerQ,n_]:=mp[j,n]=b[n-1] mm[j-1,n]+((a[n-1])^2)(∑_(l=0)^(j-2))mm[j-l-2,n]mm[l,n-1];

Table[mp[j,n],{j,3}]//Simplify

Table[mm[j,n],{j,3}]//Simplify

b[1,n_]=b[n];b[l_?IntegerQ,n_]:=b[l,n]=l g[l,n]-(∑_(j=1)^(l-1))g[l-j,n]b[j,n];

Table[b[j,n],{j,3}]//Simplify

cp[-2,n_]=1;
cp[l_?IntegerQ,n_]:=-(1/(l+2))(∑_(j=1)^(l+2))cp[l-j,n]b[j,n];

cp[-1,n]//Simplify

cp[0,n]//Simplify

Clear[g,h,k,mp,mm,b,cp];

Commutation methods

Single commutation method

In this section we want to consider the single commutation method. As first steps we load our package

Needs["DiscreteMath`JacDEqs`"]

and introduce a few solutions of our original operator.

SolutionJacobi[u]; SolutionJacobi[v];SolutionJacobi[f,g];
W[f__,n_]:=a[n] Casoratian[f,n];

Next, we declare the commuted operator (Subscript[H, σ])

(Subscript[a, σ])[n_]:=-(((Sqrt[-a[n] ])(Sqrt[-a[n+1]]) (Sqrt[u[λ,n]]) (Sqrt[u[λ,n+2]]))/(u[λ,n+1])) ;
(Subscript[b, σ])[n_]:=λ-a[n] (((u[λ,n])/(u[λ,n+1]))+((u[λ,n+1])/(u[λ,n])));
(Subscript[JacobiDE, σ])[f_,n_]:=(Subscript[a, σ])[n] S[f,{n,1}] + (Subscript[a, σ])[n-1] S[f,{n,-1}] + (Subscript[b, σ])[n] f;
(Subscript[W, σ])[f__,n_]:=(Subscript[a, σ])[n] Casoratian[f,n];

and the operator A.

(Subscript[v, σ])[λ_,n_]:=(1/((Sqrt[-a[n]]) (Sqrt[u[λ,n]]) (Sqrt[u[λ,n+1]])));
A[f_,n_]:=(Subscript[v, σ])[λ,n]W[u[λ,n],f,n];

Note that (Subscript[v, σ])is a solution of (Subscript[H, σ]) u = λ u

Simplify[(Subscript[JacobiDE, σ])[(Subscript[v, σ])[λ,n],n]==λ (Subscript[v, σ])[λ,n]]

Now we can verify that A u(z) satisfies (Subscript[H, σ]) (A u(z)) = z A u(z)

Simplify[(Subscript[JacobiDE, σ])[A[u[z,n],n],n]==z A[u[z,n],n]]

and that A f satisfies (Subscript[H, σ]) (A f) = A g.

Simplify[(Subscript[JacobiDE, σ])[A[f[n],n],n]==A[g[n],n]]

Moreover, the following relations for the Wronskians are valid:

Simplify[(Subscript[W, σ])[A[u[z,n],n],A[v[z,n],n],n]== (z-λ)W[u[z,n],v[z,n],n] ]

Simplify[(Subscript[W, σ])[A[u[z,n],n],A[v[zz,n],n],n]==((-a[n] a[n+1] Casoratian[u[λ,n],u[z,n],v[zz,n],n])/(u[λ,n+1]))]

Simplify[(Subscript[W, σ])[A[u[z,n],n],A[v[zz,n],n],n]==
	((z-λ)W[u[z,n],v[zz,n],n] - (((z-zz) v[zz,1+n]W[u[z,n],u[λ,n],n])/(u[λ,1+n])))]

For additional material see the section on the Kac-van Moerbeke lattice below.

Clear[u,v,f,A,Subscript];

Double commutation method

In this section we want to consider the double commutation method. Again we load our package

Needs["DiscreteMath`JacDEqs`"]

and introduce two solutions of our original operator.

SolutionJacobi[u];SolutionJacobi[v];
W[f__,n_]:=a[n] Casoratian[f,n];

Next, we declare the doubly commuted operator (Subscript[H, γ]).

(Subscript[u, γ])[n_]=((u[λ,n])/((Sqrt[(Subscript[c, γ])[n-1] ])(Sqrt[(Subscript[c, γ])[n]])));
(Subscript[c, γ])[n_]=(1/γ) +(Underoverscript[(∑^⋆), (j=-∞), n]) ((u[λ,j])^2);
(Subscript[a, γ])[n_]= a[n](( (Sqrt[(Subscript[c, γ])[n-1] ])(Sqrt[(Subscript[c, γ])[n+1] ]))/((Subscript[c, γ])[n])) ;
(Subscript[b, γ])[n_]= b[n] -  BDiff[((a[n] u[λ,n] u[λ,n+1])/((Subscript[c, γ])[n])),n];
(Subscript[JacobiDE, γ])[f_,n_]:=(Subscript[a, γ])[n] S[f,{n,1}] + (Subscript[a, γ])[n-1] S[f,{n,-1}] + (Subscript[b, γ])[n] f;
(Subscript[W, γ])[f__,n_]:=(Subscript[a, γ])[n] Casoratian[f,n];

As a first fact we note that (Subscript[u, γ]) satisfies (Subscript[H, γ]) (Subscript[u, γ])= λ (Subscript[u, γ]).

Simplify[SplitSum[(Subscript[JacobiDE, γ])[(Subscript[u, γ])[n],n] - λ (Subscript[u, γ])[n]]]

To obtain more solutions we  introduce the operator A

A[f_,{z_,n_}]:=  (( (Sqrt[(Subscript[c, γ])[n] ]))/(Sqrt[(Subscript[c, γ])[n-1] ]))f-(1/(z-λ))(Subscript[u, γ])[n] W[u[λ,n],f,n];

and set

(Subscript[u, γ])[z_,n_]=A[u[z,n],{z,n}];

The sequence (Subscript[u, γ])(z) satisfies (Subscript[H, γ]) (Subscript[u, γ])(z) = z (Subscript[u, γ])(z)

Simplify[ SplitSum[(Subscript[JacobiDE, γ])[(Subscript[u, γ])[z,n],n] - z (Subscript[u, γ])[z,n]] ]

and for (((Subscript[u, γ])(z,n))^2)we have

Together[SplitSum[(((Subscript[u, γ])[z,n])^2)-(((u[z,n])^2)+(1/(((z-λ))^2))BDiff[(((W[ u[λ,n],u[z,n],n])^2)/((Subscript[c, γ])[n])),n])]]

In addition, we note two equivalent form for (Subscript[u, γ])(z)

Factor[SplitSum[(Subscript[u, γ])[z,n]-( (( (Sqrt[(Subscript[c, γ])[n-1] ]))/(Sqrt[(Subscript[c, γ])[n] ]))u[z,n]-(1/(z-λ))(Subscript[u, γ])[n] S[W[u[λ,n],u[z,n],n], {n,-1}])]]

AA[f_,n_]:=(( (Sqrt[(Subscript[c, γ])[n] ]))/(Sqrt[(Subscript[c, γ])[n-1] ]))f - (Subscript[u, γ])[n] (Underoverscript[(∑^⋆), (j=-∞), n])u[λ,j](f /. n->j);

Together[ AA[u[λ,n],n]-(((Subscript[u, γ])[n])/γ)]

and the following relations for the Wronskian:

Simplify[SplitSum[(Subscript[W, γ])[(Subscript[u, γ])[n],A[u[z,n],{z,n}],n]]]

Factor[SplitSum[(Subscript[W, γ])[A[u[z,n],{z,n}],A[v[zz,n],{zz,n}],n] -(W[u[z,n],v[zz,n],n]+(((z-zz)W[ u[λ,n],u[z,n],n]W[ u[λ,n],v[zz,n],n])/((z-λ)(zz-λ)(Subscript[c, γ])[n])))]]

Finally, we turn to the solutions (Subscript[c, β])and (Subscript[s, β]).

u[λ,0]=-Sin[α];u[λ,1]=Cos[α];
SolutionJacobi[sb,-Sin[α],Cos[α],0];
SolutionJacobi[cb,Cos[α],Sin[α],0];

They are transformed into (Subscript[c, (OverHat[β]),γ]) and (Subscript[s, (OverHat[β]),γ])

cp=SplitSum[(Subscript[c, γ])[n+1]]/. n->0;
cm=SplitSum[(Subscript[c, γ])[n-1]]/. n->0;
cn=(Subscript[c, γ])[0];
β=Cot[α];(OverHat[β])=((Sqrt[cm])/(Sqrt[cp]))β;
δ=((a[0])/(cn  (z-λ)))-(β/(cp(1+((OverHat[β])^2)))) ;
(Subscript[sb, γ])[z_,n_]=(((Sqrt[cm])(Sqrt[cp]))/cn)A[sb[z,n],{z,n}];
(Subscript[cb, γ])[z_,n_]=(cp/cn)((1+((OverHat[β])^2))/(1+(β^2)))(A[cb[z,n],{z,n}]-δ A[sb[z,n],{z,n}]);

as can be seen from

Simplify[SplitSum[(((Subscript[sb, γ])[z,n+1])/((Subscript[sb, γ])[z,n])) ] /. n->0]/(OverHat[β])

and

Simplify[SplitSum[(((Subscript[cb, γ])[z,n+1])/((Subscript[cb, γ])[z,n])) ](OverHat[β]) /. n->0]

Clear[u,v,sb,cb,A,AA,β,δ,Subscript];

Dirichlet commutation method

In this section we want to consider the Dirichlet commutation method. As before we load our package

Needs["DiscreteMath`JacDEqs`"]

and introduce two solutions of our original operator.

SolutionJacobi[u]; SolutionJacobi[v];
u[(Subscript[μ, 0]),0]=0;
W[f__,n_]:=a[n] Casoratian[f,n];

Next, we introduce the operator (Subscript[H, (μ,σ)])

(Subscript[W, {μ,σ}])[n_]=((W[u[(Subscript[μ, 0]),n],u[μ,n],n])/(μ-(Subscript[μ, 0])));
(Subscript[a, {μ,σ}])[n_]:=((a[n] (Sqrt[(Subscript[W, {μ,σ}])[n-1] ])(Sqrt[(Subscript[W, {μ,σ}])[n+1]]))/((Subscript[W, {μ,σ}])[n]));
(Subscript[b, {μ,σ}])[n_]:=b[n]-BDiff[((a[n]u[(Subscript[μ, 0]),n] u[μ,n+1])/((Subscript[W, {μ,σ}])[n])),n];
(Subscript[JacobiDE, {μ,σ}])[f_,n_]:=(Subscript[a, {μ,σ}])[n] S[f,{n,1}] + (Subscript[a, {μ,σ}])[n-1] S[f,{n,-1}] + (Subscript[b, {μ,σ}])[n] f;
(Subscript[W, {μ,σ}])[f__,n_]:=(Subscript[a, {μ,σ}])[n]Casoratian[f,n];

We note the alternate expression for (Subscript[b, (μ,σ)])

Simplify[(Subscript[b, {μ,σ}])[n] -( (Subscript[μ, 0])-((a[n]u[μ,n+1] (Subscript[W, {μ,σ}])[n-1])/(u[μ,n] (Subscript[W, {μ,σ}])[n]))-((a[n-1]u[μ,n-1] (Subscript[W, {μ,σ}])[n])/(u[μ,n](Subscript[W, {μ,σ}])[n-1])))]

and (Subscript[W, (μ,σ)])

Simplify[SplitSum[ FDiff[(Underoverscript[(∑^⋆), (j=0), n])u[ (Subscript[μ, 0]),j] u[μ,j]- (Subscript[W, {μ,σ}])[n],n] ]]

Defining the quantities

(Subscript[u, μ])[n_]:=((u[(Subscript[μ, 0]),n])/((Sqrt[(Subscript[W, {μ,σ}])[n-1]])(Sqrt[(Subscript[W, {μ,σ}])[n]])));
(Subscript[u, (Subscript[μ, 0])])[n_]:=((u[μ,n])/((Sqrt[(Subscript[W, {μ,σ}])[n-1]])(Sqrt[ (Subscript[W, {μ,σ}])[n]])));
A[f_,z_,n_]:=((Sqrt[(Subscript[W, {μ,σ}])[n]])/(Sqrt[(Subscript[W, {μ,σ}])[n-1]])) f-(((Subscript[u, (Subscript[μ, 0])])[n] W[u[(Subscript[μ, 0]),n],f,n])/(z-(Subscript[μ, 0])));

we get several solutions of our new operator.

Together[(Subscript[JacobiDE, {μ,σ}])[(Subscript[u, μ])[n],n] -μ (Subscript[u, μ])[n]]

Together[(Subscript[JacobiDE, {μ,σ}])[(Subscript[u, (Subscript[μ, 0])])[n],n] -(Subscript[μ, 0]) (Subscript[u, (Subscript[μ, 0])])[n]]

Together[(Subscript[JacobiDE, {μ,σ}])[A[u[z,n],z,n],n] -z A[u[z,n],z,n]]

In addition, we note

Simplify[((A[u[z,n],z,n])^2)-(((u[z,n])^2)+((u[μ,n])/((((z-(Subscript[μ, 0])))^2)u[(Subscript[μ, 0]),n] ))BDiff[(((W[u[(Subscript[μ, 0]),n],u[z,n],n])^2)/( (Subscript[W, {μ,σ}])[n])),n])]

and the following formulas for Wronskians.

Simplify[(Subscript[W, {μ,σ}])[(Subscript[u, (Subscript[μ, 0])])[n],A[u[z,n],z,n],n]- ((W[u[μ,n],u[z,n],n])/( (Subscript[W, {μ,σ}])[n])) ]

Simplify[(Subscript[W, {μ,σ}])[(Subscript[u, μ])[n],A[u[z,n],z,n],n]-((z-μ)/(z-(Subscript[μ, 0]))) ((W[u[(Subscript[μ, 0]),n],u[z,n],n])/( (Subscript[W, {μ,σ}])[n])) ]

Simplify[(Subscript[W, {μ,σ}])[A[u[z,n],z,n],A[v[zz,n],zz,n],n]-((z-μ)/(z-(Subscript[μ, 0]))) W[u[z,n],v[zz,n],n] -((z-zz)/((z-(Subscript[μ, 0]))(zz-(Subscript[μ, 0])))) ((W[u[μ,n],u[z,n],n]W[u[(Subscript[μ, 0]),n],v[zz,n],n])/( (Subscript[W, {μ,σ}])[n]))]

Finally, we turn to the solutions c(z,n) and s(z,n) which transform into

(Subscript[s, {μ,σ}])[z_,n_]=((Sqrt[ (Subscript[W, {μ,σ}])[m+1]])/(Sqrt[ (Subscript[W, {μ,σ}])[m]]))A[s[z,n,0],z,n] /. m->0;
(Subscript[c, {μ,σ}])[z_,n_]=((z-(Subscript[μ, 0]))/(z-μ)) (A[c[z,n,0],z,n]+(((μ-(Subscript[μ, 0]))u[μ,m+1] A[s[z,n,0],z,n])/((z-(Subscript[μ, 0]))u[μ,m])))/. m->0;

as can be seen from

(Subscript[s, {μ,σ}])[z,0]

Simplify[(Subscript[s, {μ,σ}])[z,1]]

Simplify[(Subscript[c, {μ,σ}])[z,0]]

Simplify[(Subscript[c, {μ,σ}])[z,1]]

Note also

Simplify[(( (Subscript[W, {μ,σ}])[n+1])/( (Subscript[W, {μ,σ}])[n])) /. n->0]

Clear[u,v,A,Subscript];

Toda lattice

Now we turn to the Toda lattice. Our package is needed again.

Needs["DiscreteMath`JacDEqs`"]

To define the Toda hierarchy we need the recursions for (Subscript[g, j]) and (Subscript[h, j]) found before. We only consider the homogeneous case for simplicity.

g[0,n_]=1; h[0,n_]=0;
g[j_?IntegerQ,n_]:=g[j,n]=(1/2) (h[j-1,n]+h[j-1,n-1])+b[n] g[j-1,n];
h[j_?IntegerQ,n_]:=h[j,n]=2 ((a[n])^2) (∑_(l=0)^(j-1))g[j-l-1,n] g[l,n+1]-(1/2) (∑_(l=0)^(j-1))h[j-l-1,n] h[l,n];

This command will display the r'th (homogenous) Toda equation

TL[r_]:=Simplify[{(∂_t)a[n,t]==a[n](g[r+1,n+1]-g[r+1,n]),(∂_t) b[n,t]==h[r+1,n]-h[r+1,n-1]}]/. {a[n_]->a[n,t],b[n_]->b[n,t]}

For example

TL[1]//TableForm

Next, we introduce the second operator in the Lax pair  which will be denoted by TodaP.

TodaP[r_,f_,n_]:=-JacobiDE[r+1,f,n]+(∑_(j=0)^r)(2 a[n]g[j,n]S[JacobiDE[r-j,f,n],{n,1}]-h[j,n] JacobiDE[r-j,f,n]) +g[r+1,n] f;

Here is the first one applied to a sequence f(n)

TodaP[0,f[n],n]

or, alternatively, in matrix form (cut off after 2 terms in both directions).

m=2;
Table[Coefficient[TodaP[0,f[n],n],f[n+i-j]],{i,-m,m},{j,-m,m}] //MatrixForm
m=.;

Similarly, the second one applied to a sequence f(n)

TodaP[1,f[n],n] //Simplify

The following command will display the coefficients in a handy format.

CoeffTodaP[k_]:=Block[{dummy,i,n},
	dummy=TodaP[k,f[n],n];
	TableForm[Table[{f[n+i],":",Simplify[Coefficient[dummy,f[n+i]]]},{i,-k-2,k+2}]]
	];

For example for r=0

CoeffTodaP[0]

and r=1

CoeffTodaP[1]

We can also compute the commutator of H and (Subscript[P, 2r+2])

Commutator[r_,f_,n_]:= JacobiDE[TodaP[r,f,n],n] -  TodaP[r,JacobiDE[f,n],n];

CoeffCommutator[r_]:=Block[{dummy,i,n},
	dummy=Commutator[r,f[n],n];
	TableForm[Table[{f[n+i],":",Simplify[Coefficient[dummy,f[n+i]]]},{i,-2,2}]]
	];

and display its matrix coefficients

CoeffCommutator[0]

Clear[g,h,TL,TodaP];

Kac-van Moerbeke lattice

Finally we come to the Kac van Moerbeke lattice.

Needs["DiscreteMath`JacDEqs`"]

Again we start with the recursions for (Subscript[G, j]) and (Subscript[H, j]) and only consider the homogeneous case for simplicity.

G[0,n_]=1; H[0,n_]=0;
G[j_?IntegerQ,n_]:=G[j,n]=(1/2) (H[j-1,n]+H[j-1,n-2])+(((ρ[n])^2)+((ρ[n-1])^2)) G[j-1,n];
H[j_?IntegerQ,n_]:=H[j,n]=2 (((ρ[n]ρ[n+1]))^2) (∑_(l=0)^(j-1))G[j-l-1,n] G[l,n+2]-(1/2) (∑_(l=0)^(j-1))H[j-l-1,n] H[l,n];

We can easily compute the first few

Table[Simplify[{G[j,n],H[j,n]}],{j,0,1}]//TableForm

This command will display the r'th Kac-van Moerbeke equation

KM[r_]:=Simplify[(∂_t) ρ[n,t]==ρ[n](G[r+1,n+1]-G[r+1,n])] /. ρ[n_]->ρ[n,t];

For example

KM[1]

Next we need two pairs of sequences (Subscript[a, 1]), (Subscript[b, 1])and (Subscript[a, 2]), (Subscript[b, 2])

(Subscript[a, 1])[n_]:= (Subscript[ρ, e])[n] (Subscript[ρ, o])[n]; (Subscript[b, 1])[n_]:= (((Subscript[ρ, e])[n])^2) + (((Subscript[ρ, o])[n-1])^2);
(Subscript[a, 2])[n_]:= (Subscript[ρ, e])[n+1] (Subscript[ρ, o])[n]; (Subscript[b, 2])[n_]:= (((Subscript[ρ, e])[n])^2) + (((Subscript[ρ, o])[n])^2);

plus the corresponding quantities (Subscript[g, k,j]),(Subscript[h, k,j]), k=1,2.

(Subscript[g, k_])[0,n_]=1; (Subscript[h, k_])[0,n_]=0;(Subscript[g, k_])[j_?IntegerQ,n_]:=(Subscript[g, k])[j,n]=(1/2) ((Subscript[h, k])[j-1,n]+(Subscript[h, k])[j-1,n-1])+(Subscript[b, k])[n] (Subscript[g, k])[j-1,n];
(Subscript[h, k_])[j_?IntegerQ,n_]:=(Subscript[h, k])[j,n]=2 (((Subscript[a, k])[n])^2) (∑_(l=0)^(j-1))(Subscript[g, k])[j-l-1,n] (Subscript[g, k])[l,n+1]-(1/2) (∑_(l=0)^(j-1))(Subscript[h, k])[j-l-1,n] (Subscript[h, k])[l,n];

They satisfy several relations which can be checked for given j:

Simplify[((Subscript[g, 2])[j+1,n] ==  (((Subscript[ρ, o])[n])^2)(Subscript[g, 1])[j,n+1]  +(((Subscript[ρ, e])[n])^2)(Subscript[g, 1])[j,n]+ (Subscript[h, 1])[j,n]) /. j->4]

Simplify[((Subscript[h, 2])[j+1,n] -(Subscript[h, 1])[j+1,n]==  (((Subscript[ρ, o])[n])^2)( 2((((Subscript[ρ, e])[n+1])^2)(Subscript[g, 1])[j,n+1]  -(((Subscript[ρ, e])[n])^2)(Subscript[g, 1])[j,n])+((Subscript[h, 1])[j,n+1] -(Subscript[h, 1])[j,n]))) /. j->4]

Simplify[((Subscript[g, 1])[j+1,n] == (((Subscript[ρ, e])[n])^2)(Subscript[g, 2])[j,n]  +(((Subscript[ρ, o])[n-1])^2)(Subscript[g, 2])[j,n-1]+ (Subscript[h, 2])[j,n-1]) /. j->4]

Simplify[((Subscript[h, 2])[j+1,n-1] -(Subscript[h, 1])[j+1,n]==  (((Subscript[ρ, e])[n])^2)( -2((((Subscript[ρ, o])[n])^2)(Subscript[g, 2])[j,n]  -(((Subscript[ρ, o])[n-1])^2)(Subscript[g, 2])[j,n-1])-((Subscript[h, 2])[j,n] -(Subscript[h, 2])[j,n-1]))) /. j->4]

Simplify[((Subscript[h, 2])[j,n-1] -(Subscript[h, 1])[j,n]== 2(((Subscript[ρ, e])[n])^2)( (Subscript[g, 1])[j,n]  -(Subscript[g, 2])[j,n])) /. j->4]

Simplify[((Subscript[h, 1])[j+1,n]-(Subscript[h, 2])[j+1,n-1]==2(  (((Subscript[a, 1])[n])^2)(Subscript[g, 1])[j,n]- (((Subscript[a, 2])[n-1])^2) (Subscript[g, 2])[j,n-1] )-(((Subscript[ρ, o])[n])^2)((Subscript[h, 2])[j,n-1]-(Subscript[h, 1])[j,n] )+ (((Subscript[ρ, e])[n])^2)((Subscript[h, 2])[j,n]-(Subscript[h, 2])[j,n-1])) /.j->4]

Now we introduce the corresponding operators A, (SuperStar[A]), (Subscript[H, k]), (Subscript[P, k,2r+2]), D, (Subscript[Q, 2r+2]), k=1,2.

A[f_,n_]:=(Subscript[ρ, o])[n] S[f,{n,1}] + (Subscript[ρ, e])[n] f;
(SuperStar[A])[f_,n_]:=(Subscript[ρ, o])[n-1] S[f,{n,-1}] + (Subscript[ρ, e])[n] f;
(Subscript[JacobiDE, k_])[f_,n_]:= (Subscript[a, k])[n] S[f,{n,1}] + S[(Subscript[a, k])[n] f,{n,-1}] + (Subscript[b, k])[n] f;
(Subscript[JacobiDE, k_])[r_?IntegerQ,f_,n_]:= Nest[(Subscript[JacobiDE, k])[#,n]&,f,r];
DiracDE[{f_,g_},n_]:= {A[g,n],(SuperStar[A])[f,n]};(Subscript[TodaP, k_])[r_,f_,n_]:=-(Subscript[JacobiDE, k])[r+1,f,n]+(∑_(j=0)^r)(2 (Subscript[a, k])[n](Subscript[g, k])[j,n]S[(Subscript[JacobiDE, k])[r-j,f,n],{n,1}]-(Subscript[h, k])[j,n] (Subscript[JacobiDE, k])[r-j,f,n]) +(Subscript[g, k])[r+1,n] f;
KMQ[r_,{f_,g_},n_]:={(Subscript[TodaP, 1])[r,f,n],(Subscript[TodaP, 2])[r,g,n]};

Again we can use Mathematica to check some simple formulas. For example (SuperStar[A])A= (Subscript[H, 1])

Simplify[(SuperStar[A])[A[f[n],n],n] == (Subscript[JacobiDE, 1])[f[n],n] ]

or A (SuperStar[A])= (Subscript[H, 2])

Simplify[A[(SuperStar[A])[f[n],n],n] ==(Subscript[JacobiDE, 2])[f[n],n] ]

Moreover, we have ((Subscript[H, 2])^r)A = A ((Subscript[H, 1])^r)

Simplify[((Subscript[JacobiDE, 2])[r,A[f[n],n],n]== A[(Subscript[JacobiDE, 1])[r,f[n],n],n])/. r->4]

respectively (SuperStar[A])((Subscript[H, 2])^r)= ((Subscript[H, 1])^r)(SuperStar[A])

Simplify[((SuperStar[A])[(Subscript[JacobiDE, 2])[r,f[n],n],n]==(Subscript[JacobiDE, 1])[r,(SuperStar[A])[f[n],n],n])/. r->4]

and (Subscript[P, 1,2r+2])A - A (Subscript[P, 2,2r+2])= (Subscript[ρ, o])((SuperPlus[(Subscript[g, 1,r+1])])-(Subscript[g, 2,r+1]))(SuperPlus[S]) + (Subscript[ρ, e])((Subscript[g, 2,r+1])-(Subscript[g, 1,r+1]))

Simplify[((Subscript[TodaP, 2])[r,A[f[n],n],n]-A[(Subscript[TodaP, 1])[r,f[n],n],n] ==
(Subscript[ρ, o])[n] ( (Subscript[g, 1])[r+1,n+1]-  (Subscript[g, 2])[r+1,n]) f[n+1] +(Subscript[ρ, e])[n]( (Subscript[g, 2])[r+1,n]-  (Subscript[g, 1])[r+1,n]) f[n])/. r->4]

respectively   (Subscript[P, 1,2r+2])(SuperStar[A])-(SuperStar[A]) (Subscript[P, 2,2r+2]) = -(SuperMinus[(Subscript[ρ, o])])((SuperMinus[(Subscript[g, 2,r+1])])-(Subscript[g, 1,r+1]))(SuperMinus[S]) + (Subscript[ρ, e])((Subscript[g, 2,r+1])-(Subscript[g, 1,r+1]))

Simplify[((Subscript[TodaP, 1])[r,(SuperStar[A])[f[n],n],n]-(SuperStar[A])[(Subscript[TodaP, 2])[r,f[n],n],n] ==
	-(Subscript[ρ, o])[n-1] ( (Subscript[g, 2])[r+1,n-1]-  (Subscript[g, 1])[r+1,n]) f[n-1] +(Subscript[ρ, e])[n]( (Subscript[g, 2])[r+1,n]-  (Subscript[g, 1])[r+1,n]) f[n])/. r->4]

Clear[G,H,KM,Subscript,A,DiracDE,KMQ];