/****************************************************************/ /* */ /* *** M A C S Y M 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: HIR_SING.MAX */ /* */ /* 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 solitons). */ /* */ /* 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 VAX-8600, and on PC with */ /* Macsyma 417.125 and version 2.0 */ /* */ /* language: MACSYMA release 412, and PC Macsyma 417 and 2.0 */ /* */ /* authors: W. Hereman and W. Zhuang */ /* Dept. of Mathematical and Computer Sciences */ /* Colorado School of Mines, */ /* Golden, CO 80401-1887, USA */ /* */ /* updated: May 29, 1995 (7:00 p.m.) */ /* */ /****************************************************************/ /* Block 1: Friendly comment generator */ commentinter(name):=Block( print("/*********************************************************/"), print("/* WELCOME TO THE MACSYMA PROGRAM HIR_SING.MAX */"), 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 released on May 29, 1995 */"), print("/* Copyright 1995 */"), print("/*********************************************************/"))$ /* End of the Block commentinter(name) for the friendly interface */ /* Block 2: Tool to find exponents */ get_exponent(expr,var):=Block( if freeof(var,expr) then return('[]), if atom(expr) then return(1), if inpart(expr,0)="^" and inpart(expr,1)=var then return(inpart(expr,2)) else return('[]))$ /* End of the Block get_exponent(expr,var) to compute the exponent of var */ /* Block 3: Tool to determine list of powers */ pow(expr,var):=Block( [power_list:'[],ff:expand(expr),tmp1,tmp2,infun,inflag:true], declare(power_list,special), if not atom(ff) and part(ff,0) = "+" then infun:args(ff) else infun:[ff], for iii in infun do (tmp1:iii, if not atom(tmp1) and part(tmp1,0) = "-" then tmp1:substpart("+",tmp1,0), if atom(tmp1) or inpart(tmp1,0)#"*" then tmp1:[tmp1], tmp2 :maplist(lambda([u],get_exponent(u,var)),tmp1), tmp2 :delete('[],tmp2), if tmp2 = '[] then( if not member(0,power_list) then power_list:endcons(0,power_list)) else( if not member(first(tmp2),power_list) then power_list:endcons(first(tmp2),power_list))), power_list)$ /* End of the Block pow(exp,var) to compute the powers of var */ /* Block 4: Tool to determine highest power */ inthipow(exp,var):=Block([list,len,high], list:pow(exp,var), list:sublist(list,integerp), len:length(list), if len=0 then return(0) else if len=1 then return(part(list,1)) else (high:part(list,1), for i:2 thru len do if part(list,i) > high then high:part(list,i), return(high)))$ /* End of the Block inthipow(exp,var) to compute the highest power of var */ /* Block 5: Tool to strip off numerical factors */ stripper(expr):=Block([tempexpr,ne,de,le,pp], tempexpr:factor(expr), le:length(tempexpr), if le#0 and part(tempexpr,0)="-" then( tempexpr:-tempexpr, le:length(tempexpr)), if le#0 and part(tempexpr,0)="/" then( tempexpr:factor(tempexpr), ne:num(tempexpr), de:denom(tempexpr), if (scalarp(de) and is(de>0)) then de:1 else de:stripper(de), if (scalarp(ne) and is(ne>0)) then ne:1 else ne:stripper(ne), tempexpr:factor(ne/de), le:length(tempexpr)), if le#0 and part(tempexpr,0)="*" then( for i:1 thru le-1 do( pp:part(tempexpr,i), if (scalarp(pp) and is(pp>0)) then tempexpr:tempexpr/pp)), return(tempexpr))$ /* End of the Block stripper(expr) to strip off numerical factors */ /* Block 6: Dispersion Relation */ dispersion(P):=Block([K,L,M,i], sol:solve(P(K[i],-omega[i],L[i]),omega[i]), len:length(sol), om[i]:=ev(rhs(part(sol,1))), if len>1 then if om[i]=0 then (kill(om),om[i]:=ev(rhs(part(sol,2)))))$ /* End of the Block dispersion(P) for the dispersion law */ /* Block 7: Hirota Condition Tester for 3 and 4 solitons */ condition(P,N,MYRANDOM):=Block([ii,iii,j,ff,test,exptest,num], if MYRANDOM=true then ( if test3count>=2 then ( print("Starting the random test(s) for the existence of a", N,"soliton solution.")), n_random(N)) else print("Starting the symbolic test for the existence of a", N,"soliton solution."), term1(N):=P(sum(K[I]*S[I],I,1,N),-sum(OM[I]*S[I],I,1,N),sum(L[I]*S[I],I,1,N)), term2(N):=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,N),I,1,N), term(N):=term1(N)*term2(N), t2(I,II,F):=ev(F,[S[I]=II]), test:term(N), for J:1 thru N do test:t2(J,-1,test)+t2(J,1,test), exptest:test, if numberp(test) #true then test:factor(test), kill(K,L,OM), dispersion(P), if test #0 then( if MYRANDOM=true then ( print("The equation did not pass the random number test(s) for "), print("the existence of a", N,"soliton solution."), if (length(exptest)>=1 and part(exptest,0)="-") then exptest:-exptest, if (numberp(exptest)=true and exptest>10^25) then( print("The test with random numbers lead to the value "), print(exptest), print("which is so large that it may not be reliable on some computers"), print("(due to storage problems with large integers)."), print("Set the random-number test(s) equal to zero, and run"), print("the symbolic test(s). "), return(false)), num:numfactor(ev(test)), if (numberp(exptest)=true and exptest<10^25) then( test:float(test/num)), if numberp(test) #true then print("The condition ",test," = 0 must be satisfied.")) else ( 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)) else ( if MYRANDOM=true then ( if test3count>=2 then ( print("The equation passed the random number test(s) for the existence"), print("of a",N,"soliton solution."))) else ( print("The equation passed the symbolic test for the existence"), print("of a",N,"soliton solution.")), return(true)) )$ /* End of the Block condition(P,N,RANDOM), for n soliton test with n>=3. */ /* Block 8: Hirota Condition Tester for 1 and 2 solitons */ condition_1_2soliton(P):=Block([cond1,cond2,kk,mm,ll], cond1:ratsimp(P(0,0,0)), cond2:ratsimp(P(kk,mm,ll)-P(-kk,-mm,-ll)), if cond1 #0 or cond2 #0 then(print("There is no one- or two-soliton solution."), if cond1 #0 then( print("The bilinear operator has a constant term."), print("The condition ",factor(cond1)," = 0 should be satisfied.")), if cond2 #0 then(print("The bilinear operator is not even."), print("The condition ",factor(cond2)," = 0 should be satisfied.")), return(false)) else( print("The equation has at least a one- and two-soliton solution."), return(true)))$ /* End of the Block condition_1_2soliton(P), test for one and two-solitons */ /* Block 9: Coupling coefficients a[i,j] */ a_ij(P):=Block( a[i,j]:=-P(k[i]-k[j],-om[i]+om[j],l[i]-l[j])/P(k[i]+k[j], -om[i]-om[j],l[i]+l[j]), 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(a[i,j])) )$ /* End of the Block a_ij(P), two-soliton coupling coefficients */ /* Block 10: Coupling coefficient b[1,2,3] */ b_123(P):=Block(b[1,2,3]:a[1,2]*a[1,3]*a[2,3], print("The coefficient b[1,2,3] is calculated via the polynomial form."), print("The coefficient b[1,2,3] = ", factor(b[1,2,3])) )$ /* End of the Block b_123(P), three-soliton coupling coefficient */ /* Block 11: Hirota bilinear operators */ hirota_op(B):=Block([f,g,x,y,t,m,n], depends([f,g],[x,y,t]), Dt[n](f,g):=sum((-1)^(n-j)*n!*diff(f,t,j)*diff(g,t,n-j)/(j!*(n-j)!),j,0,n), Dx[n](f,g):=sum((-1)^(n-j)*n!*diff(f,x,j)*diff(g,x,n-j)/(j!*(n-j)!),j,0,n), Dy[n](f,g):=sum((-1)^(n-j)*n!*diff(f,y,j)*diff(g,y,n-j)/(j!*(n-j)!),j,0,n), Dxt[m,n](f,g):=sum((-1)^(m-j)*m!*sum((-1)^(n-i)*n!* diff(diff(f,x,j),t,i)*diff(diff(g,x,m-j),t,n-i) /(i!*(n-i)!),i,0,n)/(j!*(m-j)!),j,0,m), print("The equation in f corresponding to the given bilinear operator is "), print(stripper(factor(B(f,f))), " = 0"), P(KKK,MMM,LLL):=ratcoef(ev(B(%E^(KKK*X+MMM*T+LLL*Y+CST),1),diff), %E^(KKK*X+MMM*T+LLL*Y+CST)), print("For this equation the polynomial P(K,-OMEGA,L) = ", P(K,-'OMEGA,L)) )$ /* End of the Block hirota_op(B), with the Hirota operators */ /* Block 12: Conditional random number generator */ n_random(n):=Block([i,j,count,memb,memb_d,tmemb,d], memb:[0,1],memb_d:[0,1],tmemb:memb_d,i:1,count:0, while i<=n do (count:count+1,k[i]:random(25-2*hipow), if member(k[i],memb)#true then (j:i-1, while j>0 do (d:abs(k[i]-k[j]), if member(d,memb_d)#true then(memb_d:append([d],memb_d),j:j-1) else (memb_d:tmemb,j:-1 )), if j#-1 then (memb:append([k[i]],memb),tmemb:memb_d,i:i+1)), if count>=10 then( if hipow<=5 then( k[1]:2+random(3),k[2]:6+random(3),k[3]:10+random(3), k[4]:14+random(3),i:5) else( k[1]:2,k[2]:4,k[3]:7,k[4]:10,i:5))), if lfree#true then(memb:[0,1],memb_d:[0,1],tmemb:memb_d,i:1,count:0, while i<=n do (count:count+1,l[i]:random(25-2*hipow), if member(l[i],memb)#true then (j:i-1, while j>0 do (d:abs(l[i]-l[j]), if member(d,memb_d)#true then(memb_d:append([d],memb_d),j:j-1) else (memb_d:tmemb,j:-1)), if j#-1 then (memb:append([l[i]],memb),tmemb:memb_d,i:i+1 ) ), if count>=10 then( if hipow<=5 then( l[1]:2+random(3),l[2]:6+random(3),l[3]:10+random(3), l[4]:14+random(3),i:5) else( l[1]:2,l[2]:5,l[3]:7,l[4]:10,i:5)))), if test3count>=2 then ( print("Wavenumbers k[i] selected for the random number test(s): "), for ii:1 thru n do print('k[ii]," = ",k[ii])), if lfree#true then( if test3count>=2 then ( print("Wavenumbers l[i] selected for the random number test(s): "), for ii:1 thru n do print('l[ii]," = ",l[ii])) ))$ /* End of the Block n_random(n), selection of random integers */ /* for k[i], and also l[i] in the 2D case (2 space variables) */ /* Block 13: Hirota routine, main program */ hirota(B,name,n,Rtest_3soliton,Rtest_4soliton,test_3soliton,test_4soliton):= Block([nt,result2,result3,r3,r4,lfree], showtime:all, commentinter(name), if n>1 and (((Rtest_3soliton=0 and Rtest_4soliton=0) and (test_3soliton=False and test_4soliton=False))) then( print("WARNING! All the random number tests, as well as the symbolic tests"), print("for the existence of soliton solutions are ignored!"), print("Therefore, you must be sure that the equation indeed admits solitons.")), hirota_op(B),result2:condition_1_2soliton(P), if result2=true then( if n=1 then print("Starting with the construction of the one-soliton solution."), dispersion(P), print("For the ", name," equation, "), if len>1 then( print("there are ", len, " different dispersion relations."), print(sol), print("We use the dispersion relation "), print(" OMEGA[I] = ",om[i])) else( print("We use the dispersion relation "), print(" OMEGA[I] = ", om[i])), hipow:inthipow(om[1],k[1]), if hipow>=5 then ( print("WARNING!"), print("The degree of K[I] in the dispersion law is ", hipow), print("The random number test(s) may create very large integers."), print("This renders the test unreliable on some computers,"), print("(due to storage problems with large integers)."), print("Set the random-number test(s) equal to zero, and run"), print("the symbolic test(s). ")), lfree:freeof(l[i],om[i]), if lfree then print("In the expansion of f we use THETA = K X - OMEGA T + CST.") else print("In the expansion of f we use THETA = K X - OMEGA T + L Y + CST."), if Rtest_3soliton>0 then ( nt:1, test3count:1, while nt<=Rtest_3soliton+1 do ( r3:condition(p,3,true), if r3=true then (nt:nt+1,test3count:test3count+1) else nt:Rtest_3soliton+10) ) else r3:true, if test_3soliton=true then ( if r3=true then result3:condition(P,3,false) else result3:false) else ( if r3=true then result3:true else result3:false), if Rtest_4soliton>0 and result3=true then ( nt:1, test3count:1, while nt<=Rtest_4soliton+1 do ( r4:condition(p,4,true), if r4=true then (nt:nt+1, test3count:test3count+1) else nt:Rtest_4soliton+10) ) else ( if result3=true then r4:true else r4:false), if test_4soliton=true and result3=true then ( if r4=true then result4:condition(P,4,false)), if n=3 then ( if result3=true then ( print("Starting the construction of the three-soliton solution."), a_ij(P), b_123(P)) else n:2), if n=2 then ( print("Starting the construction of the two-soliton solution."), a_ij(P)), output(N)))$ /* End of the Block */ /* Hirota(B,name,N,Rtest_3soliton,Rtest_4soliton,test_3soliton,test_4soliton) */ /* for the Hirota method */ /* Block 14: Output generator */ output(N):=( if N=1 then f:1+exp(THETA[1]), if N=2 then f:1+exp(THETA[1])+exp(THETA[2])+'a[1,2]*exp(THETA[1]+THETA[2]), if N=3 then 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), if N=1 then ( print("At the end of the computations the form of the function f"), print("is available."), print("The form of f can be obtained by typing f; ."), print("The explicit form of f can be obtained by typing expression(f); .")), if N=2 then ( print("At the end of the computations the form of the function f"), print("and the coefficient a[1,2] are explicitly available."), print("The explicit factored form of a[1,2] can be obtained"), print("by typing factor(a[1,2]); "), print("The explicit forms of theta[i] and omega[i] are also available."), print("The form of f can be obtained by typing f; ."), print("The explicit form of f can be obtained by typing expression(f); .")), if N=3 then ( print("At the end of the computations the form of the function f"), print("and the coefficients a[i,j] and b[1,2,3] are available."), print("The explicit factored form of a[1,2] and b[1,2,3] can be obtained"), print("by entering factor(a[1,2]); and factor(b[1,2,3]); "), print("The explicit forms of theta[i] and omega[i] are also available."), print("The form of f can be obtained by typing f; ."), print("The explicit form of f can be obtained by typing expression(f); .")), if freeof(L[I],OM[I]) then( for i:1 thru N do (omega[i]:factor(om[i]), theta[i]:k[i]*x-om[i]*t+cst[i])) else (for i:1 thru N do (omega[i]:factor(om[i]), theta[i]:k[i]*x-om[i]*t+l[i]*y+cst[i])), if freeof(L[I],OM[I]) then (expression(f):=ev(subst( ['a[1,2]=factor(a[1,2]),'a[1,3]=factor(a[1,3]),'a[2,3]=factor(a[2,3]), 'b[1,2,3]=factor(b[1,2,3]),theta[1]=K[1]*X-OM[1]*T+CST[1], theta[2]=K[2]*X-OM[2]*T+CST[2],theta[3]=K[3]*X-OM[3]*T+CST[3]],f))) else (expression(f):=ev(subst( ['a[1,2]=factor(a[1,2]),'a[1,3]=factor(a[1,3]),'a[2,3]=factor(a[2,3]), 'b[1,2,3]=factor(b[1,2,3]),theta[1]=K[1]*X-OM[1]*T+L[1]*Y+CST[1], theta[2]=K[2]*X-OM[2]*T+L[2]*Y+CST[2], theta[3]=K[3]*X-OM[3]*T+L[3]*Y+CST[3]],f))))$ /* End of the Block output(n), returning f, its explicit form, all omega[i], */ /* and theta[i], and the coupling coefficients a[i,j] and b[1,2,3] */ /* ****************** end of the program HIROTA_SINGLE.MAX *************** */