/**************************************************************/ /* */ /* *** 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_SYS.MAX */ /* */ /* purpose: calculate the DISPERSION LAW, the CONDITIONS for */ /* the existence of up to a four soliton solution and the */ /* explicit CONSTRUCTION of the one, two and the three */ /* SOLITON SOLUTION for single evolution or wave equations, */ /* of MKdV type. */ /* */ /* computers: tested on VAX-8600 and on PC */ /* language: MACSYMA release 412 */ /* and PC Macsyma 417.125 and 2.0 */ /* */ /* authors: W. Hereman and W. Zhuang */ /* Dept. of Mathematical and Computer Sciences */ /* Colorado School of Mines, */ /* Golden, CO 80401-1886, USA */ /* */ /* updated: May 29, 1995 */ /* */ /**************************************************************/ COMMENTINTER(name):=BLOCK( PRINT("/*********************************************************/"), PRINT("/* WELCOME TO THE MACSYMA PROGRAM HIR_SYS */"), PRINT("/* BY WILLY HEREMAN AND WUNING ZHUANG */"), PRINT("/* FOR THE CALCULATION OF SOLITON SOLUTIONS */"), PRINT("/* OF THE ", NAME," EQUATION */"), PRINT("/* WITH HIROTA'S METHOD */"), PRINT("/* Version 1.0 released on May 29, 1995 */"), PRINT("/* Copyright 1995 */"), PRINT("/*********************************************************/"))$ dispersion(P1,check_coefficients,n,name):=block( m:rhs(part(solve(p1(k,m,l),m),1)), m1:ratsubst(l1,l,ratsubst(k1,k,m)), m2:ratsubst(l2,l,ratsubst(k2,k,m)), m3:ratsubst(l3,l,ratsubst(k3,k,m)), m4:ratsubst(l4,l,ratsubst(k4,k,m)), if check_coefficients=true then( for i:1 thru 4 do om[i]:ratsubst(l[i],l,ratsubst(k[i],k,m)), f1:sum(%i*h(i,x,y,t),i,1,n), g1:sum(-%i*h(i,x,y,t),i,1,n), if freeof(l,m) then gradef(h(i,x,y,t),0,k[i]*h(i,x,y,t),0,-om[i]*h(i,x,y,t)) else gradef(h(i,x,y,t),0,k[i]*h(i,x,y,t),l[i]*h(i,x,y,t), -om[i]*h(i,x,y,t)) ), if freeof(l,m) then print("In the expansions THETA = K X - OMEGA T + CONSTANT.") else print("In the expansions THETA = K X - OMEGA T + L Y + CONSTANT."), print("For the ", name," equation, the dispersion relation is "), print(" OMEGA =",m) )$ /* End of the block for the dispersion law */ condition_4soliton(P1,P2):= block([cond41,cond42,temp,s1,s2,s3,s4], print("Starting the test for the existence of a four soliton solution"), cond41:0, cond42:0, for s1:-1 thru 1 step 2 do( for s2:-1 thru 1 step 2 do( for s3:-1 thru 1 step 2 do( for s4:-1 thru 1 step 2 do( temp:P2(s2*k2-s1*k1,s2*m2-s1*m1,s2*l2-s1*l1) *P2(s3*k3-s1*k1,s3*m3-s1*m1,s3*l3-s1*l1) *P2(s3*k3-s2*k2,s3*m3-s2*m2,s3*l3-s2*l2) *P2(s4*k4-s1*k1,s4*m4-s1*m1,s4*l4-s1*l1) *P2(s4*k4-s2*k2,s4*m4-s2*m2,s4*l4-s2*l2) *P2(s4*k4-s3*k3,s4*m4-s3*m3,s4*l4-s3*l3), cond41:ratsimp(cond41) +P1(s1*k1+s2*k2+s3*k3+s4*k4,s1*m1+s2*m2+s3*m3+s4*m4, s1*l1+s2*l2+s3*l3+s4*l4) *sin((s1+s2+s3+s4)*%pi/2)*temp, cond42:ratsimp(cond42) +P2(s1*k1+s2*k2+s3*k3+s4*k4,s1*m1+s2*m2+s3*m3+s4*m4, s1*l1+s2*l2+s3*l3+s4*l4) *cos((s1+s2+s3+s4)*%pi/2)*temp )))), cond41:ratsimp(cond41), cond42:ratsimp(cond42), if cond41=0 and cond42=0 then( print("There also exists a four soliton solution!"), return(true)) else( print("There is no four soliton solution!"), if cond41 # 0 then print("The condition",factor(cond41),"=0 should be satisfied"), if cond42 # 0 then print("The condition",factor(cond42),"=0 should be satisfied"), return(false)) )$ condition_3soliton(P1,P2):= block([cond31,cond32,temp,s1,s2,s3], print("Starting the test for the existence of a three soliton solution!"), /* Three Soliton Condition */ cond31:0, cond32:0, for s1:-1 thru 1 step 2 do( for s2:-1 thru 1 step 2 do( for s3:-1 thru 1 step 2 do( temp:P2(s2*k2-s1*k1,s2*m2-s1*m1,s2*l2-s1*l1) *P2(s3*k3-s1*k1,s3*m3-s1*m1,s3*l3-s1*l1) *P2(s3*k3-s2*k2,s3*m3-s2*m2,s3*l3-s2*l2), cond31:cond31 +P1(s1*k1+s2*k2+s3*k3,s1*m1+s2*m2+s3*m3,s1*l1+s2*l2+s3*l3) *sin((s1+s2+s3)*%pi/2)*temp, cond32:cond32 +P2(s1*k1+s2*k2+s3*k3,s1*m1+s2*m2+s3*m3,s1*l1+s2*l2+s3*l3) *cos((s1+s2+s3)*%pi/2)*temp ))), cond31:ratsimp(cond31), cond32:ratsimp(cond32), if cond31=0 and cond32=0 then( print("There also exists a three soliton solution!"), return(true)) else( print("There is no three soliton solution!"), if cond31 # 0 then print("The condition",factor(cond31),"=0 should be satisfied"), if cond32 # 0 then print("The condition",factor(cond32),"=0 should be satisfied"), return(false)) )$ /* End of the block for the condition_3soliton */ condition_1_2soliton(P1,P2):=block([cond1,cond21,cond22,kk,mm,ll], cond1:ratsimp(p2(0,0,0)), cond21:ratsimp(p1(kk,mm,ll)+p1(-kk,-mm,-ll)), cond22:ratsimp(p2(kk,mm,ll)-p2(-kk,-mm,-ll)), if cond1 #0 or cond21 #0 or cond22 #0 then( print("There is not a one or two soliton solution!"), if cond1 #0 then( print("The bilinear operator contains a constant term!"), print("The condition",factor(cond1),"=0 should be satisfied")), if cond21 #0 then( print("The bilinear operator is not odd!"), print("The condition",factor(cond21),"=0 should be satisfied")), if cond22 #0 then( print("The bilinear operator is not even!"), print("The condition",factor(cond22),"=0 should be satisfied")), return(false)) else( print("There is at least a one or two soliton solutions!"), return(true)) )$ /* End of the block for the condition_single */ construct_2soliton(B2,n):=block([b2onl2,tf2,tg2,f2,g2,i,j], f2:0, g2:0, for i:1 thru n do ( tf2:sum(aa[i,j]*%i*h(i,x,y,t)*%i*h(j,x,y,t),j,i+1,n), f2:ev(f2+tf2), tg2:sum(aa[i,j]*(-%i)*h(i,x,y,t)*(-%i)*h(j,x,y,t),j,i+1,n), g2:ev(g2+tg2) ), b2onl2:ratsimp(ev(B2(f2,1)+B2(f1,g1)+B2(1,g2),diff,eval)), eqtwo[1,2]:ratcoef(b2onl2,h(1,x,y,t)*h(2,x,y,t),1), aa[1,2]:ratsimp(rhs(part(solve(eqtwo[1,2],aa[1,2]),1))), if n=3 then( aa[1,3]:ratsubst(k[3],k[2],ratsubst(l[3],l[2],aa[1,2])), aa[2,3]:ratsubst(k[2],k[1],ratsubst(l[2],l[1],aa[1,3]))), print("The coefficient", aa[i,j]," is calculated via the bilinear operator."), print("The coefficient", aa[i,j]," = ",factor( ratsubst(l[i],l[1],ratsubst(l[j],l[2], ratsubst(k[i],k[1],ratsubst(k[j],k[2],aa[1,2])))))) )$ /* end of the block for the construct_2soliton */ construct_3soliton(B1,n):=block( [b1onl3,teq_right,eq_right,eq_left,i,j,s,i1,j1,s1], /* construct f3 */ b1onf1g2: ev(b1(%i*h(i,x,y,t),aa[j,s]*(-%i)*h(j,x,y,t)*(-%i)*h(s,x,y,t)),diff,eval), b1onf2g1: ev(b1(aa[j,s]*%i*h(j,x,y,t)*%i*h(s,x,y,t),(-%i)*h(i,x,y,t)),diff,eval), teq_right:ratcoef(b1onf1g2+b1onf2g1,h(i,x,y,t)*h(j,x,y,t)*h(s,x,y,t),1), b1onf3:ev(b1(%i*h(1,x,y,t)*%i*h(2,x,y,t)*%i*h(3,x,y,t),1),diff,eval), b1ong3:ev(b1(1,(-%i)*h(1,x,y,t)*(-%i)*h(2,x,y,t)*(-%i)*h(3,x,y,t)),diff,eval), eq_left:ratsimp(ratcoef(b1onf3+b1ong3,h(1,x,y,t)*h(2,x,y,t)*h(3,x,y,t),1)), eq_right:0, for i1:1 thru n do ( for j1:1 thru n do ( for s1 : j1+1 thru n do ( if i1 # j1 and i1 # s1 then ( eq_right:eq_right+ratsimp(ev(teq_right,[i=i1,j=j1,s=s1])) )))), bb[1,2,3]:-eq_right/eq_left, bb[1,2,3]:factor(bb[1,2,3]), print("The coefficient",'bb[1,2,3]," is calculated via the bilinear operator."), print("The coefficient",'bb[1,2,3]," = ",bb[1,2,3]) )$ /* end of the block for construct_3soliton */ a_ij(P2):=block( a[1,2]:P2(k1-k2,m1-m2,l1-l2)/P2(k1+k2,m1+m2,l1+l2), print("The coefficient", a[i,j]," is calculated via the polynomial form."), print("The polynomial is P2 =", p2(k,'m,l)), print("The coefficient", a[i,j]," = ", factor(ratsubst(l[i],l1,ratsubst(l[j],l2, ratsubst(k[i],k1,ratsubst(k[j],k2,a[1,2])))))) )$ /* end of the block a_ij */ check_a():=block([check], check:ratsimp(aa[1,2]- ratsubst(l[1],l1,ratsubst(l[2],l2, ratsubst(k[1],k1,ratsubst(k[2],k2,a[1,2]))))), if check=0 then ( print("The coefficients", a[i,j]," and", aa[i,j]," are the same!"), return(true)) else( print("The coefficients", a[i,j]," and", aa[i,j]," are different!"), return(false)) )$ /* end of the block check_a */ b_123(P2):=block( a[1,3]:P2(k1-k3,m1-m3,l1-l3)/P2(k1+k3,m1+m3,l1+l3), a[2,3]:P2(k2-k3,m2-m3,l2-l3)/P2(k2+k3,m2+m3,l2+l3), 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 */ check_b():=block([check], check:ratsimp(bb[1,2,3]-ratsubst(l[1],l1,ratsubst(l[2],l2,ratsubst(l[3],l3, ratsubst(k[1],k1,ratsubst(k[2],k2,ratsubst(k[3],k3,b[1,2,3]))))))), if check=0 then ( print("The coefficients ",'b[1,2,3],"and ",'bb[1,2,3],"are the same!"), return(true)) else( print("The coefficients ",'b[1,2,3],"and ",'bb[1,2,3],"are different!"), return(false)) )$ /* end of the block check_b */ Hirota_op(B1,B2):=block( depends([f,g,h],[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), P1(K,M,L):=RATCOEF(EV(B1(1,%E^(L*Y+K*X-M*T+cst)),DIFF),%E^(L*Y+K*X-M*T+cst)), P2(K,M,L):=RATCOEF(EV(B2(1,%E^(L*Y+K*X-M*T+cst)),DIFF),%E^(L*Y+K*X-M*T+cst)))$ /* end of the block for the Hirota operation */ /* Begin the block of the Hirota Method */ Hirota(B1,b2,name,n,test_for_3soliton,check_coefficients,test_for_4soliton):= block( showtime:all, commentinter(name), Hirota_op(B1,B2), result2:condition_1_2soliton(P1,P2), /* check the even and no constant term bilinear operator */ if result2=true then( if n=1 then print("Starting the construction of the one soliton solution!"), dispersion(P1,check_coefficients,n,name), if test_for_3soliton=true then result3:condition_3soliton(P1,P2) else result3:true, if test_for_4soliton=true and result3=true then condition_4soliton(P1,P2), if n=3 then ( if result3=true then ( print("Starting the construction of the three soliton solution!"), a_ij(P2), b_123(P2)) else n:2 ), if n=2 then ( print("Starting the construction of the two soliton solution!"), a_ij(P2)), output(n), if check_coefficients=true then( print("Starting the verification of the coefficients"), if n>1 then( construct_2soliton(B2,n), check_a()), if n>2 then ( construct_3soliton(B1,n), check_b()) ) ) )$ /* end of the block for the Hirota Method */ output(n):=([f,g], if n=1 then ( f:1+%i*exp(THETA[1]), g:1-%i*exp(THETA[1])), if n=2 then ( f:1+%i*exp(THETA[1])+%i*exp(THETA[2])-'a[1,2]*exp(THETA[1]+THETA[2]), g:1-%i*exp(THETA[1])-%i*exp(THETA[2])-'a[1,2]*exp(THETA[1]+THETA[2])), if n=3 then ( f:1+%i*exp(THETA[1])+%i*exp(THETA[2])+%i*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]*%i*exp(THETA[1]+THETA[2]+THETA[3]), g:1-%i*exp(THETA[1])-%i*exp(THETA[2])-%i*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]*%i*exp(THETA[1]+THETA[2]+THETA[3])), print("The function f =", f), print("The function g =", g) )$