######################################################################
## HDequation.txt Save this file as   HDequation.txt to use it,      #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `HDequation.txt` <Enter>                     #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################



with(plots):

print(`First Written: July 2024: tested for Maple 2020 `):
print(`This Version  : July 28, 2024  `):
print():
print(`This is  HDequation.txt, a Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` Numeric and Symbolic Experiments with the Harry Dym Equation `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/HDequation.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the Story functions type: ezraSt();`):
print():




ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: EstC, GFrnA, Mnr, M2nr, NuCap, Prt, TP, Stat, STATrnA  `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` HDequation.txt: A Maple package for  experimenting with the Harry Dym Equation`):

print(`The MAIN procedures are:  HD, HDnice, HDniceS, HDt, KdV, KdVt, KdVnice, KdVniceS, PlotHD, PlotKdV  `):
print(``):

elif nargs=1 and args[1]=HD then
print(`HD(N,a0,a1,a2,t,x): Inputs the functions of t, a0,a1,a2,outputs the list of the first N+1 coefficients from x^0 to x^N`):
print(`of the solution of the Harry Dym equation f_t=f^3*f_{xxx} with initial conditions`):
print(` u(0,t)=a0(t), u_x(0,t)=a1(t), u_{xx}(0,t)=a2(t). Try:`):
print(`HD(6,t,t,t,t,x);`):

elif nargs=1 and args[1]=HDnice then
print(`HDnice(N,b): The compact way to express the first N terms of the general solutions of the Harry Dym equation with u(0,t)=a0(t),u_x(0,t)=a1(t),u_xx(0,t)=a2(t)`):
print(`where diff(ai(t),t$r) is denoted by b[i,r] for i=0,1,2. Try:`):
print(`HDnice(10,b):`):


elif nargs=1 and args[1]=HDniceS then
print(`HDniceS(N,b): The compact way to express the first N terms of the general solutions of the Harry Dym equation with u(0,t)=1,u_x(0,t)=a1(t),u_xx(0,t)=a2(t)`):
print(`where diff(ai(t),t$r) is denoted by b[i,r] for i=1,2. Try:`):
print(`HDniceS(10,b):`):


elif nargs=1 and args[1]=HDt then
print(`HDt(N,a0,a1,a2,t,x): Inputs the functions of t, a0,a1,a2,outputs the truncated power series coefficients from x^0 to x^N`):
print(`of the solution of the Harry Dym equation f_t=f^3*f_{xxx} with initial conditions`):
print(` u(0,t)=a0(t), u_x(0,t)=a1(t), u_{xx}(0,t)=a2(t). Try:`):
print(`HDt(6,t,t,t,t,x);`):



elif nargs=1 and args[1]=KdV then
print(`KdV(N,a0,a1,a2,t,x): Inputs the functions of t, a0,a1,a2,outputs the list of the first N+1 coefficients from x^0 to x^N`):
print(`of the solution of the KdV quation f_t=6*f*f_x-f_xxx with initial conditions`):
print(` u(0,t)=a0(t), u_x(0,t)=a1(t), u_{xx}(0,t)=a2(t). Try:`):
print(`KdV(6,t,t,t,t,x);`):



elif nargs=1 and args[1]=KdVnice then
print(`KdVnice(N,b): The compact way to express the first N terms of the general solutions of the KdV equation with u(0,t)=a0(t),u_x(0,t)=a1(t),u_xx(0,t)=a2(t)`):
print(`where diff(ai(t),t$r) is denoted by b[i,r] for i=0,1,2. Try:`):
print(`KdVnice(10,b):`):


elif nargs=1 and args[1]=KdVniceS then
print(`KdVniceS(N,b): The compact way to express the first N terms of the general solutions of the KdV equation with u(0,t)=1,u_x(0,t)=a1(t),u_xx(0,t)=a2(t)`):
print(`where diff(ai(t),t$r) is denoted by b[i,r] for i=1,2. Try:`):
print(`KdVniceS(10,b):`):

elif nargs=1 and args[1]=KdVt then
print(`KdVt(N,a0,a1,a2,t,x): Inputs the functions of t, a0,a1,a2,outputs the truncated power series coefficients from x^0 to x^N`):
print(`of the solution of the   f_t=6*f*f_x-f_xxx  with initial conditions`):
print(` u(0,t)=a0(t), u_x(0,t)=a1(t), u_{xx}(0,t)=a2(t). Try:`):
print(`KdVt(6,t,t,t,t,x);`):

elif nargs=1 and args[1]=PlotHD then
print(`PlotHD(a0,a1,a2,c1,c2,t,eps): inputs specific functions of t, a0,a1,a2, and numbers c1 and c2`):
print(`finds the truncated forms to 20 and 30 terms for the Harry Dym equation with initial conditions`):
print(`u(0,t)=a0, u_x(0,t)=a1, u_{xx}(0,t)=a2, and if the differce in the region {-c1<=x<=c1, -c2<=t<=c2}`):
print(`is less than epsilon, plots the one with 30 terms. Otherwise returns FAIL. Try:`):
print(`PlotHD(1,1,t,t,1/2,1/2,10^(-5));`):


elif nargs=1 and args[1]=PlotKdV then
print(`PlotKdV(a0,a1,a2,c1,c2,t,eps): inputs specific functions of t, a0,a1,a2, and numbers c1 and c2`):
print(`finds the truncated forms to 20 and 30 terms for the KdV equation with initial conditions`):
print(`u(0,t)=a0, u_x(0,t)=a1, u_{xx}(0,t)=a2, and if the differce in the region {-c1<=x<=c1, -c2<=t<=c2}`):
print(`is less than epsilon, plots the one with 30 terms. Otherwise returns FAIL. Try:`):
print(`PlotKdV(1,1,t,t,1/2,1/2,10^(-5));`):


else
 print(`There is no such thing as`, args):

fi:


end:



#HD(N,a0,a1,a2,t,x): Inputs the functions of t, a0,a1,a2,outputs the list of the first N+1 coefficients from x^0 to x^N
#of the solution of the Harry Dym equation f_t=f^3*f_{xxx} with initial conditions
#a0(t),a1(t),a2(t). Try:
#HD(6,1+t,t,t,t,x);
HD:=proc(N,a0,a1,a2,t,x) local gu,a,b,i,lu,n,gu0,kha:
gu:=a0+a1*x+a2*x^2:
lu:=[a0,a1,a2]:


for n from 3 to N do

gu:=add(diff(lu[i+1],t)*x^i,i=0..n-1)+a*x^n-(add(lu[i+1]*x^i,i=0..n-1)+a*x^n)^3*diff(add(lu[i+1]*x^i,i=0..n-1)+a*x^n,x$3):
gu0:=coeff(gu,x,n-3):
kha:=subs(solve({gu0},{a}),a):
lu:=[op(lu),kha]:
od:
lu:

end:


#HDt(N,a0,a1,a2,t,x): The first N terms of the formal solution of the Harry Dym equation with initial conditions
#f(0,t)=a0, f_x(0,t)=a1, f_{xx}(0,t)=a2, where a0,a1,a2 depend on t. Try:
#HDt(5,1+t,1,1,t,x);
HDt:=proc(N,a0,a1,a2,t,x) local lu,i,f,gu:
lu:=HD(N,a0,a1,a2,t,x):
f:=add(lu[i+1]*x^i,i=0..N):
gu:=expand(diff(f,t)-f^3*diff(f,x$3)):
gu:=add(normal(coeff(gu,x,i))*x^i,i=0..N):
if ldegree(gu,x)<N-2 then
 FAIL:
else
f:
fi:
end:



#PlotHD(a0,a1,a2,t,c1,c2,eps): inputs specific functions of t, a0,a1,a2, and numbers c1 and c2
#finds the truncated forms to 20 and 30 terms for the Harry Dym equation with initial conditions
#u(0,t)=a0, u_x(0,t)=a1, u_{xx}(0,t)=a2, and if the differce in the region {-c1<=x<=c1, -c2<=t<=c2}
#is less than epsilon, plots the one with 30 terms. Otherwise returns FAIL. Try:
#PlotHD(1,1,t,t,1/2,1/2,10^(-5));
PlotHD:=proc(a0,a1,a2,t,c1,c2,eps) local f20,f30,i,j,x,di:
f20:=HDt(20,a0,a1,a2,t,x):
f30:=HDt(30,a0,a1,a2,t,x):
di:=normal(f30-f20):
if subs(t=0,denom(di))=0 then
  RETURN(FAIL):
fi:

if max(seq(seq(evalf(abs(subs({x=c1*i/10,t=c2*j/10},f20-f30))),i=-10..10),j=-10..10))>eps then
  RETURN(FAIL):
fi:

plot3d(f30,x=-c1..c1,t=-c2..c2):


end:


#HDnice(N,b): The compact way to express the first N terms of the general solutions of the Harry Dym equation with u(0,t)=a0(t),u_x(0,t)=a1(t),u_xx(0,t)=a2(t)
#where diff(ai(t),t$r) is denoted by b[i,r] for i=0,1,2. Try:
#HDnice(10,b):
HDnice:=proc(N,b) local t,x,a0,a1,a2,gu,i:
gu:=HD(N,a0(t),a1(t),a2(t),t,x):

for i from N by -1 to 1 do
 gu:=subs(diff(a0(t),t$i)=b[0,i],gu):
od:
gu:=subs(a0(t)=b[0,0],gu):



for i from N by -1 to 1 do
 gu:=subs(diff(a1(t),t$i)=b[1,i],gu):
od:
gu:=subs(a1(t)=b[1,0],gu):



for i from N by -1 to 1 do
 gu:=subs(diff(a2(t),t$i)=b[2,i],gu):
od:
gu:=subs(a2(t)=b[2,0],gu):

gu:

end:




#HDniceS(N,b): The compact way to express the first N terms of the general solutions of the Harry Dym equation with u(0,t)=1,u_x(0,t)=a1(t),u_xx(0,t)=a2(t)
#where diff(ai(t),t$r) is denoted by b[i,r] for i=0,1,2. Try:
#HDniceS(10,b):
HDniceS:=proc(N,b) local t,x,a1,a2,gu,i:
gu:=HD(N,1,a1(t),a2(t),t,x):



for i from N by -1 to 1 do
 gu:=subs(diff(a1(t),t$i)=b[1,i],gu):
od:
gu:=subs(a1(t)=b[1,0],gu):



for i from N by -1 to 1 do
 gu:=subs(diff(a2(t),t$i)=b[2,i],gu):
od:
gu:=subs(a2(t)=b[2,0],gu):

gu:

end:



##start KdV

#KdV(N,a0,a1,a2,t,x): Inputs the functions of t, a0,a1,a2,outputs the list of the first N+1 coefficients from x^0 to x^N
#of the solution of the KdV equation f_t=6*f*f_x-f_xxxx with initial conditions
#a0(t),a1(t),a2(t). Try:
#KdV(6,1+t,t,t,t,x);
KdV:=proc(N,a0,a1,a2,t,x) local gu,a,b,i,lu,n,gu0,kha:
gu:=a0+a1*x+a2*x^2:
lu:=[a0,a1,a2]:


for n from 3 to N do

gu:=add(diff(lu[i+1],t)*x^i,i=0..n-1)+a*x^n-6*(add(lu[i+1]*x^i,i=0..n-1)+a*x^n)*diff(add(lu[i+1]*x^i,i=0..n-1)+a*x^n,x)+diff(add(lu[i+1]*x^i,i=0..n-1)+a*x^n,x$3):

gu:=expand(gu):

gu0:=coeff(gu,x,n-3):
kha:=subs(solve({gu0},{a}),a):
lu:=[op(lu),kha]:
od:
lu:

end:


#KdVt(N,a0,a1,a2,t,x): The first N terms of the formal solution of the KdV  equation with initial conditions
#f(0,t)=a0, f_x(0,t)=a1, f_{xx}(0,t)=a2, where a0,a1,a2 depend on t. Try:
#KdVt(5,1+t,1,1,t,x);
KdVt:=proc(N,a0,a1,a2,t,x) local lu,i,f,gu:
lu:=KdV(N,a0,a1,a2,t,x):
f:=add(lu[i+1]*x^i,i=0..N):
gu:=expand(diff(f,t)-6*f*diff(f,x)+diff(f,x$3)):
gu:=add(normal(coeff(gu,x,i))*x^i,i=0..N):
if ldegree(gu,x)<N-2 then
 FAIL:
else
f:
fi:
end:


#PlotKdV(a0,a1,a2,t,c1,c2,eps): inputs specific functions of t, a0,a1,a2, and numbers c1 and c2
#finds the truncated forms to 20 and 30 terms for the KdV equation with initial conditions
#u(0,t)=a0, u_x(0,t)=a1, u_{xx}(0,t)=a2, and if the differce in the region {-c1<=x<=c1, -c2<=t<=c2}
#is less than epsilon, plots the one with 30 terms. Otherwise returns FAIL. Try:
#PlotKdV(1,1,t,t,1/2,1/2,10^(-5));
PlotKdV:=proc(a0,a1,a2,t,c1,c2,eps) local f20,f30,i,j,x,di:
f20:=KdVt(20,a0,a1,a2,t,x):
f30:=KdVt(30,a0,a1,a2,t,x):
di:=normal(f30-f20):
if subs(t=0,denom(di))=0 then
  RETURN(FAIL):
fi:

if max(seq(seq(evalf(abs(subs({x=c1*i/10,t=c2*j/10},f20-f30))),i=-10..10),j=-10..10))>eps then
  RETURN(FAIL):
fi:

plot3d(f30,x=-c1..c1,t=-c2..c2):


end:


#KdVnice(N,b): The compact way to express the first N terms of the general solutions of the KdV equation with u(0,t)=a0(t),u_x(0,t)=a1(t),u_xx(0,t)=a2(t)
#where diff(ai(t),t$r) is denoted by b[i,r] for i=0,1,2. Try:
#KdVniceice(10,b):
KdVnice:=proc(N,b) local t,x,a0,a1,a2,gu,i:
gu:=KdV(N,a0(t),a1(t),a2(t),t,x):

for i from N by -1 to 1 do
 gu:=subs(diff(a0(t),t$i)=b[0,i],gu):
od:
gu:=subs(a0(t)=b[0,0],gu):



for i from N by -1 to 1 do
 gu:=subs(diff(a1(t),t$i)=b[1,i],gu):
od:
gu:=subs(a1(t)=b[1,0],gu):



for i from N by -1 to 1 do
 gu:=subs(diff(a2(t),t$i)=b[2,i],gu):
od:
gu:=subs(a2(t)=b[2,0],gu):

gu:

end:


#KdVniceS(N,b): The compact way to express the first N terms of the general solutions of the KdV equation with u(0,t)=1,u_x(0,t)=a1(t),u_xx(0,t)=a2(t)
#where diff(ai(t),t$r) is denoted by b[i,r] for i=0,1,2. Try:
#KdVniceS(10,b):
KdVniceS:=proc(N,b) local t,x,a1,a2,gu,i:
gu:=KdV(N,1,a1(t),a2(t),t,x):



for i from N by -1 to 1 do
 gu:=subs(diff(a1(t),t$i)=b[1,i],gu):
od:
gu:=subs(a1(t)=b[1,0],gu):



for i from N by -1 to 1 do
 gu:=subs(diff(a2(t),t$i)=b[2,i],gu):
od:
gu:=subs(a2(t)=b[2,0],gu):

gu:

end: