Lattice  integrability  tester

• continuous case,solving of equation  D(y)=b :  
          integration by parts vs .  
              homotopy operator
• discrete case,solving of equation ((T^m)-a)(y)=b  by generalized summation by parts algorithm
• sample code for computation of formal symmetry
• examples
 
17 August 2014  • V . E .  Adler
29 August 2014  • Version for ITP talk

• main definitions

vars[f_] := 
      Union[Cases[f, _u, {0, Infinity}]]
\(\(\(df[f_] := 
      Plus @@ \((\(D[f, #] dif[#] &\) /@ 
            vars[f])\)\)\(
\)\(
\)
    \) (*\ for\ discrete\ *) \)
T[f_, k_] := f /. u[j_] :> u[j + k]
dt[f_, g_] := df[f] /. dif[u[j_]] :> T[g, j]
ords[
        f_] := If[# == {}, {Infinity, -Infinity}, {#[[1, 
                  1]], #[[-1, 1]]}] &[
        vars[f]]


      
Dx[f_] := df[f] /. dif[u[j_]] :> u[j + 1]
Dx[f_, k_] := Nest[Dx, f, k]
ord[
        f_] := If[# == {}, -1, #[[-1]][[1]]] &[
        vars[f]]

• pint[f] : 
        integration by parts .  Returns (y, 
            z)  where  D (y) = f   if   z = 
        0
• H[f] : 
        homotopy operator .  Returns  y  if  E (f) = 
      0

pint[f_] := 
      Module[{a = 0, b = Together[f], c, z, 
          m}, 
While[
And[
            b =!= 0, 
m = ord[b]; 
c = 
              D[b, u[m]]; 
z = 
              Together[D[c, u[m]]]; 
z == 
              0
], 
c = 
            Integrate[c, u[m - 1]]; 
a = 
            Expand[a + c]; 
b = 
            Together[b - Dx[c]]
]; 
{a, 
          Factor[z]}
]

    
iu[f_] := 
      Sum[u[i] Sum[(((-1))^(k - i - 1)) 
              Dx[D[f, u[k]], k - i - 1], {k, i + 1, ord[f]}], {i, 0, 
          ord[f] - 1}]
H[f_] := 
      Integrate[
        iu[f]/λ /. {u[
                k_] :> λ u[k]}, {λ, 0, 
          1}]

(* 
      pint is much faster than H *) 

f = 
        Expand[Dx[
            Sum[u[j]^2 u[j + 1], {j, 0, 15}]]];

    Timing[pint[f]]

    Timing[H[f]]

f = 
        Dx[u[10]^2 u[2] + u[4]^2 Exp[u[7]] + 
            g[u[8]] u[4]^2];
Timing[pint[f]]
Timing[H[f]]

(* 
      and sometimes H does not work *) 

f = u[1]/u[0]^2;

    pint[f]

    H[f]

• generalized summation by parts algorithm

psum[m_, a_, b_] /; 
        m > 0 := 
Block[{$RecursionLimit = Infinity}, 
        psu[m, a, b]]

    
psu[m_, aa_, bb_] := 
      Module[
{a = Together[aa], 
          b = Together[bb], 
A, B, p1, p2, q1, q2, r, 
          y}, 
q2 = ords[a]; q1 = q2[[1]]; 
        q2 = q2[[2]]; 
p2 = ords[b]; p1 = p2[[1]]; 
        p2 = p2[[2]]; 

Catch[
If[a === 0, 
            Throw[{T[b, -m], 
                0}]]; 

 (* 
            Case a = const *) 
If[
            q1 == Infinity, 
If[
              p2 < p1 + m, 
y = 
                If[a === 1, const, b/(1 - a)]; 
Throw[
                Together[{y, 
                    T[y, m] - a y - 
                      b}]]
]; 
B = 
              Together[
                D[b, u[p2]]]; 
r = ords[
                  B][[1]]; 
If[r < p1 + m, 
              Throw[{0, D[B, u[r]]}]]; 
B = 
              Integrate[B, u[p2]]; 
Throw[
              psu[m, a, b - B + a T[B, -m]] + {T[B, -m], 
                  0}]
]; 

 (* Reduction of a *) 
A = 
            Together[
              D[a, u[q1]]/a]; 
r = ords[
                A][[2]]; 
If[
            And[q1 ≤ q2 - m, 
              r ≤ q2 - m], 
A = 
              Exp[Together[Integrate[A, u[q1]]]]; 
Throw[
              Together[
                psu[m, a T[A, m]/A, b T[A, m]]/
                  A]]
]; 

 (* Reduction of b *) 
If[
            p2 > q2, 
B = 
              Integrate[D[b, u[p2]], u[p2]]; 
Throw[
              psu[m, a, b - B + a T[B, -m]] + {T[B, -m], 
                  0}]
]; 
 
If[
            p1 < q1, 
B = 
              Together[
                D[b, u[p1]]/a]; 
r = ords[
                  B][[2]]; 
If[r > q2 - m, 
              Throw[{0, D[B, u[r]]}]]; 
B = 
              Integrate[B, u[p1]]; 
Throw[
              psu[m, a, b + T[B, m] - a B] - {B, 
                  0}]
]; 

 (* 
            Solving of a linear system *) 
y = 
            If[q1 ≤ 
                q2 - m, 
-(D[b, u[q1]] D[a, u[r2]] - 
                      D[b, u[q1], u[r2]] a)/
(D[a, 
                        u[q1]] D[a, u[r2]] - 
                    D[a, u[q1], u[r2]] a), 
b/(1 - 
                    a)]; 
Throw[
            Together[{y, 
                T[y, m] - a y - 
                  b}]]
]
]

• Discrete homotopy operator (for T - 1)

diu[f_] := 
      Sum[u[i] D[Sum[T[f, -k + i], {k, i + 1, ords[f][[2]]}], 
            u[i]], {i, 
          0, ords[f][[2]] - 1}]
dh[f_] := 
      Integrate[
        diu[f]/λ /. {u[
                k_] :> λ u[k]}, {λ, 0, 
          1}]
dH[f_] := T[dh[T[f, -#]], #] &[ords[
            f][[1]]]

f = 
        Expand[(T[#, 1] - #) &[
            Product[u[j], {j, 0, 
                100}]]];

Timing[psum[1, 1, f]]
Timing[dH[f]]

f = 
        Expand[(T[#, 1] - #) &[
            Sum[j u[j], {j, 0, 100}]]];
Timing[psum[1, 1, f]]
Timing[dH[f]]

• computation of formal symmetry

Clear[m, a, b, c, f, F, 
      g]

dt[f_] := dt[f, F]
f[j_] := D[F, u[j]]
g[j_] /; j > 0 := f[j]
a[j_] := T[f[m], j]/f[m]
b[j_] := 
      1/f[m] (dt[g[j + m]] - 
            Sum[f[s] T[g[j + m - s], s] - 
g[j + m - s] 
                  T[f[s], j + m - s], {s, -m, 
                m - 1}])

mytest[k_] := 
      Do[
ps = 
          Factor[psum[m, a[j], b[j]]]; 
obst = 
          ps[[2]]; 
Print[obst]; 
If[
          Not[obst === 0], Break[]]; 
g[j] = 
          ps[[1]] /. const -> c[j], 
{j, 
          0, -k, -1}]

• Example : 
        Volterra and Bogoyavlensky lattices

F = u[0] (u[1] - u[-1]);
F2 = 
        u[0] (u[1] (u[2] + u[1] + u[0]) - 
              u[-1] (u[0] + u[-1] + 
                    u[-2]));
Expand[dt[F, F2] - dt[F2, F]]

    
m = 1;
mytest[4]
Table[
      Factor[g[j] /. c[j_] :> 0], {j, 
        m, -4, -1}]

    
m = 2;
F = F2;
mytest[4]
Table[
      Factor[g[j] /. c[j_] :> 0], {j, 
        m, -4, -1}]

m = 2;
F = 
        u[0] (u[2] + u[1] - u[-1] - 
              u[-2]);
mytest[6]
Table[
      Factor[g[j] /. c[j_] :> 0], {j, 
        m, -6, -1}]

• Example : sample classification

m = 2;
F = 
        u[0] (u[2] + k1 u[1] + k2 u[0] + k3 u[-1] + 
              k4 u[-2]);
mytest[6]
Collect[-obst, _u]

F = 
        u[0] (u[2] + k1 u[1] - (1 + k4) u[0] - k1 u[-1] + 
              k4 u[-2]);
mytest[6]
Collect[-obst, _u]

m = 2;
F = 
        u[0] (u[2] - (1 + k4) u[0] + 
              k4 u[-2]);
mytest[6]
obst

• Example : 
        Lattices with implicit rhs

mytest1[k_, ru_] := 
      Do[
ps = 
          Factor[psum[m, a[j] /. ru, b[j] /. ru]]; 
obst = 
          ps[[2]]; 
Print[obst]; 
If[
          Not[obst === 0], Break[]]; 
g[j] = 
          ps[[1]] /. const -> c[j], 
{j, 
          0, -k, -1}]

    
m = 1;
F = 
        h[u[1] - u[0]] + h[u[0] - u[-1]];
P[x_] := α x^2 + β x + γ
mytest1[
      4, {h'[x_] :> 
          P[h[x]], h''[x_] :> P'[h[x]] 
            P[h[x]]}]

m = 1;
F = 
        1/(Y[u[1] + u[0]] - 
              Y[u[0] + u[-1]]);
Q[x_] := α x^4 + β x^3 + γ x^2 + δ x + ϵ
mytest1[
      1, {Y'[x_] :> 
          Q[Y[x]], Y''[x_] :> Q'[Y[x]] 
            Q[Y[x]]}]

• Example : Modified B and V

m = 2;
F = 
        u[0] (u[2] u[1] - 
              u[-1] u[-2]);
mytest[6]

c[-1] := 0
mytest[6]

c[j_] /; OddQ[j] := 0
mytest[6]
Clear[c]

F = (u[0]^2 + 
              1) ((u[1]^2 + 1) (u[2] - 
                    u[0]) + 
(u[-1]^2 + 
                    1) (u[0] - 
                    u[-2]));
mytest[4]

c[j_] /; OddQ[j] := 0
mytest[4]
Clear[c]

F = (u[0]^2 + 
              1) ((u[1]^2 + 1) (u[2] + 
                    u[0]) - 
(u[-1]^2 + 
                    1) (u[0] + 
                    u[-2]));
mytest[4]

ru2 = {u[j_] :> (((-1))^j) 
              u[j]};
ff[j_] := (((-1))^j) T[F, j] /. 
        ru2

    
ff[0];
Expand[{ff[1] - T[ff[0], 1], 
        ff[2] - T[ff[0], 2]}]

F1 = (u[0]^2 + 1) (u[1] - u[-1]);
Expand[dt[F, F1] - dt[F1, F]]

• Example : dSK and dKK

m = 2;
F = u[0]^2 (u[2] u[1] - u[-1] u[-2]) - 
          u[0] (u[1] - u[-1]);
mytest[4]

m = 3;
F = 
        u[0] (u[3] u[1] + u[2] u[0] - u[0] u[-2] - 
              u[-1] u[-3] - u[2] - u[1] + u[-1] + 
              u[-2]);
mytest[4]

m = 3;
F = 
        u[0] (u[3] u[2] + u[1] u[0] - u[0] u[-1] - 
              u[-2] u[-3] - u[2] - u[1] + u[-1] + 
              u[-2]);
mytest[4]

m = 5;
F = u[0] (T[#, 3] + T[#, 2] - T[#, -2] - T[#, -3]) &[
        u[2] u[0] u[-2] - u[1] u[-1]]
mytest[6]

m = 5;
F = 
      u[0] ((T[#, 8] + T[#, 5] - T[#, 3] - #) &[
                  u[-3] u[-4] 
                    u[-5]] - 
(T[#, 7] + T[#, 5] - 
                    T[#, 2] - #) &[
            u[-3] u[-4]])
mytest[6]