
(* Last updated: January 31, 2007 at 23:50 by Willy Hereman at home *)

(****************************************************************)
(*                                                              *)
(*      *** M A T H E M A T I C A  P R O G R A M ***            *)
(*                                                              *)
(*   COMPUTER ASSISTED CALCULATION OF SOLITON SOLUTIONS         *)
(*      OF NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS             *)
(*                 WITH HIROTA'S METHOD                         *)
(*                                                              *)
(*  program name: hirota.m                                      *)
(*                                                              *)
(*  purpose: calculation of the DISPERSION LAW, verIfication    *)
(*  of the CONDITIONS for the EXISTENCE of two, three, and four *)
(*  soliton solutions, and the Explicit SOLITON SOLUTIONS       *)
(*  (up to three-soliton solutions).                            *)
(*                                                              *)
(*  input: single bilinear equation of KdV type (type I),       *)
(*         and a set of parameters to control the execution.    *)
(*                                                              *)
(*  output: the function f and the coupling coefficients        *)
(*          Note that f is related to u via a transformation    *)
(*          of type u = c*d^2 log(f)/dx^2, where c is a         *)
(*          scaling constant.                                   *)
(*                                                              *)
(*  computers: tested on IBM RS 6000, and PC                    *)
(*                                                              *)
(*  language: MATHEMATICA 2.2                                   *)
(*                                                              *)
(*  authors: W. Hereman and W. Zhuang                           *)
(*           Dept. of Mathematical and Computer Sciences        *)
(*           Colorado School of Mines,                          *)
(*           Golden, CO 80401-1887, USA                         *)
(*                                                              *)
(*  last updated: January 31, 2007    (23:50)                   *)
(*  previously updated: May 29, 1995  (21:30)                   *)
(*                                                              *)
(*  Copyright 1995                                              *)
(*                                                              *)
(****************************************************************)

BeginPackage["hirota`", "Global`"]

hirota::usage = " hirota[B,name,n,Rtest3soliton,Rtest4soliton,test3soliton,
                         test4soliton] gives a n-soliton solution of the
                         bilinear equation, after testing its existence.
                         B stands for the bilinear form B(f,g) for the given
                         PDE. The user can assign a name to the PDE, e.g. KdV.
                         The parameter n refers to the n-soliton solution one 
                         wants to calculate. The n is either 1, 2 or 3.
                         Random number tests as well as symbolic tests for
                         the existence of n-solitons can be performed."

Dx::usage = "Dx[n,x,f,g] is a n_th order Hirota bilinear operator in the
             variable x. Similar definitions for operators Dy and DT.
             Note that DT is used instead of Dt (which was defined in 
             Mathematica for total derivative)."

Dxt::usage = "Dxt[m,n,x,t,f,g] is a general form of Hirota's bilinear operator,
              m_th order in the variable x and n_th order in the variable t."

P::usage = "P[k,m,l] gives the polynamial corresponding to the bilinear 
            operators."

B::usage = "B[f,g] is a given bilinear form of the PDE.
            For instance, B[f_,g_]:=Dx[4,x,f,g]+Dxt[1,1,x,t,f,g] 
            for the KdV equation."

(* Setup for collecting the data in a log file, transcript of computations *)

logfile="";
OpenLog[filename_String]:=(logfile=OpenWrite[filename];
  If[ logfile===$Failed, Return[]];
  AppendTo[$Echo,logfile];AppendTo[$Output,logfile];);
CloseLog[]:=($Echo=Complement[$Echo,{logfile}];
  $Output=Complement[$Output,{logfile}];
  Close[logfile]);

(* Block 1: Tool to strip off numerical factors *)
stripper[expr_]:=Block[{tempexpr,ne,de,le,pp,i}, 
tempexpr=Factor[expr]; le=Length[tempexpr];
If[(le != 0 && Head[tempexpr]===Times),
   For[i=1,i<=le-1,i++,pp=Part[tempexpr,i];
       If[NumberQ[pp],tempexpr=tempexpr/pp]]];
Return[tempexpr]];
(* End of the Block 1: stripper[expr] to strip off numerical factors *)

(* Block 2: Friendly comment generator *)
COMMENTINTER[name_]:=Block[{},
Print["/*********************************************************/"];
Print["/*    WELCOME TO THE MATHEMATICA PROGRAM HIROTA.M        */"];
Print["/*            BY WILLY HEREMAN AND WUNING ZHUANG         */"];
Print["/*              FOR THE CALCULATION OF SOLITONS          */"];
Print["/*           OF THE ", name," EQUATION                   */"];
Print["/*              WITH HIROTA'S METHOD                     */"];
Print["/*     Version 1.0 firts released on May 29, 1995        */"];
Print["/*      Last updated on January 25, 2007                 */"];
Print["/*               Copyright 1995-2007                     */"];
Print["/*********************************************************/"]]
(* End of the Block 2: Commentinter[name] for the friendly interface *)

(* Block 3: Dispersion Relation *)
dispersion[P_]:=Block[{K,L,M,i},
  orgsol = M[i]/. Solve[P[K[i],-M[i], L[i]]==0, M[i]];  
  lenorgsol = Length[orgsol];
  OM[i_]:= M[i] /. Part[Solve[{P[K[i],-M[i],L[i]]==0,M[i]!=0},M[i]],1]];

(* End of the Block 3: dispersion[P] for the dispersion law *)

(* Block 4: Hirota Condition Tester for 3 and 4 solitons *)
condition[P_,N_,MYRANDOM_]:=Block[{ii,iii,j,ff,test,num},
If[ MYRANDOM,
   Print["Starting the random test(s) for the existence of a "];
   Print[ N," soliton solution."];
   NRandom[N],
   Print["Starting the symbolic test for the existence of a "];
   Print[ N," soliton solution."]];
term1[NN_]:=P[ Sum[K[I]*S[I],{I,1,NN}],-Sum[OM[I]*S[I],{I,1,NN}],
               Sum[L[I]*S[I],{I,1,NN}]];
term2[NN_]:=Product[Product[
            P[K[I]*S[I]-K[J]*S[J],-OM[I]*S[I]+OM[J]*S[J],L[I]*S[I]-L[J]*S[J]]
            *S[I]*S[J],{J,I+1,NN}],{I,1,NN}];
term[NN_]:=term1[NN]*term2[NN];
t2[I_,II_,F_]:=F /.{ S[I] -> II };
test=term[N];
Do[ test=t2[J,-1,test]+t2[J,1,test], {J,1,N}];
test=Expand[test];
test=Factor[test];
(* Put back in as before *)
Clear[K,L,OM];  dispersion[P]; 
If[ test =!= 0, 
   If[ MYRANDOM==True,  
       Print["The equation did not pass the random number test(s) for "];
       Print["the existence of a ", N," soliton solution."];
       If[ NumberQ[test] != True,
           Print["The condition ",test," = 0 must be satisfied."]],
        Print["The equation did not pass the symbolic test for "];
        Print["the existence of a ", N," soliton solution."];
        Print["The condition ",test," = 0 must be satisfied."] ];
   Return[False],
   If[ MYRANDOM==True,
      Print["The equation passed the random number test(s) for the existence"],
      Print["The equation passed the symbolic test for the existence"]];
   Print["of a ",N," soliton solution."]; 
   Return[True] ] ] (* End of If test =!= 0 and End of block *)
(* End of the Block 4: condition(P,N,RANDOM), for n soliton test with n>=3. *)

(* Block 5: Hirota Condition Tester for 1 and 2 solitons *)
condition1AND2soliton[P_]:=Block[{cond1,cond2,kk,mm,ll},
cond1=Simplify[P[0,0,0]];
cond2=Simplify[P[kk,mm,ll]-P[-kk,-mm,-ll]];
If [(cond1 =!= 0) || (cond2 =!= 0), 
   Print["There is no one- or two-soliton solution."];
   If [cond1 =!= 0, 
      Print["The bilinear operator contains a constant term."];
      Print["The condition ",Factor[cond1]," = 0 should be satisfied."] ];
   If [cond2 =!= 0,        
      Print["The bilinear operator is not even."] ];
   Return[False], 
   Print["The equation has at least a one- and two-soliton solution."];
   Return[True] ]  ]
(* End of the Block 5: condition_1_2soliton, test for one and two-solitons *)

(* Block 6: Coupling coefficients a_ij(P) *)
aij[P_]:=Block[{},
aa[ii_,jj_]:=Factor[-P[K[ii]-K[jj],-OM[ii]+OM[jj],L[ii]-L[jj]]/P[K[ii]+K[jj],
        -OM[ii]-OM[jj],L[ii]+L[jj]]];
Print["The coefficient a[I,J] is calculated via the polynomial form."];
Print["The polynomial is P[K,-OMEGA,L] = ", P[K,-OMEGA,L] ];
Print["The coefficient a[I,J] = ", Factor[aa[I,J]] ] ] 
(* End of the Block 6: a_ij(P), two-soliton coupling coefficients *)

(* Block 7: Coupling coefficients b_123[P] *)
b123[P_]:=Block[{}, bb[1,2,3]=aa[1,2]*aa[1,3]*aa[2,3];
Print["The coefficient b[1,2,3] is calculated via the polynomial form."];
Print["The coefficient b[1,2,3] = ", Factor[bb[1,2,3]] ] ]
(* End of the Block 7: b_123[P], three-soliton coupling coefficient *)

(* Block 8: Definition of Hirota's bilinear operators and polynomial P *)
DT[n_,t_,f_,g_]:=Sum[(-1)^(n-j)*n!*D[f,{t,j}]*D[g,{t,n-j}]
                     /(j!*(n-j)!),{j,0,n}];
Dx[n_,x_,f_,g_]:=Sum[(-1)^(n-j)*n!*D[f,{x,j}]*D[g,{x,n-j}]
                     /(j!*(n-j)!),{j,0,n}];
Dy[n_,y_,f_,g_]:=Sum[(-1)^(n-j)*n!*D[f,{y,j}]*D[g,{y,n-j}]
                     /(j!*(n-j)!),{j,0,n}];
Dxt[m_,n_,x_,t_,f_,g_]:=Sum[(-1)^(m-j)*m!*Sum[(-1)^(n-i)*n!*
                            D[f,{x,j},{t,i}]*D[g,{x,m-j},{t,n-i}]
                            /(i!*(n-i)!),{i,0,n}]/(j!*(m-j)!),{j,0,m}];
(* *)
(* old code 1995 works for v. 4.0 of Mathematica but not under v. 5.0 *)
(*
P[k_,m_,l_]:=Coefficient[B[1,Exp[l*y+k*x+m*t+cst]], Exp[l*y+k*x+m*t+cst]]; 
*)
(* *)
(* new code Hereman, January 31, 2007 works for Mathematica 4.0 and 5.0 *)
P[k_,m_,l_] := Coefficient[B[1,Exp[l*y+k*x+m*t+cst]] /. Exp[__] -> zz, zz, 1];
(* End of Block 8: Hirota bilinear operators and polynomial P *)

(* Block 9: Conditional random number generator *)
NRandom[n_]:=Block[{i,j,count,memb,membd,tmemb,d,iij},
memb={0,1}; membd={0,1}; tmemb=membd; i=1; count=0;
While[ i <= n, count=count+1; K[i]=Random[Integer, 20];
       If[ MemberQ[memb,K[i]] != True,
           j=i-1;
           While[ j>0, d=Abs[K[i]-K[j]];
                  If[ MemberQ[membd,d] != True, Append[membd,{d}]; j=j-1,
                      membd=tmemb;j=-1 ]];
           If[ j != -1, Append[memb,{K[i]}]; tmemb=membd; i=i+1]];
       If[ count>=10,
           K[1]=2+Random[Integer,{0,3}];  K[2]=6+Random[Integer,{0,3}];
           K[3]=11+Random[Integer,{0,3}]; K[4]=15+Random[Integer,{0,3}];
           i=5 ]];  (* End While i < n *)
  Print["Wavenumbers K[I] selected for the random number test(s): "]; 
  Do[Print["for this test K[",iij,"] = ", K[iij]], {iij,1,n}];
If[ lfree != True, i=1; count=0;
    While[ i<=n, count=count+1; L[i]=Random[Integer, 20];
           If[ MemberQ[memb,L[i] ] != True, j=i-1;
               While[ j>0, d=Abs[L[i]-L[j]];
                      If[ MemberQ[membd,d] != True, membd=Append[membd, {d} ];
                          j=j-1, membd=tmemb; j=-1 ]];
               If[ j != -1, memb=Append[memb, {L[i]}]; tmemb=membd; i=i+1] ];
           If[ count>=10,
           L[1]=2+Random[Integer,{0,3}]; L[2]=6+Random[Integer,{0,3}];
           L[3]=10+Random[Integer,{0,3}];L[4]=14+Random[Integer,{0,3}];
           i=5] ]; 
           Print["Wavenumbers L[I] selected for the random number test(s): "]; 
           Do[Print["for this test L[",iij,"] = ", L[iij]], {iij,1,n}] ]] 
(* End of the Block 9: NRandom[n], selection of random integers *)
(* for K[i], and also L[i] in the 2D case *)

(* Block 10: Hirota routine, main program *)
hirota[B_,name_,N_,Rtest3soliton_,Rtest4soliton_,test3soliton_,test4soliton_]:=
Block[{NN,nt,result2,result3,r3,r4,lfree}, NN = N;
COMMENTINTER[name]; 
Print["The equation in f corresponding to the given bilinear operator is "];
Print[stripper[Factor[B[f[x,y,z,t], f[x,y,z,t]]]], " = 0"]; 
P[k,m,l];
Print["For this equation the polynomial P(K,-OMEGA,L)= ", P[K,-OMEGA,L] ];
result2=condition1AND2soliton[P];
If[ result2==True,
   If[ N==1,
      Print["Starting the construction of the one-soliton solution."] ];
   dispersion[P]; Print["For the ", name," equation, "];
   If[ lenorgsol > 1, 
      Print["there are ", lenorgsol, " different dispersion relations;"];
      Print["the first one is: ", Part[orgsol,1]] ];
   Print["we use the dispersion relation: "]; Print[" OMEGA[I] = ", OM[I] ];
   lfree=FreeQ[ OM[II],L[II] ];
   If[ lfree,
     Print["In the Expansion of f we use THETA = K X - OMEGA T + CST."],
     Print["In the Expansion of f we use THETA = K X - OMEGA T + L Y + CST."]];
   If[ Rtest3soliton>0, nt=1;
      While[ nt<=Rtest3soliton,
            r3=condition[P,3,True];
            If[ r3==True, nt=nt+1, nt=Rtest3soliton+1] ], r3=True ];
   If[ test3soliton==True,
       If[ r3==True, result3=condition[P,3,False], result3=False ],
       If[ r3==True, result3=True, result3=False] ];
   If[ (Rtest4soliton>0) && (result3==True), nt=1;
      While[ nt<=Rtest4soliton, r4=condition[P,4,True];
             If[ r4==True, nt=nt+1, nt=Rtest4soliton+1 ] ],
      If[ result3==True, r4=True, r4=False ] ];
   If[ (test4soliton==True) && (result3==True),
      If[ r4==True, result4=condition[P,4,False] ] ];
   If[ NN==3,
      If[ result3==True,
          Print["Starting the construction of the three-soliton solution."];
          aij[P]; b123[P], NN=2 ] ];
   If[ NN==2,Print["Starting the construction of the two-soliton solution."];
       aij[P] ]; output[NN] ] ]  
(* End of the Block 10: *)
(* 
Hirota(B,name,N,Rtest_3soliton,Rtest_4soliton,test_3soliton,test_4soliton) 
for the Hirota method 
*)

(* Block 11: Hirota routine, main program *)
output[N_]:=Block[{},
  If[ N==1, f=1+Exp[THETA[1]] ];
  If[ N==2, f=1+Exp[THETA[1]]+Exp[THETA[2]]+a[1,2]*Exp[THETA[1]+THETA[2]] ];
  If[ N==3,  
      f=1+Exp[THETA[1]]+Exp[THETA[2]]+Exp[THETA[3]]
        +a[1,2]*Exp[THETA[1]+THETA[2]]+a[1,3]*Exp[THETA[1]+THETA[3]]
	+a[2,3]*Exp[THETA[2]+THETA[3]]
	+b[1,2,3]*Exp[THETA[1]+THETA[2]+THETA[3]] ];
  Print["The function f = ", f];
  F = f;
 Print["At the end of the computations the form of the function f"];
 If[ FreeQ[ OM[I], L[I] ], 
             THETA1[i_]:=K[i]*X-OM[i]*T+CST[i];
             Do[ THETA[j]=THETA1[j], {j,1,N}]; THETA[i]=THETA1[i],
             THETA2[i_]:=K[i]*X-OM[i]*T+L[i]*Y+CST[i]; THETA[i]=THETA2[i];
             Do[ THETA[j]=THETA2[j], {j,1,N}] ]; 
 If[ N==1, Print["is available."]];
 If[ N==2, a[1,2]=aa[1,2];
           Print["and the coefficient a[1,2] are explicitly available."]];
 If[ N==3, a[1,2]=aa[1,2];a[1,3]=aa[1,3];a[2,3]=aa[2,3];b[1,2,3]=bb[1,2,3];
           OMEGA[3]=OM[3];
 Print["and the coefficients a[i,j] and b[1,2,3] are explicitly available."]];
 Print["The explicit forms of OMEGA[i] and THETA[i] are also available."];
   Print["The Explicit form of f can be obtained by typing EXPRF."];
   Do[ OMEGA[j]=OM[j], {j,1,N}]; OMEGA[i]=OM[i];   EXPRF=f ]
(* End of the Block 11: output(n), returning f and the coupling coefficients *)

EndPackage[];
Print["Package hirota.m was successfully loaded."];

(* ****************** End of the program HIROTA.M *************** *)

