© Copyright Notice

Copyright © by Unal Goktas ((Department of Computer Science and Engineering, Texas A&M University, College Station, Texas, USA) and Willy Hereman (Department of Mathematical and Computer Sciences, Colorado School of Mines, Golden, Colorado, USA). No part of the PDESolitonSolutions package can be sold or reproduced without written consent of the authors. Wolfram Research, Inc. (Champaign, Illinois, USA) is the holder of the copyright to the Mathematica software system.

About the Package: PDESolitonSolutions

Nonlinear partial differential equations (PDEs) possess a number of remarkable properties reflecting their rich mathematical structure. Such properties include the Painleve property, nontrivial prolongation and bi-Hamiltonian structures, soliton solutions, Lax pairs, Backlund transformations, symmetries and conservation laws. These properties reveal the adequacy of PDEs as models for physical relevant phenomena. The search for these intrinsic properties becomes more attractive due to the availability of CAS such as Mathematica. Indeed, Mathematica is well suited and often necessary tool to perform the computations inherent in the investigation of integrability.

PDESolitonSolutions is a collection of Mathematica functions addressing computations of soliton solutions in analyzing nonlinear PDEs.  A description of the simplified Hirota method and code is given in the accompanying paper:  Willy Hereman and Unal Goktas, Symbolic computation of solitary wave solutions and solitons through homogenization of degree, in `Nonlinear and Modern Mathematical Physics,’ Proceedings of the 6th International Workshop on Nonlinear and Modern Mathematical Physics (NMMP-2022), Eds.: S. Manukure and W.-X. Ma,     Springer Verlag, New York (2024), 61pp.  This book chapter is also available on the ArXiv.

A complete set of mathematical definitions and the concepts used in the package and this manual are beyond the scope of this document. For the formal definitions of the concepts and the theory behind the functions in the PDESolitonSolutions package, consult the references given in the References section (below). In this manual, a basic knowledge of Mathematica is assumed. For further information on Mathematica, consult the standard Mathematica references by Wolfram Research, Inc.

Loading the Package

To access all of the functions in PDESolitonSolutions, all you need to do is to load the file PDESolitonSolutions.m by typing the following line in a Mathematica notebook:

Set the directory so that Mathematica knows where to get  the code. For example:

(* SetDirectory["Dropbox\\Homogenization-Paper-2022\\Homogenization-code-UG"] *)

SetDirectory["C:\\data\\talks\\talks22\\Homogenization-code-UG"]

Load the package. The command below will generate two lines of messages.

Get["PDESolitonSolutions-March24-2024.m"]

Usage

??PDESolitonSolutions

??PDESolitonSolutionsTimeConstraint

Examples

If you haven't already loaded the package, load it now.

The Burgers Equation

burgerseq:=(∂_t)u[x,t]+2 u[x,t] (∂_x)u[x,t]-(∂_({x,2}))u[x,t]==0

burgerseq

Compute 1-soliton solutions:

burgers1sol = PDESolitonSolutions[burgerseq,u,{x},t,{1},PrintInformation->True]

burgers1sol = PDESolitonSolutions[burgerseq,u,{x},t,{1},PrintInformation->False]

(* Testing the solution symbolically *)

burgerseq/.DSolve`DSolveToPureFunction[burgers1sol]//Simplify

PDESolitonSolutions[burgerseq,u,{x},t,{1},UseTruncatedPainleveExpansion->False,PrintInformation->True]

PDESolitonSolutions[burgerseq,u,{x},t,{1},UseTruncatedPainleveExpansion->False,PrintInformation->False]

Compute 1 and 2-soliton solutions:

burgers1and2sol= PDESolitonSolutions[burgerseq,u,{x},t,2, PrintInformation -> True]

(* Testing both solutions symbolically *)

burgerseq/.DSolve`DSolveToPureFunction[burgers1and2sol]//Simplify

Compute only the 2-soliton solution:

burgerseqsoliton2=PDESolitonSolutions[burgerseq,u,{x},t,{2}]

burgerseqsoliton2evaluated=burgerseqsoliton2[[1,1,2]]/.{k[1]->-2,k[2]->3,delta[1]->2,delta[2]->-3}

With[{sol=burgerseqsoliton2evaluated},Manipulate[Plot[sol,{x,-25,25},PlotRange->All],{t,-10,2,0.1},SaveDefinitions->True,ControlPlacement->{Top}]]

Compute 1, 2 and 3-soliton solutions:

burgers3sol =PDESolitonSolutions[burgerseq,u,{x},t,3, PrintInformation -> True]

Compute only the 3-soliton solution:

burgers3solagain = PDESolitonSolutions[burgerseq,u,{x},t,{3}]

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(burgerseq[[1]]/.DSolve`DSolveToPureFunction[burgers3solagain [[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{10}]

The Korteweg-de Vries (KdV) Equation

kdveq:=(∂_t)u[x,t]+6 u[x,t] (∂_x)u[x,t]+(∂_({x,3}))u[x,t]==0

kdveq

Compute 1-soliton solutions:

kdvsoliton1 = PDESolitonSolutions[kdveq,u,{x},t,{1},PrintInformation->True]

(* Testing the solution symbolically *)

kdveq/.DSolve`DSolveToPureFunction[kdvsoliton1]//Simplify

Compute 1 and 2-soliton solutions:

kdvsoliton1and2 = PDESolitonSolutions[kdveq,u,{x},t,2]

(* Testing the solution symbolically *)

kdveq/.DSolve`DSolveToPureFunction[kdvsoliton1and2]//Simplify

Compute only the 2-soliton solution:

kdvsoliton2=PDESolitonSolutions[kdveq,u,{x},t,{2},PrintInformation->True]

(* Testing the solution numerically ten times *)

Table[PossibleZeroQ[Chop[Quiet[(kdveq[[1]]/.DSolve`DSolveToPureFunction[kdvsoliton2[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{10}]

(* Testing the solution symbolically *)

kdveq/.DSolve`DSolveToPureFunction[kdvsoliton2]//Simplify

kdvsoliton2evaluated=kdvsoliton2[[1,1,2]]/.{k[1]->1,k[2]->-2,delta[1]->2,delta[2]->-3}

(* To animate the collision of soliton click on the "+" sign on the time indicator bar, and the "run" arrow *)

With[{sol=kdvsoliton2evaluated},Manipulate[Plot[sol,{x,-30,30},PlotRange->All],{t,-5,5,0.01},SaveDefinitions->True,ControlPlacement->{Top}]
]

Compute only the 3-soliton solution:

kdvsoliton3=PDESolitonSolutions[kdveq,u,{x},t,{3}, PrintInformation -> True]

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(kdveq[[1]]/.DSolve`DSolveToPureFunction[kdvsoliton3[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{10}]

kdvsoliton3evaluated=kdvsoliton3[[1,1,2]]/.{k[1]->2,k[2]->3/2,k[3]->1,delta[1]->1,delta[2]->-2,delta[3]->3}

With[{sol=kdvsoliton3evaluated},
Manipulate[Plot[sol,{x,-25,25},PlotRange->All],{t,-10,10,0.01},SaveDefinitions->True,ControlPlacement->{Top}]
]

kdvsolitons1and2=PDESolitonSolutions[kdveq,u,{x},t,2]

kdvsolitons1through3=PDESolitonSolutions[kdveq,u,{x},t,3]

Test the first solution 10 times:

Table[PossibleZeroQ[Chop[Quiet[(kdveq[[1]]/.DSolve`DSolveToPureFunction[kdvsolitons1through3[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{10}]

Test the second solution a dozen times:

Table[PossibleZeroQ[Chop[Quiet[(kdveq[[1]]/.DSolve`DSolveToPureFunction[kdvsolitons1through3[[2]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{12}]

Test the third solution fifteen times:

Table[PossibleZeroQ[Chop[Quiet[(kdveq[[1]]/.DSolve`DSolveToPureFunction[kdvsolitons1through3[[3]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{15}]

(* The computation would be aborted if, e.g., PDESolitonSolutionsTimeConstraint -> 1/10, because it could not be finished in 1/10 of a second! *)

kdvsolitons1through3=PDESolitonSolutions[kdveq,u,{x},t,3, PrintInformation -> True, PDESolitonSolutionsTimeConstraint -> 1/10]

The Lax Equation

Laxfifthorderpde:=(∂_t)u[x,t]+30 ((u[x,t])^2) (∂_x)u[x,t]+20 (∂_x)u[x,t] (∂_({x,2}))u[x,t]+10 u[x,t] (∂_({x,3}))u[x,t]+(∂_({x,5}))u[x,t]==0

Laxfifthorderpde

Compute 1-soliton solutions:

laxonesoliton = PDESolitonSolutions[Laxfifthorderpde,u,{x},t,{1},PrintInformation->True]

(* Testing the solution symbolically *)

Laxfifthorderpde/.DSolve`DSolveToPureFunction[laxonesoliton]//Simplify

Compute only the 2-soliton solution:

laxtwosoliton= PDESolitonSolutions[Laxfifthorderpde,u,{x},t,{2}, PrintInformation -> True]

laxtwosoliton= PDESolitonSolutions[Laxfifthorderpde,u,{x},t,{2}, UseTruncatedPainleveExpansion->False, PrintInformation -> True]

(* Testing the solution numerically a dozen times *)

Table[PossibleZeroQ[Chop[Quiet[(Laxfifthorderpde[[1]]/.DSolve`DSolveToPureFunction[laxtwosoliton[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{12}]

PDESolitonSolutions[Laxfifthorderpde,u,{x},t,{2}]//FullSimplify

Compute only the 3-soliton solution:

Laxfifthorderpde3soliton=PDESolitonSolutions[Laxfifthorderpde,u,{x},t,{3}, PrintInformation -> True]

(* Testing the solution numerically a dozen times *)

Table[PossibleZeroQ[Chop[Quiet[(Laxfifthorderpde[[1]]/.DSolve`DSolveToPureFunction[Laxfifthorderpde3soliton[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{12}]

Laxfifthorderpdesoliton1through3=PDESolitonSolutions[Laxfifthorderpde,u,{x},t,3]

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(Laxfifthorderpde[[1]]/.DSolve`DSolveToPureFunction[Laxfifthorderpdesoliton1through3])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{10}]

The Sawada-Kotera Equation

SKfifthorderpde:=(∂_t)u[x,t]+5 ((u[x,t])^2) (∂_x)u[x,t]+5 (∂_x)u[x,t] (∂_({x,2}))u[x,t]+5 u[x,t] (∂_({x,3}))u[x,t]+(∂_({x,5}))u[x,t]==0

SKfifthorderpde

Compute 1-soliton solutions:

SKonesoliton= PDESolitonSolutions[SKfifthorderpde,u,{x},t,{1},PrintInformation->True]//Simplify

(* Testing the solution symbolically *)

SKfifthorderpde/.DSolve`DSolveToPureFunction[SKonesoliton]//Simplify

Compute only the 2-soliton solution:

PDESolitonSolutions[SKfifthorderpde,u,{x},t,{2}]//Simplify

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(SKfifthorderpde[[1]]/.DSolve`DSolveToPureFunction[SKonesoliton[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{10}]

SKtwosoliton = PDESolitonSolutions[SKfifthorderpde,u,{x},t,{2},UseTruncatedPainleveExpansion->False, PrintInformation -> True]

Compute only the 3-soliton solution:

SKfifthorderpde3soliton=PDESolitonSolutions[SKfifthorderpde,u,{x},t,{3}]

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(SKfifthorderpde[[1]]/.DSolve`DSolveToPureFunction[SKfifthorderpde3soliton[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{10}]

SKsolitons1through3 = PDESolitonSolutions[SKfifthorderpde,u,{x},t,3(* ,UseTruncatedPainleveExpansion->False, *) (*, PrintInformation -> True*)]

The Kaup-Kupershmidt Equation

KKfifthorderpde:=(∂_t)u[x,t]+20 ((u[x,t])^2) (∂_x)u[x,t]+25 (∂_x)u[x,t] (∂_({x,2}))u[x,t]+10 u[x,t] (∂_({x,3}))u[x,t]+(∂_({x,5}))u[x,t]==0

KKfifthorderpde

Compute 1-soliton solutions:

KKonesoliton = PDESolitonSolutions[KKfifthorderpde,u,{x},t,{1},PrintInformation->True](* //Simplify *)

(* Testing the solution symbolically *)

KKfifthorderpde/.DSolve`DSolveToPureFunction[KKonesoliton]//Simplify

Compute only the 2-soliton solution:

KKtwosoliton = PDESolitonSolutions[KKfifthorderpde,u,{x},t,{2}, PrintInformation -> True]

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(KKfifthorderpde[[1]]/.DSolve`DSolveToPureFunction[KKtwosoliton[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{10}]

Compute only the 3-soliton solution (this computation is very long and will take some time, that is, several minutes on laptops). The computation can be finished under the default setting PDESolitonSolutionsTimeConstraint -> 900 (15 minutes).

KKfifthorderpde3soliton=PDESolitonSolutions[KKfifthorderpde,u,{x},t,{3}]

(* Testing the solution numerically 3 times. This will also take several minutes. Note that 10 tests would take about 15 minutes. *)

Table[PossibleZeroQ[Chop[Quiet[(KKfifthorderpde[[1]]/.DSolve`DSolveToPureFunction[KKfifthorderpde3soliton[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{3}]

{True,True,True}

Compute only the 3-soliton solution again showing the details of the computation (this computation is very long and will take some time, that is, several minutes on laptops):

KKfifthorderpde3soliton=PDESolitonSolutions[KKfifthorderpde,u,{x},t,{3}, PrintInformation -> True]

(* Note that the computation could not be finished if, for example. PDESolitonSolutionsTimeConstraint -> 90 (about 3 minutes) *)

KKfifthorderpde3soliton=PDESolitonSolutions[KKfifthorderpde,u,{x},t,{3}, PrintInformation->True, PDESolitonSolutionsTimeConstraint -> 90]

The Modified KdV Equation

mkdveq:=(∂_t)u[x,t]+24 ((u[x,t])^2) (∂_x)u[x,t]+(∂_({x,3}))u[x,t]==0

mkdveq

Try to compute 1-soliton solutions (you will notice that the current version of the software fails to do this):

PDESolitonSolutions[mkdveq,u,{x},t,{1},PrintInformation->True]

The Fisher Equation with convection

fisherwithconvection:=(∂_t)u[x,t]+α u[x,t] (∂_x)u[x,t]-(∂_({x,2}))u[x,t]-u[x,t] (1-u[x,t])==0

fisherwithconvection

Compute 1-soliton solutions:

fisherconv1sol = PDESolitonSolutions[fisherwithconvection,u,{x},t,{1},PrintInformation->True]//Simplify

(* testing the solution symbolically *)

(fisherwithconvection/.DSolve`DSolveToPureFunction[fisherconv1sol[[1]]])//Simplify

FullSimplify[fisherconv1sol//TrigToExp]

Compute only the 2-soliton solution:

fisherconv2sol =PDESolitonSolutions[fisherwithconvection,u,{x},t,{2}, PrintInformation-> True]

So the Fisher equation with convection has only a 1-soliton solution.

The Fisher Equation

fisher:=(∂_t)u[x,t]-(∂_({x,2}))u[x,t]-u[x,t] (1-u[x,t])==0

fisher

Compute 1-soliton solutions:

fisher1sol = PDESolitonSolutions[fisher,u,{x},t,{1},PrintInformation->True] (* //Simplify *)

(* Testing both solutions symbolically *)

(fisher/.DSolve`DSolveToPureFunction[fisher1sol])//Simplify

FullSimplify[fisher1sol//TrigToExp]

Try  to compute a 2-soliton solution:

PDESolitonSolutions[fisher,u,{x},t,{2}, PrintInformation -> True]

So the Fisher equation has only a 1-soliton solution.

The FitzHugh-Nagumo equation with convection

Replace α by ((mm-2)/(Sqrt[mm])) in the FitzHugh-Nagumo equation:

fitzhughnagumowithm:=(∂_t)u[x,t]+α u[x,t] (∂_x)u[x,t]-(∂_({x,2}))u[x,t]+u[x,t] (1-u[x,t]) (a-u[x,t])==0/. {α->((mm-2)/(Sqrt[mm]))}

fitzhughnagumowithm

Compute 1-soliton solutions:

fitzhughnagumowithm1sol = PDESolitonSolutions[fitzhughnagumowithm,u,{x},t,{1}, PrintInformation -> True] (* //FullSimplify *)

(* Testing the solutions symbolically *)

(fitzhughnagumowithm/.DSolve`DSolveToPureFunction[fitzhughnagumowithm1sol])//Simplify

Compute 2-soliton solutions (this is quite cumbersome and will take a long time):

fitzhughnagumowithm2sol = PDESolitonSolutions[fitzhughnagumowithm,u,{x},t,{2}, PrintInformation -> True] (* //FullSimplify *)

Note that the solutions returned above are bi-soliton solutions. Test each of these solutions numerically 10 times (no worries this is fast):

Table[PossibleZeroQ[Chop[Quiet[(fitzhughnagumowithm[[1]]/.DSolve`DSolveToPureFunction[fitzhughnagumowithm2sol])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{10}]

Return back to the original FitzHugh-Nagumo equation:

fitzhughnagumo:=(∂_t)u[x,t]+α u[x,t] (∂_x)u[x,t]-(∂_({x,2}))u[x,t]+u[x,t] (1-u[x,t]) (a-u[x,t])==0

fitzhughnagumo

Compute 1-soliton solutions:

fitzhughnagumo1sol = PDESolitonSolutions[fitzhughnagumo,u,{x},t,{1},PrintInformation->True]

(* Testing the solutions symbolically *)

(fitzhughnagumo/.DSolve`DSolveToPureFunction[fitzhughnagumo1sol])//Simplify

A Burgers equation with a cubic source term

burgerswithcubicwithm:=(∂_t)u[x,t]+α u[x,t] (∂_x)u[x,t]-(∂_({x,2}))u[x,t]-3 u[x,t] (2-u[x,t]) (u[x,t]+1)==0 /. α->((3 mm-2)/(Sqrt[mm]))

burgerswithcubicwithm

Compute 1-soliton solutions:

burgerscubicwithm1sol = PDESolitonSolutions[burgerswithcubicwithm,u,{x},t,{1},PrintInformation->True]

(* Testing the solutions symbolically *)

(burgerswithcubicwithm/.DSolve`DSolveToPureFunction[burgerscubicwithm1sol ])//Simplify

burgerscubicwithm1sollogderivtf = PDESolitonSolutions[burgerswithcubicwithm,u,{x},t,{1},UseTruncatedPainleveExpansion->False, PrintInformation->True]

(* Testing the solutions symbolically *)

(burgerswithcubicwithm/.DSolve`DSolveToPureFunction[burgerscubicwithm1sollogderivtf])//Simplify

Compute 2-soliton solutions (computing these bi-soliton solutions will take several minutes):

burgerscubicwithm2sol = PDESolitonSolutions[burgerswithcubicwithm,u,{x},t,{2},PrintInformation->True]

(* Testing the bi-soliton solutions numerically 10 times (which will be fast) *)

Table[PossibleZeroQ[Chop[Quiet[(burgerswithcubicwithm[[1]]/.DSolve`DSolveToPureFunction[burgerscubicwithm2sol])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{10}]

Return back to the original equation:

burgerswithcubicorg:=(∂_t)u[x,t]+α u[x,t] (∂_x)u[x,t]-(∂_({x,2}))u[x,t]-3 u[x,t] (2-u[x,t]) (u[x,t]+1)==0 (* /.  α->((3 mmm-2)/(Sqrt[mmm])) *)

burgerswithcubicorg

Compute 1-soliton solutions:

burgerscubicorg1sol = PDESolitonSolutions[burgerswithcubicorg,u,{x},t,{1},PrintInformation->True]

(* Testing these solutions symbolically *)

(burgerswithcubicorg/.DSolve`DSolveToPureFunction[burgerscubicorg1sol ])//Simplify

Compute 2-soliton solutions (again this computation will take a long time):

burgerscubicorg2sol=PDESolitonSolutions[burgerswithcubicorg,u,{x},t,{2}, PrintInformation -> True]

The solutions returned above are bi-soliton solutions. Now, testing these solutions numerically 10 times (this is quite fast).

Table[PossibleZeroQ[Chop[Quiet[(burgerswithcubicorg[[1]]/.DSolve`DSolveToPureFunction[burgerscubicorg2sol])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{10}]

A wave equation with cubic source term

waveeqwithcubicterm:=(1/8) (∂_({t,2}))u[x,t]+(∂_t)u[x,t]+u[x,t] (∂_x)u[x,t]-(1/8) (∂_({x,2}))u[x,t]-u[x,t] (u[x,t]-1) (u[x,t]+2)==0

waveeqwithcubicterm

Compute 1-soliton solutions:

waveeqwithcubicterm1sol = PDESolitonSolutions[waveeqwithcubicterm,u,{x},t,{1},PrintInformation->True]

(* Testing the solution symbolically *)

(waveeqwithcubicterm/.DSolve`DSolveToPureFunction[waveeqwithcubicterm1sol])//Simplify

FullSimplify[waveeqwithcubicterm1sol]

Compute only the 2-soliton solution:

waveeqwithcubicterm2sol= PDESolitonSolutions[waveeqwithcubicterm,u,{x},t,{2}, PrintInformation -> True]//FullSimplify

The solution returned above is a bi-soliton solution.

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(waveeqwithcubicterm[[1]]/.DSolve`DSolveToPureFunction[waveeqwithcubicterm2sol[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{10}]

Compute only the 3-soliton solution (this is more cumbersome and slower, especially when fully simplified):

waveeqwithcubicterm3sol= PDESolitonSolutions[waveeqwithcubicterm,u,{x},t,{3}, PrintInformation -> True]//FullSimplify

The solution returned above is a tri-soliton solution.

(* Testing the solution symbolically *)

(waveeqwithcubicterm/.DSolve`DSolveToPureFunction[waveeqwithcubicterm3sol])//Simplify

A combined KdV-Burgers equation

combinedkdvburgers:=(∂_t)u[x,t]+6 u[x,t] (∂_x)u[x,t]+(∂_({x,3}))u[x,t]-5 β (∂_({x,2}))u[x,t]==0

combinedkdvburgers

Compute 1-soliton solutions:

combinedkdvburgers1sol = PDESolitonSolutions[combinedkdvburgers,u,{x},t,{1},PrintInformation->True] (* //FullSimplify *)

The first solution above is from the branch k[1] = β, the second from the branch k[1] = - β  Testing these solutions symbolically:

(combinedkdvburgers/.DSolve`DSolveToPureFunction[combinedkdvburgers1sol])//FullSimplify

An equation due to Calogero

calogeroeq:=(∂_t)u[x,t]-3 (3 u[x,t] ((((∂_x)u[x,t]))^2)+((u[x,t])^4) (∂_x)u[x,t]+((u[x,t])^2) (∂_({x,2}))u[x,t])-(∂_({x,3}))u[x,t]==0

calogeroeq

Compute 1-soliton solutions:

calogeroeq1sol = PDESolitonSolutions[calogeroeq,u,{x},t,{1}, PrintInformation-> True] (* // FullSimplify *)

transformedcalogeroeq=Numerator[Together[calogeroeq[[1]]/.u->Function[{x,t},v[x,t]^(1/2)]]]==0

Compute 1-soliton solutions for the transformed Calogero equation:

transformedcalogeroeq1sol = PDESolitonSolutions[transformedcalogeroeq,v,{x},t,{1},PrintInformation->True] (* //FullSimplify *)

(* Testing this solution symbolically *)

(transformedcalogeroeq/.DSolve`DSolveToPureFunction[transformedcalogeroeq1sol[[1]]])//Simplify

transformedcalogeroeq2sol = PDESolitonSolutions[transformedcalogeroeq,v,{x},t,{2}, PrintInformation -> True]

(* No two-soliton solution could be computed. *)

An equation with two but not three solitons

eqtwonothreesolitons:=(∂_t)u[x,t]+(15/784)*((u[x,t])^3)*(∂_x)u[x,t]+(15/28)*u[x,t]*(∂_x)u[x,t]*(∂_({x,2}))u[x,t]+(15/56)*u[x,t]^2*(∂_({x,3}))u[x,t]+(5/2)*(∂_({x,2}))u[x,t]*(∂_({x,3}))u[x,t]+(∂_x)u[x,t]*(∂_({x,4}))u[x,t]+u[x,t]*(∂_({x,5}))u[x,t]+(∂_({x,7}))u[x,t]==0

eqtwonothreesolitons

eqtwonothreesolitons1sol =PDESolitonSolutions[eqtwonothreesolitons,u,{x},t,{1},PrintInformation->True]

(* Testing the 1-soliton solution symbolically *)

eqtwonothreesolitons/.DSolve`DSolveToPureFunction[eqtwonothreesolitons1sol]//Simplify

(* Computing the 2-soliton solution. This is still very fast. *)

eqtwonothreesolitons2sol= PDESolitonSolutions[eqtwonothreesolitons,u,{x},t,{2},PrintInformation->True]

(* Testing the 2-soliton solution five times. *)

Table[PossibleZeroQ[Chop[Quiet[(eqtwonothreesolitons[[1]]/.DSolve`DSolveToPureFunction[eqtwonothreesolitons2sol[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],k[3],delta[1],delta[2],delta[3],x,t})]]],{5}]

(* Trying to compute a 3-soliton solution. This is cumbersome and takes a lot of time because it does not exist! *)

(* The time constraint setting $PDESolitonSolutionsTimeConstraint -> 900 inside the code will prevent that the execution would have to be aborted by the user. *)

eqtwonothreesolitons3sol= PDESolitonSolutions[eqtwonothreesolitons,u,{x},t,{3},PrintInformation->True]

PDESolitonSolutionsTimeConstraint has played a significant role in the above output. To make the code run longer than the default 900 seconds (15 min) but it might not help because the perturbation scheme does not terminate since a three soliton solution does not exist. Running the code with less time allocated will make it abort quicker but some of the output might be useful (in some cases) For example, the equation in homogenized form will be available, etc.

eqtwonothreesolitons3sol= PDESolitonSolutions[eqtwonothreesolitons,u,{x},t,{3},PrintInformation->True, PDESolitonSolutionsTimeConstraint->100]

Kadomtsev-Petviashvili (KP) equation (KP-II case)

Equation (1):

kppde:=(∂_x)((∂_t)u[x,y,t]+6 u[x,y,t] (∂_x)u[x,y,t]+(∂_({x,3}))u[x,y,t])+3(∂_({y,2}))u[x,y,t]==0

kppde

Compute 1-soliton solutions:

kppde1sol = PDESolitonSolutions[kppde,u,{x,y},t,{1},PrintInformation->True,UseTruncatedPainleveExpansion->False] (* //FullSimplify *)

(* Testing the solution symbolically *)

kppde/.DSolve`DSolveToPureFunction[kppde1sol]//Simplify

Compute 2-soliton solutions:

kppdesoliton2=PDESolitonSolutions[kppde,u,{x,y},t,{2}, PrintInformation -> True](*//Simplify*)

(* Testing the solution numerically 10 times *)

Table[PossibleZeroQ[Chop[Quiet[(kppde[[1]]/.DSolve`DSolveToPureFunction[kppdesoliton2[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->400]&)/@{k[1],k[2],l[1],l[2],delta[1],delta[2],x,y,t})],10^-6]],{10}]

kppdesoliton2evaluated=kppdesoliton2[[1,1,2]]/.{k[1]->1,k[2]->-3/2,l[1]->2,l[2]->7/2,delta[1]->1,delta[2]->3/2}

With[{sol=kppdesoliton2evaluated},
Manipulate[Plot3D[sol,{x,-30,10},{y,-30,10},PlotPoints -> 200,PlotRange->All],{t,-1.5,1.5,0.0001},SaveDefinitions->True,ControlPlacement->{Top}]
]

Compute 3-soliton solutions (output suppressed here because it is very long):

kppdesoliton3=PDESolitonSolutions[kppde,u,{x,y},t,{3}];

(* Testing the solution numerically 3 times. Note the 10^(-8) for chopping, any value smaller than 10^(-8) will be replaced by zero. The testing will take several minutes. *)

Table[PossibleZeroQ[Chop[Quiet[(kppde[[1]]/.DSolve`DSolveToPureFunction[kppdesoliton3[[1]]])/.((#->RandomReal[{-3,3},WorkingPrecision->300]&)/@{k[1],k[2],k[3],l[1],l[2],l[3],delta[1],delta[2],delta[3],x,y,t})],10^(-8)]],{3}]

kppdesoliton3evaluated=kppdesoliton3[[1,1,2]]/.{k[1]->1.1,k[2]->1.7,k[3]->2.5,l[1]->1.2,l[2]->1.8,l[3]-> 2.3,delta[1]->0,delta[2]->0,delta[3]->0};

(* Showing the 3-soliton collision. If needed, click on the "+" on the t-bar, then click on play to run the demo. *)

With[{sol=kppdesoliton3evaluated},
Manipulate[Plot3D[sol,{x,-15,15},{y,-15,15},PlotPoints -> 200,PlotRange->All],{t,-2.3,2.5,0.1},SaveDefinitions->True,ControlPlacement->{Top}]
]

A (3+1)-dimensional evolution equation (due to X. Geng and Y. Ma)

Reference: Geng, X., Ma, Y.: N-soliton solution and itsWronskian form of a (3+1)-dimensional nonlinear evolution equation. Phys. Lett. A 369(4), 285–289 (2007).

gengpdeinu:=3 (∂_({x,2},z))u[x,y,z,t]-(∂_y)(2 (∂_(x,t))u[x,y,z,t]+(∂_({x,4}))u[x,y,z,t]-2  (∂_x)u[x,y,z,t] (∂_({x,2}))u[x,y,z,t])+2 (∂_({x,1}))((∂_({x,2}))u[x,y,z,t] (∂_y)u[x,y,z,t])==0

gengpdeinu

Compute a 1-soliton solution:

gengpdeinu1sol = PDESolitonSolutions[gengpdeinu,u,{x,y,z},t,{1},PrintInformation->True] (* //FullSimplify *)

(* Testing the 1-soliton solution symbolically *)

gengpdeinu/.DSolve`DSolveToPureFunction[gengpdeinu1sol]//Simplify

Compute 2-soliton solutions:

gengpdeinu2sol=PDESolitonSolutions[gengpdeinu,u,{x,y,z},t,{2},PrintInformation->True]

(* Testing the solution numerically 10 times. Takes only a bit of time. *)

Table[PossibleZeroQ[Chop[Quiet[(gengpdeinu[[1]]/.DSolve`DSolveToPureFunction[gengpdeinu2sol[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->300]&)/@{k[1],k[2],l[1],l[2],m[1],m[2],delta[1],delta[2],x,t})]]],{10}]

(* The computation below is very cumbersome and does not finish in the allocated default time of 15 minutes. Run the computation first for 90 seconds of CPU time to see some of the steps *)

gengpdeinu3sol = PDESolitonSolutions[gengpdeinu,u,{x,y,z},t,{3},PrintInformation->True, PDESolitonSolutionsTimeConstraint -> 90]

PDESolitionSolutionsTimeConstraint has played a significant role in the above output. Even with the default setting (PDESolitonSolutionsTimeConstraint --> 900 ) the computation was aborted because it could not be finished within 15 minutes of CPU time. Increasing the time by setting PDESolitionSolutionsTimeConstraint --> 1200 does not appear to help. The software can not find a 3-soliton solution which might not exist for this example. This needs further investigation.

Another equation with two but not three solitons (due to Il’in et al.)

Reference: Il’in, I.A., Noshchenko, D.S., Perezhogin, A.S.: On classification of higher-order integrable nonlinear partial differential equations. Chaos, Solitons & Fractals 76, 278-281 (2015). doi: 10.1016/j.chaos.2015.04.004.

eqtwonothreesolitons:=(∂_t)u[x,t]+(15/784)*((u[x,t])^3)*(∂_x)u[x,t]+(15/28)*u[x,t]*(∂_x)u[x,t]*(∂_({x,2}))u[x,t]+(15/56)*u[x,t]^2*(∂_({x,3}))u[x,t]+(5/2)*(∂_({x,2}))u[x,t]*(∂_({x,3}))u[x,t]+(∂_x)u[x,t]*(∂_({x,4}))u[x,t]+u[x,t]*(∂_({x,5}))u[x,t]+(∂_({x,7}))u[x,t]==0

eqtwonothreesolitons

Compute a 1-soliton solution:

eqtwonothreesolitons1sol =PDESolitonSolutions[eqtwonothreesolitons,u,{x},t,{1},PrintInformation->True]  (* //FullSimplify *)

(* Testing the solution symbolically *)

eqtwonothreesolitons/.DSolve`DSolveToPureFunction[eqtwonothreesolitons1sol]//Simplify

Compute a 2-soliton solution:

eqtwonothreesolitons2sol= PDESolitonSolutions[eqtwonothreesolitons,u,{x},t,{2},PrintInformation->True]

(* Testing the solution numerically 5 times. This is fast. *)

Table[PossibleZeroQ[Chop[Quiet[(eqtwonothreesolitons[[1]]/.DSolve`DSolveToPureFunction[eqtwonothreesolitons2sol[[1]]])/.((#->RandomReal[{-2,2},WorkingPrecision->200]&)/@{k[1],k[2],delta[1],delta[2],x,t})]]],{5}]

Compute 3-soliton solution (which in this example does not exist):

eqtwonothreesolitons3sol= PDESolitonSolutions[eqtwonothreesolitons,u,{x},t,{3},PrintInformation->True, PDESolitonSolutionsTimeConstraint -> 300]

A three-soliton solution can not be computed (actually, it does not exist as explained in our paper). PDESolitionSolutionsTimeConstraint has played a significant role in the above output. Even with PDESolitionSolutionsTimeConstraint --> 900 or higher the code aborts.

Scope and Limitations

PDESolitonSolutions can handle equations that are polynomial in the (one) dependent variable. Up to three space variables and one time variable are allowed. Currently no systems of PDEs can be handled. The code can not handle equations such as the mKdV, sine-Gordon, nonlinear Schrodinger (NLS), Benjamin-Ono (BO) equations.

Acknowledgements

We gratefully acknowledge Wuning Zhuang (former Master student) who helped with early implementations of the Hirota method in Macsyma and Mathematica, and Ameina Nuseir (former Ph.D. student) who first computed the soliton solutions of the Kaup-Kuperschmidt equation. We also thank Andrew Cook (former REU student) who helped write Mathematica code for the simplified Hirota method. Thanks to M. Trott (Wolfram Research, Inc.) for many valuable comments and suggestions. The National Science Foundation (USA) is thanked for the support over the many years this package was developed.

References

W. Hereman and A. Nuseir, Symbolic Methods to Construct Exact Solutions of Nonlinear Partial Differential Equations, Mathematics and Computers in Simulation, 43, 13-27 (1997).
W. Hereman and U. Goktas, Symbolic computation of solitary wave solutions and solitons through homogenization of degree, in: Proceedings of the 6th International Workshop on Nonlinear and Modern Mathematical Physics (NMMP-2022), Eds. S. Manukure and W.-X. Ma, Springer Proceedings in Mathematics and Statistics, Springer Verlag, New York (2024) 61pp (in press).

Contact Information

For information about the features of PDESolitonSolutions, or in case of trouble, contact Unal Goktas (ugoktas@tamu.edu) or Willy Hereman (whereman@mines.edu). Authors will greatly appreciate hearing your comments about the PDESolitonSolutions package. Please acknowledge the paper and code if  you use it in research projects. 
Refer to the paper as follows:  W. Hereman and U. Goktas, Symbolic computation of solitary wave solutions and solitons through homogenization of degree, in `Nonlinear and Modern Mathematical Physics,’ Proceedings of the 6th International Workshop on Nonlinear and Modern Mathematical Physics (NMMP-2022), Eds.: S. Manukure and W.-X. Ma, Springer Verlag, New York (2024), 61pp. 
Refer to the code as follows:  U. Goktas and W. Hereman,  PDESolitonSolutions.m: A Mathematica package for the symbolic computation of solitary wave and soliton solutions of polynomial nonlinear PDEs using a simplified version of Hirota’s method. Dept. Appl. Math. Stat., Colorado School of Mines, Golden, Colorado (2023).   https://inside.mines.edu/~whereman