Mathematica Code from "Glimpses of Soliton Theory"

Chapter 1 MyAnimate  and SimpleEvolver

This actually defines TWO commands: one for animating the graph of a function of x and t as a curve that evolves and one for animating a function of x, y and t as a surface that evolves:

MyAnimate[f_,{x,x0_,x1_},{y,y0_,y1_},{t,t0_,t1_},n_]:=ListAnimate[Table[Plot[f,{x,x0,x1},PlotStyle->AbsoluteThickness[2],PlotRange->{y0,y1}],{t,t0,t1,(t1-t0)/n}]]
MyAnimate[f_,{x,x0_,x1_},{y,y0_,y1_},{z,z0_,z1_},{t,t0_,t1_},n_]:=ListAnimate[Table[Plot3D[f,{x,x0,x1},{y,y0,y1},PlotRange->{z0,z1}],{t,t0,t1,(t1-t0)/n}]]

This animates the APPROXIMATE dynamics of a solution to the equation u_t = u*(u - u_x) with a given initial profile by showing what it would look like tsize units of time later.  (It is more accurate for smaller values of tsize, and you can alter this in the definition.)    Moreover, changing the line that starts RHS =. .. will allow you to change it to do the same thing to any other equation of the form u_t = (stuff not involving t derivatives).

SimpleEvolver[profile_, {x, x0_, x1_}] := (RHS = profile (profile - D[profile, x]);
    numpts = 1000;
    tsize = .1;
    profiletable = Table[{xi, (profile /. x -> xi)}, {xi, x0, x1, (x1 - x0)/numpts}];
    evolvedtable = Table[{xi, (profile + tsize RHS /. x -> xi)}, {xi, x0, x1, (x1 - x0)/numpts}];
    Show[Graphics[{Line[profiletable], Dashing[.01], Line[evolvedtable]}], Axes -> True])

Chapter 2 AnimBurgers (Burgers Equation)

AnimBurgers[f_,{x,a_,b_}]:=Module[{liszt,i,j,n,min,max},liszt[0]=Table[{i,N[f/.x->i]},{i,a-(b-a),b+(b-a),(b-a)/100}];
min=Min[Table[N[f/.x->i],{i,a,b,(b-a)/100}]];
max=Max[Table[N[f/.x->i],{i,a,b,(b-a)/100}]];
For[i=1,i≤100,i=i+1,liszt[i]=Table[{liszt[i-1][[j,1]]+.1 liszt[i-1][[j,2]],liszt[i-1][[j,2]]},{j,1,Length[liszt[0]]}]];
ListAnimate[Table[Show[Graphics[{AbsoluteThickness[2],RGBColor[0,0,1],Line[liszt[i]]}],Axes->True,AspectRatio->.5,PlotRange->{{a,b},{min-.2 (max-min),max+.2 (max-min)}}],{i,1,99,1}]]]

The command above takes a function of x (and only x) and a range of x values as its argument.  It then illustrates the dynamics of a solution with that initial profile under the evolution induced by the IVB equation by moving each point to the right at a speed proportional to its height.  (Unfortunately, we cannot write a closed formula for this in most situaions, and it often stops being the graph of a function after a short amount of time.)

AnimBurgers[1+.5 E^(-x^2),{x,-2,8}]

Chapter 3 usol and u2sol (the general 1-solution solution and a particular 2-soliton solution to the KdV equation)

As part of a homework problem, you are supposed to write a command which takes in a function of x and t and outputs the number zero if it is a solution.  I am putting a definition here just as a placeholder for this command:

KdV[v_]:=Print["You are supposed to create your own KdV[] command.  Please write one. -Alex"];

usol[x_,t_,k_]:=8k^2/(Exp[k x + k^3 t]+Exp[-k x - k^3 t])^2

The formula above solves the KdV equation for every value of k.  It is a solution that looks like a bell of height 2k^2 travelling to the left with speed k^2, as these two examples illustrate:

MyAnimate[usol[x,t,1],{x,-3,3},{y,-1,8},{t,-2,2},10]

MyAnimate[usol[x,t,2],{x,-3,3},{y,-1,8},{t,-2,2},10]

Here is a formula for a 2-soliton solution in which we see peaks that look very much like the two illustrated above, but now they appear TOGETHER in a solution...and it is NOT just a sum of the previous two solutions!

u2sol[x_,t_]:=((24 ((E^(2 (t+x)))+6 (E^(6 (3 t+x)))+4 (E^(4 (4 t+x)))+4 (E^(20 t+8 x))+(E^(34 t+10 x))))/(((1+3 (E^(2 (t+x)))+(E^(6 (3 t+x)))+3 (E^(4 (4 t+x)))))^2))

Chapter 4 Elliptic Solutions to KdV: uell

For any values of the constant parameters c, omega, k1 and k2 this is a solution to the KdV equation made out of the Weirstrass p-function for the curve with parameters k1 and k2.

uell[x_,t_,c_,omega_,k1_,k2_]:=-2 WeierstrassP[x+c t+omega,{k1,k2}]+2c/3

Chapter 5 maketau and makeu

This command just turns a function into 2 times its second logarithmic derivative:

makeu[tau_]:=Simplify[2 D[Log[tau],{x,2}]]

If tau[x,t] is a function which turns into a solution of KdV when you apply makeu, then we call it "a tau-function of the KdV equation".  As it turns out, it is often easier to write the tau-function than the solution.  For example, the solution -2/x^2 that we have encountered before comes from tau=x:

makeu[x]

maketau[vec_]:=Simplify[Wronskian[Table[Exp[vec[[i,2]] x+vec[[i,2]]^3 t]+vec[[i,1]] Exp[-vec[[i,2]] x-vec[[i,2]]^3 t],{i,1,Length[vec]}],x]]

maketau[{{a_,k_}}]:=Exp[k x+k^3 t]+a Exp[-k x-k^3 t]

The command above takes a list of the form { {c1,k1}, {c2,k2}, {c3,k3},....,{cn,kn} } and turns it into a tau-function for a KdV solution (or, if you do not choose wisely, into the number zero).  The resulting solution will have n different things moving to the left (either soliton-like peaks or vertical asymptotes that go to negative infinity).  The parameters k1,..,kn determine the SPEEDS of those objects (they will have speeds k1^2, k2^2, ..., kn^2) and the c1,...,cn parameters will control both the positions and whether they are smooth soliton peaks or "singularities' with vertical asymptotes.

For example:

makeu[maketau[{{1,1},{1,2}}]]

is the solution shown in the book as part of example 6.1.  It has one smooth peak and one vertical asymptote.  However, changing c2 from 1 to -1 makes it into the pure 2-soliton solution that we saw back in 4.5.

Chapter 6 Differential Operators

applying differential operators

odoapply[ll_,f_]:=Module[{i},
Simplify[Sum[ Coefficient[ll,Dd,i] D[f,{x,i}] ,{i,0,Exponent[Collect[ll,Dd],Dd]}]]]

Type the differential operator L as a polynomial in Dd with coefficients that are functions of x.  Then it can be applied to the function f using odoapply[L,f]

LL=Dd^2-2/x^2

odoapply[LL,x Exp[2x]]

Multiplying Differential Operators

Execute these definitions:

odomult[L_,M_]:=Module[{i,f,rem},
rem=odoapply[L,odoapply[M,f'[x]]];
Sum[Simplify[Coefficient[rem,D[f'[x],{x,i}]]] Dd^i,{i,0,Exponent[Collect[L*M,Dd],Dd]}]]

odosimp[LL_]:=Module[{n,i,L,outL},(
L=Collect[LL,Dd];
n=Exponent[L,Dd];
For[i=0;outL=0,i<n+1,i=i+1,outL=outL+Simplify[Coefficient[L,Dd,i]] Dd^i];
outL)]

Now you can define differential operators using "Dd" for "partial":

LL=Dd^2+Sin[x]
MM=Dd^3-x Dd + Exp[x]

and you can compute their products with "odomult[ , ]":

odomult[LL,MM]

odomult[MM,LL]

Then, of course, their commutator is

comm=odomult[LL,MM]-odomult[MM,LL]

But, notice that this is not written in a nice form (there is more than one Dd^2 term!) and so you can simplify it with "odosimp[]":

odosimp[comm]

making a differential operator with specified kernel

findK[vec_]:=odosimp[Wronskian[Append[vec,ff'[x]],x]/Wronskian[vec,x]/.{D[ff[x],{x,i_}]->Dd^(i-1),Wronskian[{f_},x]->f}]

The command above takes a list of functions of x as an argument and gives you back the lowest order monic differential operator which has all of them in its kernel.  For instance:

K=findK[{x,E^x}]

This means that d^2 -x/(x-1) d + 1/(x-1) is a differential operator that gives zero when applied to x or to E^x.  We can check this as follows:

odoapply[K,x]

odoapply[K,Exp[x]]

Aha, it worked!

Chapter 10: KP

This tests to see if u is a solution to the KP equation, returning "0" if it is:

KP[u_]:=Simplify[-3/4*D[u,{y,2}]+D[D[u,t]-1/4*(6*u*D[u,x]+D[u,{x,3}]),x]]

This tests to see if tau is a solution to the bilinear KP equation, returning "0" if it is :

bilinearKP[tau_]:=Module[{newtau1,newtau2,prod,x1,x2,y1,y2,z1,z2},Print["Applying bilinear KP equation to ",tau];
newtau1=tau /. {x->x1-y1,y->x2-y2,t->x3-y3};
newtau2=tau /. {x->x1+y1,y->x2+y2,t->x3+y3};
prod=newtau1*newtau2;
Simplify[D[prod,{y1,4}]-4D[D[prod,y3],y1]+3D[prod,{y2,2}] /.{x1->x,
x2->y,x3->t,y1->0,y2->0,y3->0}]]

The command phi[lambda,n] makes a nicely weighted function by differentiating Exp[xz+yz^2+tz^3] n times with respect to z and evaluating it at z=lambda:

phi[lambda_,n_]:=Simplify[D[Exp[z x + z^2 y + z^3 t],{z,n}]/.z->lambda]

Chapter 11: Gamma_24 and the Geometric Origins of Bilinear KP

Could the Bilinear KP Equation be an algebraic equation "in disguise"?

Omega[i_,j_]:=Wronskian[{phi[0,i-1],phi[0,j-1]},x]

gentau=c12 Omega[1,2]+c13 Omega[1,3]+ c14 Omega[1,4]+c23 Omega[2,3]+ c24 Omega[2,4]+ c34 Omega[3,4]

When we plug the linear combination gentau into Bilinear KP we see the equation for a Grassmann Cone:

bilinearKP[gentau]

Can we isolate the coefficients of gentau using derivatives?  Yes, if we make differential operators out of the polynomials Omega[i,j] themselves:

P[i_,j_]:=Omega[i,j]/.{x->ddx,y->1/2 ddy,t->1/3 ddt}

If you apply the differential operator:

P[2,3]

to gentau and then set (x,y,t)=(0,0,0) you get

D[gentau,{x,2}]-D[gentau,y]/.{x->0,y->0,t->0}

This isolates the coefficient c23 by applying P[2,3]!  Can we get c24 the same way?

P[2,4]

-2 D[gentau,t]+2 D[gentau,{x,3}]/.{x->0,y->0,t->0}

Yes!

We can define  a command that will give us these differential expressions that isolate the coefficients

c[i_,j_]:=Expand[P[i,j] ddx^2 ddy^2 ddt^2 tau[x,y,t]]/.ddx^ii_ ddy^jj_ ddt^kk_ tau[x,y,t]->Derivative[ii-2,jj-2,kk-2][tau][x,y,t]

c[2,3]

c[3,4]

And we substitute those into the Plucker Relation (just as we substituted D[f/2,{x,2}], D[f/6,{y,3}],... into the equation for the sphere)

plucker=c[1,2]c[3,4]-c[1,3]c[2,4]+c[1,4]c[2,3]

COULD THIS BE BILINEAR KP?!?

bilinkp=bilinearKP[tau[x,y,t]]

Simplify[2plucker-bilinkp]

Chapter 13: Gamma_kn and the Bilinear KP Hierarchy

The differential operators do pull off the coordinates

phi[i_]:=D[Exp[t1 z + t2 z^2 + t3 z^3 + t4 z^4 +t5 z^5],{z,i-1}]/.z->0

phi[2]

phi[3]

phi[4]

omega[i1_,i2_,i3_]:=Wronskian[{phi[i1],phi[i2],phi[i3]},t1]

omega[1,2,4]

Expand[843 omega[1,2,4]+953 omega[1,3,5] + 8018 omega[3,4,5]]

linearcomb[t1_,t2_,t3_,t4_,t5_]:=5058 t1+15248 (t1^3)+16036 (t1^6)-96216 (t1^4) t2+577296 (t1^2) (t2^2)+1154592 (t2^3)-45744 t3-384864 (t1^3) t3-2309184 t1 t2 t3+2309184 (t3^2)+1154592 (t1^2) t4-2309184 t2 t4

Pappliedtotau[i1_,i2_,i3_]:=Expand[((i1-1)!(i2-1)!(i3-1)!)^(-2)(omega[i1,i2,i3]/.{t1->ddt1,t2->1/2 ddt2,t3->1/3 ddt3,t4->1/4 ddt4,t5->1/5 ddt5}) ddt1^2 ddt2^2 ddt3^2 ddt4^2 ddt5^2 tau[t1,t2,t3,t4,t5]]/.
ddt1^j1_ ddt2^j2_ ddt3^j3_ ddt4^j4_ ddt5^j5_ tau[t1,t2,t3,t4,t5]->Derivative[j1-2,j2-2,j3-2,j4-2,j5-2][tau][t1,t2,t3,t4,t5]

Pappliedtotau[1,2,5]

%/.tau->linearcomb

Simplify[%/.{t1->0,t2->0,t3->0,t4->0,t5->0}]

Pappliedtotau[1,2,4]

%/.tau->linearcomb

Simplify[%/.{t1->0,t2->0,t3->0,t4->0,t5->0}]

Pappliedtotau[1,3,5]

%/.tau->linearcomb

Simplify[%/.{t1->0,t2->0,t3->0,t4->0,t5->0}]

Pappliedtotau[3,4,5]

%/.tau->linearcomb

Making Bilinear PDEs from Plucker Relations with k=3 and n=6

The command below applies the operator P[i1,i2,i3] to the arbitrary function tau[t1,t2,t3,t4,t5]:

c[i1_,i2_,i3_]:=Expand[((i1-1)!(i2-1)!(i3-1)!)^(-2)(omega[i1,i2,i3]/.{t1->ddt1,t2->1/2 ddt2,t3->1/3 ddt3,t4->1/4 ddt4,t5->1/5 ddt5}) ddt1^2 ddt2^2 ddt3^2 ddt4^2 ddt5^2 tau[t1,t2,t3,t4,t5]]/.
ddt1^j1_ ddt2^j2_ ddt3^j3_ ddt4^j4_ ddt5^j5_ tau[t1,t2,t3,t4,t5]->Derivative[j1-2,j2-2,j3-2,j4-2,j5-2][tau][t1,t2,t3,t4,t5]

c[2,3,4]

Now we substitute those into the two (of many) Plucker relations we found for Gamma_36:

longplucker=Expand[-c[1,2,3]c[4,5,6]-c[1,3,4]c[2,5,6]+c[1,3,5]c[2,4,6]-c[1,3,6]c[2,4,5]]

shortplucker=Expand[-c[1,2,3]c[2,4,5]+c[1,2,4]c[2,3,5]-c[1,2,5]c[2,3,4]]

producing solutions as Wronskian determinants

newphi[lambda_,n_]:=Simplify[D[Exp[t1 z + t2 z^2 + t3 z^3 + t4 z^4 + t5 z^5],{z,n}]/.z->lambda]

trytau=Wronskian[{newphi[1,0]+newphi[-1,0],newphi[0,5],newphi[0,0]+newphi[1,0]},t1]

tau[t1_,t2_,t3_,t4_,t5_]:=2 (E^(t2+t4)) ((E^(t2+t4)) ((t1^5)+20 (t1^3) (-1+t2)+60 (t1^2) t3+60 t1 ((-2+t2) t2+2 t4)+120 ((-1+t2) t3+t5))-5 ((t1^4)+12 (t1^2) t2+24 t1 t3+12 ((t2^2)+2 t4)) Cosh[t1+t3+t5]+20 ((t1^3)+6 t1 t2+6 t3) Sinh[t1+t3+t5])

Simplify[longplucker]

Simplify[shortplucker]