(*****************************************************************************) (* *) (* *** M A T H E M A T I C A P R O G R A M *** *) (* *) (* SYMBOLIC COMPUTATION of CONSERVED DENSITIES for SYSTEMS of *) (* NONLINEAR EVOLUTION EQUATIONS *) (* *) (* program name: condens.m *) (* *) (* purpose: computation of the conserved densities and corresponding flux *) (* with possible compatibility conditions, verification of the *) (* conservation law. *) (* *) (* input to condens.m : system of nonlinear evolution equations of any order,*) (* any degree, polynomial type, in variables x and t, *) (* only constant parameters, functions in x and t are *) (* NOT allowed as parameters *) (* *) (* u[i]_t = f(u[1],...,u[N],u[1]_{nx},...,u[N]_{nx}) *) (* with i=1,...,N; and n=1,2,... *) (* *) (* output : density and flux of desired rank (if it exists), *) (* and compatibility conditions on parameters, if applicable *) (* *) (* tested on : IBM RISC 6000, IBM Compatible PC 486, SGI Indigo2 XL *) (* *) (* language : Mathematica 3.0 and 2.2 (also versions 2.0 and 2.1) *) (* *) (* authors: Unal Goktas and Willy Hereman *) (* Department of Mathematical and Computer Sciences *) (* Colorado School of Mines *) (* Golden, CO 80401-1887, USA *) (* *) (* Version 3.2: May 4, 1998 *) (* *) (* Copyright 1998 *) (* *) (*****************************************************************************) Clear["Global`*"]; (* ------------------------------------------------------------------------- *) (*****************************************************************************) (* commentinter[]: prints a welcome message *) (*****************************************************************************) commentinter[] := Block[{}, Print["*********************************************************"]; Print[" WELCOME TO THE MATHEMATICA PROGRAM "]; Print[" by UNAL GOKTAS and WILLY HEREMAN "]; Print[" FOR THE COMPUTATION OF CONSERVED DENSITIES. "]; Print[" Version 3.2 released on May 4, 1998 "]; Print[" Copyright 1998 "]; Print["*********************************************************"] ]; (* end commentinter *) (* ------------------------------------------------------------------------- *) (*****************************************************************************) (* cls: clears the screen *) (*****************************************************************************) cls := Print["\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"]; (*****************************************************************************) (* printpage[ , ]: a subroutine for menu *) (*****************************************************************************) printpage[n_,page_] := Module[ {i,choice,lenpage}, cls; lenpage = Length[page[n]]; Print[" "]; Print[" *** MENU INTERFACE *** (page: ",n,")"]; Print["-------------------------------------------"]; For[i=1,i <= lenpage,i++, Print[Part[Part[page[n],i],1]]]; Print[" nn) Next Page"]; Print[" tt) Your System"]; Print[" qq) Exit the Program"]; Print["------------------------------------------"]; choice = Input["ENTER YOUR CHOICE: "]; Return[choice]; ]; (* end printpage *) (*****************************************************************************) (* menu: creates the menu *) (*****************************************************************************) menu := Module[ {counterpage = 1,menulist,numpages,page, choice1,control,choice2,choice3,lenmenulist,i}, menulist = { {" 1) KdV Equation (d_kdv.m)"}, {" 2) Modified KdV Equation (d_mkdv.m)"}, {" 3) Fifth Order KdV Equation-parameterized (d_5kdv.m)"}, {" 4) Fifth Order KdV Equation-Lax Case (d_5lax.m)"}, {" 5) Fifth Order KdV Equation-SK Case (d_5sk.m)"}, {" 6) Fifth Order KdV Equation-KK Case (d_5kk.m)"}, {" 7) Fifth Order KdV Equation-Ito Case (d_5ito.m)"}, {" 8) Seventh Order KdV Equation-parameterized (d_7kdv.m)"}, {" 9) Seventh Order KdV Equation-Lax Case (d_7lax.m)"}, {" 10) Seventh Order KdV Equation-SK-Ito Case (d_7ski.m)"}, {" 11) Seventh Order KdV Equation KK Case (d_7kk.m)"}, {" 12) Schamel Equation (d_scham.m)"}, {" 13) Rosenau-Hyman Equation (d_roshy.m)"}, {" 14) Hirota-Satsuma System-parameterized (d_phrsat.m)"}, {" 15) Hirota-Satsuma System (d_hirsat.m)"}, {" 16) Ito System-parameterized (d_pito.m)"}, {" 17) Ito System (d_ito.m)"}, {" 18) NLS System (real and imaginary parts) (d_nls.m)"}, {" 19) NLS System (equation and conjugate) (d_nlsalt.m)"}, {" 20) DNLS System (d_dnls.m)"}, {" 21) MVDNLS System (d_mvdnls.m)"}, {" 22) Kaup System-parameterized (d_pkaup.m)"}, {" 23) Kaup System (d_kaup.m)"}, {" 24) Kaup-Broer System (d_broer.m)"}, {" 25) Drinfel'd-Sokolov System (d_soko.m)"}, {" 26) Dispersive Long Wave System (d_disper.m)"}, {" 27) Non-dispersive Long Wave System (d_nodisp.m)"}, {" 28) 3-Component KdV System (d_3ckdv.m)"}, {" 29) 2-Component Nonlinear Schrodinger Equation (d_2cnls.m)"}, {" 30) Boussinesq System (d_bous.m)"} }; (* closes menulist *) lenmenulist = Length[menulist]; numpages = Ceiling[lenmenulist/10]; For[i = 1,i <= numpages,i++, page[i] = If[lenmenulist >= (i*10), menulist[[Table[k,{k,(i-1)*10+1,i*10}]]], menulist[[Table[k,{k,(i-1)*10+1,lenmenulist}]]]]]; choice1 = printpage[counterpage,page]; control := ( Switch[choice1, tt,Print["Make sure that you have prepared the data file for the system"]; Print["you want to test (similar to the data files we supplied)."]; choice2 = Input["If your file is ready, press 1, else 2: "]; If[choice2 === 1, choice3 = InputString["Enter the name of your data file: "]; Get[choice3], Print["Aborting the computations!"]; Print["Prepare your data file first, then start over."]; Abort[] ], nn,If[counterpage < numpages,counterpage++; choice1 = printpage[counterpage,page]; control, counterpage = 1; choice1 = printpage[1,page]; control], qq,Print["Aborting the computations!"];Abort[], 1,< List; Expand[D[f, y[x, r]] + (ReleaseHold[Thread[dummyfunc[(D[f, #1] & ) /@ Dfuncs, (Hold[Apply[Sequence, #1]] & ) /@ (Thread[{w, #1}] & ) /@ Dtimes]]] /. dummyfunc -> D) . ((-1)^#1 & ) /@ (Apply[Plus, #1] & ) /@ Dtimes] ]; (* end EulerD *) EulerD[f_, v:{(y_)[x_, r___], ___}, w:{x_, r___}] := (EulerD[f, #1, w] & ) /@ v /; If[Apply[And, (MatchQ[#1, _[Apply[Sequence, w]]] & ) /@ v], True, Message[EulerD::argx, w]]; EulerD[f_, (y_)[x_], x_] := EulerD[f, y[x], {x}]; EulerD[f_, v:{(y_)[x_], ___}, x_] := EulerD[f, v, {x}]; (*****************************************************************************) (* pdeform[]: takes an expression and prints it with subscripted derivatives *) (* and removes the functional dependence (i.e. f[x,t] -> f) *) (* Written by Tracy Otto & Tony Miller (CSM-1995). *) (* Modified by U. Goktas. *) (*****************************************************************************) pdeform[expres_] := expres /. { Derivative[n__][u[k_]][x__] :> SequenceForm[u, Subscript[k],Subscript[","], Subscript[ SequenceForm @@ Flatten[Table[#[[1]], {#[[2]]}]& /@ Transpose[{{x}, {n}}]]]], u[n_][x__] :> SequenceForm[u,Subscript[n]] }; (* end pdeform *) (*****************************************************************************) (* length[]: a modified length function *) (*****************************************************************************) length[a_] := Module[{l}, If[Head[a] === Plus,l = Length[a], If[Head[a] === Times,l = Length[a],l = 1 ] ] ;Return[l] ]; (* end length *) (*****************************************************************************) (* split[]: if free coefficients c[i] remain in the density, the density *) (* is further split in independent pieces corresponding to each *) (* coefficient c[i] *) (*****************************************************************************) split[expr_] := Module[{listofc,lenlistofc,ci,rhoi,rest}, rest=Expand[expr]; If[FreeQ[rest,c[__]], Print[" "]; Print["The density has no free coefficients."], listofc=Union[Cases[rest,c[__],12]]; lenlistofc=Length[listofc]; Print[" "]; Print["There is/are ", lenlistofc, " free coefficient(s) in the density."]; Print["These free coefficients are: ", listofc, "."]; Print["Splitting the density in independent pieces."]; Do[ ci=Part[listofc,i]; rhoi=Expand[Coefficient[rest,ci]]; rest=Expand[rest-ci*rhoi]; Print[" "]; Print["Part of the density with coefficient ", ci, ":"]; Print[" "]; Print[rhoi],{i,1,lenlistofc}]; (* end do *) Print[" "]; Print["Part of the density that has no free coefficient:"]; Print[" "]; Print[rest]; Print[" "] ] (* end if *) ]; (* end module *) (*****************************************************************************) (* normal[]: normalizes the density according to the coefficient of the *) (* first term in rho[x,t], returns normalized forms of rho[x,t] *) (* and j[x,t]. Does not divide by the coefficient if the *) (* coefficient is not a number. Returns a list with rho in the *) (* first place, and j in the second place. To be used on the list *) (* that comes out of stripper[rho[x,t]] *) (*****************************************************************************) normal[mylist_List] := Module[{leadingcoef,facexpr,denrho,rhojlist}, (* Cancel common denominators first, keep track of these with parameters *) facexpr = Factor[rho[x,t]]; denrho = Denominator[facexpr]; If[Not[NumberQ[denrho]], denrho = Factor[denrho]]; If[Not[FreeQ[Head[denrho],Times]], If[NumberQ[Part[denrho,1]] && Part[denrho,1] =!= 0, denrho = denrho/Part[denrho,1]]]; If[Not[NumberQ[denrho]], Print["WARNING!"]; Print[" "]; Print["Both rho[x,t] and j[x,t] were multiplied by: ",denrho,"."]; Print[" "]; rho[x,t] = Expand[Factor[rho[x,t]*denrho]]; j[x,t] = Expand[Factor[j[x,t]*denrho]]]; leadingcoef = Coefficient[rho[x,t],Part[mylist,1]]; If[NumberQ[leadingcoef] && leadingcoef =!= 0, rho[x,t] = Expand[rho[x,t]/leadingcoef]; j[x,t] = Expand[j[x,t]/leadingcoef]]; rhojlist = {rho[x,t],j[x,t]}; Return[rhojlist]; ]; (* end normal *) (*****************************************************************************) (* stripper[]: cancels the numerical factors of each term of the argument *) (* and puts these into a list after cancellation *) (*****************************************************************************) stripper[givenexpr_] := Module[ {lenexpr1,list = {},expr1,expr2,expr3,expr3s,lenexpr3,i,k,s}, expr1 = Expand[givenexpr]; If[Head[expr1] === Plus,lenexpr1 = Length[expr1],lenexpr1 = 1]; For[k = 1,k <= lenexpr1,k++, If[lenexpr1 == 1,expr2 = expr1,expr2 = Part[expr1,k]]; expr2 = expr2 /. c[i_] :> 1; expr3 = FactorList[expr2]; lenexpr3 = Length[expr3]; For[s = 1,s <= lenexpr3,s++, expr3s = Part[expr3,s]; If[FreeQ[expr3s,u], If[FreeQ[weightpars,Part[expr3s,1]], expr2 = expr2/Part[expr3s,1]^Part[expr3s,2]; ], s = lenexpr3+1; ]; ]; list = Append[list,expr2]; ]; Return[list]; ]; (* end stripper *) (*****************************************************************************) (* highestpair[]: a subroutine for highest order function *) (*****************************************************************************) highestpair[pairlist_] := Module[ {length,k,highestlist}, If[pairlist != {}, length = Length[pairlist]; highestlist = {Part[pairlist,1]}; If[length >= 2, For[k = 2,k <= length,k++, If[Part[pairlist,k][[1]] > Part[highestlist,Length[highestlist]][[1]], highestlist = {Part[pairlist,k]}, If[Part[pairlist,k][[1]] === Part[highestlist,Length[highestlist]][[1]], highestlist = Union[Append[highestlist,pairlist[[k]]]] ]; ]; ]; ]; Return[highestlist], If[debug, Print["The expression (given to the function that computes the"]; Print["highest order) does not have any derivatives!"]]; Return[{{0,0}}]; ] ]; (* end highestpair *) (*****************************************************************************) (* highestorder[]: returns the highest order (or the term with the highest *) (* order) of the argument *) (*****************************************************************************) highestorder[gexpr_,debug11_,debug12_,debug21_,debug22_,debug31_,debug32_] := Module[ {length1,i,newexpr1,length2,newexpr6,temporder1,temporder2,newexpr5,k, newexpr2,newexpr3,temporder3,newexpr4,orderlist1 ={},orderlist2 = {}, orderlist3 = {},highestlist1,highestlist2,highestlist3,expr}, expr = Expand[gexpr]; If[Not[FreeQ[expr,Derivative]], length1 = length[expr]; temporder1 = {}; temporder2 = {}; temporder3 = {}; newexpr5 = 1; For[i = 1,i <= length1,i++, If[Not[FreeQ[expr,Plus]],temporder1 = {}; temporder2 = {}; temporder3 = {};newexpr5 = 1]; If[length1 != 1,newexpr1 = expr[[i]],newexpr1 = expr]; If[Not[FreeQ[newexpr1,Derivative]], length2 = length[newexpr1]; For[k = 1,k <= length2,k++, If[length2 != 1,newexpr2 = newexpr1[[k]],newexpr2 = newexpr1]; If[Not[FreeQ[newexpr2,Derivative]], If[FreeQ[newexpr2,Power], newexpr3 = FullForm[newexpr2]; newexpr6 = newexpr3[[1]][[0]][[0]]; temporder1 = Append[temporder1,newexpr6[[1]]]; temporder2 = Append[temporder2,newexpr6[[2]]]; temporder3 = Append[temporder3,newexpr6[[1]]+newexpr6[[2]]]; newexpr4 = newexpr2; newexpr5 = newexpr5*newexpr2, newexpr3 = FullForm[newexpr2]; newexpr6 = newexpr3[[1]][[1]][[0]][[0]]; temporder1 = Append[temporder1,newexpr6[[1]]]; temporder2 = Append[temporder2,newexpr6[[2]]]; temporder3 = Append[temporder3,newexpr6[[1]]+newexpr6[[2]]]; newexpr4 = newexpr3[[1]]; newexpr5 = newexpr5*newexpr3[[1]]; ]; ]; ]; temporder1 = Max[temporder1]; temporder2 = Max[temporder2]; temporder3 = Max[temporder3]; orderlist1 = Append[orderlist1,{temporder1,newexpr4}]; orderlist2 = Append[orderlist2,{temporder2,newexpr4}]; orderlist3 = Append[orderlist3,{temporder3,newexpr5}]; ]; ]; ]; Which[ debug11,highestlist1 = highestpair[orderlist1]; If[debug12,Return[Table[highestlist1[[k]][[2]],{k,1,Length[highestlist1]}]], Return[highestlist1[[1]][[1]]]], debug21,highestlist2 = highestpair[orderlist2]; If[debug22,Return[Table[highestlist2[[k]][[2]],{k,1,Length[highestlist2]}]], Return[highestlist2[[1]][[1]]]], debug31,highestlist3 = highestpair[orderlist3]; If[debug32,Return[Table[highestlist3[[k]][[2]],{k,1,Length[highestlist3]}]], Return[highestlist3[[1]][[1]]]] ]; ]; (* end highestorder *) (*****************************************************************************) (* setuniformrank[]: a subroutine for the scaling function *) (*****************************************************************************) setuniformrank[list_] := Module[ {i,lengthlist,syslist={}}, lengthlist = Length[list]; For[i = 1,i < lengthlist,i++, syslist = Append[syslist,Part[list,i] == Part[list,lengthlist]]; ]; Return[syslist]; ]; (* end setuniformrank *) (*****************************************************************************) (* scaling[]: determines the scaling properties, if the equations are *) (* incompatible it will print the appropriate messages *) (*****************************************************************************) scaling[eqlist_List] := Module[ {pointsym,expr,list,syslist,msyslist = {},i,j,k,scalesol, lenmsyslist,tempscalesol = {},tempmsyslist,trouble,troubleleft, troubleright,posleft,posright,troublelist}, If[weightpars =!= {}, pointsym = Union[Table[weightpars[[i]] -> E^weight[weightpars[[i]]], {i,1,Length[weightpars]}], {u[i_Integer][x,t] :> E^weightu[i], D[u[i_Integer][x,t],{x,j_Integer}] :> E^(weightu[i]+j), D[u[i_Integer][x,t],{t,j_Integer}] :> E^(weightu[i]+j*weight[d/dt]), x -> E, t -> E^weight[d/dt]}], pointsym = {u[i_Integer][x,t] :> E^weightu[i], D[u[i_Integer][x,t],{x,j_Integer}] :> E^(weightu[i]+j), D[u[i_Integer][x,t],{t,j_Integer}] :> E^(weightu[i]+j*weight[d/dt]), x -> E, t -> E^weight[d/dt]}; ]; For[i = 1, i <= noeqs, i++, expr[i] = Part[eqlist,i]; list[i] = stripper[expr[i]]; list[i] = list[i] /. pointsym; list[i] = PowerExpand[Log[list[i]]]; syslist[i] = setuniformrank[list[i]]; msyslist = Union[msyslist,syslist[i]]; ]; scalesol = Flatten[Solve[msyslist]]; If[MemberQ[msyslist,False], Print[" "]; Print["Fatal Error! Weights are incorrect."]; Print["Check your choice(s) for the weights in the data file."]; Print["Aborting the computations!"]; CloseLog[]; Abort[], If[scalesol === {} && msyslist =!= {True}, Print[" "]; Print["In the given system there is at least one equation with terms"]; Print["of unequal rank. Scaling properties can not be determined for"]; Print["this system. The program will try to find the conflict, and,"]; Print["if successful, provide suggestions to help resolve the conflict."]; lenmsyslist = Length[msyslist]; i = 1; While[tempscalesol === {} && i <= lenmsyslist, tempmsyslist = Complement[msyslist,{Part[msyslist,i]}]; tempscalesol = Flatten[Solve[tempmsyslist]]; i++; ]; i--; If[tempscalesol =!= {}, trouble = Part[msyslist,i]; troubleleft = Part[trouble,1]; troubleright = Part[trouble,2]; For[j = 1, j <= noeqs, j++, posleft = Flatten[Position[list[j],troubleleft,{1}]]; posright = Flatten[Position[list[j],troubleright,{1}]]; If[posleft =!= {} && posright =!= {}, troublelist = Union[ Table[Part[Expand[eq[j][x,t]],posright[[k]]], {k,1,Length[posright]}], Table[Part[Expand[eq[j][x,t]],posleft[[k]]], {k,1,Length[posleft]}] ]; Print["The terms"]; Print[troublelist]; Print["in equation ",j," are incompatible. Try to introduce an"]; Print["auxiliary parameter with weight as coefficient of one of"]; Print["these terms. Aborting the computations! "]; ]; ], Print["Try to introduce auxiliary parameters with weight as coefficients"]; Print["into the system. Aborting the computations!"]; ]; CloseLog[]; Abort[], Return[scalesol]; ]; ]; ]; (* end scaling *) (*****************************************************************************) (* subroutine4[]: a subroutine for the function evaluate *) (*****************************************************************************) subroutine4[givenexpr_] := Module[ {expr1,nexpr1,freerule = {},part1,lennexpr1,i}, expr1 = Expand[givenexpr]; nexpr1 = Numerator[Factor[expr1]]; lennexpr1 = length[nexpr1]; For[i = 1,i<= lennexpr1,i++, part1 = Part[nexpr1,i]; If[MemberQ[unknownlist,part1], freerule = Union[freerule,{part1 -> 1}]]; ]; If[freerule =!= {}, Print[" "]; Print["Caution! There is a common (free) factor in the density."]; Print[" "]; Print["Setting ",freerule,"."]; ]; Return[freerule]; ]; (* end subroutine4 *) (*****************************************************************************) (* subroutine5[]: a subroutine for the function constructformrho *) (*****************************************************************************) subroutine5[nodims_,givenlist_] := Module[ {m,boundary,scale,var,len,tempwei,temppar,list,pair,ctm = 0,k,s}, list[0] = {}; For[m = 1,m <= nodims,m++, pair = Part[givenlist,m]; var = Part[pair,1]; scale = Part[pair,2]; If[scale =!= 0, If[m === 1, boundary = Floor[rhorank/scale]; list[m] = Table[{i*scale,var^i},{i,0,boundary}]; ctm = m, len = Length[list[m-1]]; list[m] = {}; For[k = 1,k <= len,k++, tempwei = Part[Part[list[m-1],k],1]; temppar = Part[Part[list[m-1],k],2]; boundary = Floor[(rhorank-tempwei)/scale]; For[s = 0,s <= boundary,s++, list[m] = Union[list[m],{{tempwei+s*scale,temppar*var^s}}]; ctm = m; ]; ]; ]; ]; ]; Return[list[ctm]]; ]; (* end subroutine5 *) (*****************************************************************************) (* constructformrho[]: returns the form of rho *) (*****************************************************************************) constructformrho[nodims_,varscalelist_] := Module[ {list,len,tempwei,temppar,inttest,formrholist = {},triviallist = {}, difflist1,difflist2,lendifflist1,i,j}, Print[" "]; Print["Starting the construction of the form of the density."]; list = subroutine5[nodims,varscalelist]; len = Length[list]; For[i = 1,i <= len,i++, tempwei = Part[Part[list,i],1]; temppar = Part[Part[list,i],2]; inttest = rhorank-tempwei; If[IntegerQ[inttest] && Not[FreeQ[temppar,u]], If[inttest === 0, formrholist = Union[formrholist,{temppar}], difflist1 = stripper[D[temppar,{x,inttest-1}]]; lendifflist1 = Length[difflist1]; For[j = 1, j <= lendifflist1,j++, difflist2 = stripper[D[Part[difflist1,j],x]]; triviallist = Union[triviallist,{Last[difflist2]}]; formrholist = Union[formrholist,Complement[difflist2,triviallist]]; ]; ]; ]; ]; lenformrho = Length[formrholist]; formrho[x,t] = 0; For[i = 1,i <= lenformrho,i++, formrho[x,t] = formrho[x,t]+c[i]*Part[formrholist,i]; ]; Print[" "]; Print["For RANK = ",rhorank,", this is the form of density rho[x,t]: "]; Print[" "]; Print[pdeform[Expand[formrho[x,t]]]]; Print[" "]; If[formrho[x,t] === 0, Print["The only density is the trivial density rho[x,t] = 0."]; Print[" "]; Print["Aborting the computations! "]; CloseLog[]; Abort[] ] ]; (* end constructformrho *) (*****************************************************************************) (* evaluate[]: computes the density and flux, with verification *) (*****************************************************************************) evaluate[expr1_,expr2_,solc_,rules_] := Module[ {test,factortest,newrhot,newformj,freerule,rhojlist}, rho[x,t] = expr1 /. solc; rho[x,t] = Expand[rho[x,t]]; newrhot[x,t] = expr2 /. Flatten[{rules,solc}]; newrhot[x,t] = Expand[newrhot[x,t]]; newformj[x,t] = (-1)*Integrate[newrhot[x,t],x]; j[x,t] = Expand[newformj[x,t]]; Print[" "]; If[FreeQ[j[x,t],Integrate] && rho[x,t] =!= 0, If[rho[x,t] =!= 0,freerule = subroutine4[rho[x,t]],freerule = {}]; rho[x,t] = rho[x,t] /. freerule; j[x,t] = j[x,t] /. freerule; test = Expand[D[rho[x,t],t]+D[j[x,t],x] /. rules]; factortest = Expand[test]; If[factortest =!= 0,factortest = Factor[factortest]]; If[factortest =!= 0,Print[" "]; Print["Automatic verification of the conservation law FAILED!"]]; Print[" "]; rhojlist = normal[stripper[rho[x,t]]]; rho[x,t] = Part[rhojlist,1]; rho[x,t] = Factor[rho[x,t]]; rho[x,t] = Expand[rho[x,t]]; j[x,t] = Part[rhojlist,2]; j[x,t] = Factor[j[x,t]]; j[x,t] = Expand[j[x,t]]; Print["*******************************************************"]; Print[" "]; Print["The normalized density rho[x,t] is:"]; Print[" "]; Print[pdeform[rho[x,t]]]; split[pdeform[rho[x,t]]]; Print[" "]; If[Length[j[x,t]] <= 20 || debug, Print["*******************************************************"]; Print[" "]; Print["The corresponding flux j[x,t] is:"]; Print[" "]; Print[pdeform[j[x,t]]], Print["*******************************************************"]; Print[" "]; Print["Flux j[x,t] has been computed. It has ",Length[j[x,t]]," terms."]; Print["Since there are more than 20 terms, j[x,t] will not be shown."]; Print["To see the flux at the end, type: pdeform[j[x,t]] or j[x,t]"] ]; Print[" "]; Print["*******************************************************"]; Print["Result of explicit verification: (rho_t + J_x) = ",factortest]; Print["*******************************************************"], Print["*******************************************************"]; Print["There is NO DENSITY of this form!"]; Print["*******************************************************"]; ] ]; (* end evaluate *) (*****************************************************************************) (* constructeqlist[]: forms the system of equations for coefficients c[i] *) (*****************************************************************************) constructeqlist[expr_] := Module[ {expr1,rules1,rules2,eqlist = {},tempexpr,expr2,templist,lentemplist,s, tempcoef,parttemplist,expr3}, rules1 = D[u[i_][x,t],{x,j_}] :> 0; rules2 = u[i_][x,t] :> 0; expr1 = Expand[expr]; tempexpr = expr1 //. rules1; tempexpr = tempexpr //. rules2; eqlist = Union[eqlist,{tempexpr}]; expr2 = expr1-tempexpr; If[expr2 =!= 0, templist = Union[stripper[expr2]]; lentemplist = Length[templist]; expr3 = expr2; For[s = 1,s <= lentemplist,s++, parttemplist = Part[templist,s]; tempcoef = Coefficient[Expand[expr2],parttemplist]; tempcoef = tempcoef //. rules1; tempcoef = tempcoef //. rules2; If[debug,expr3 = expr3 - (tempcoef*parttemplist)]; eqlist = Union[eqlist,{tempcoef}]; ]; If[debug, expr3 = Factor[expr3]; Print["Remainder of the expression during the construction of"]; Print["the system for the coefficients c[i]: ",expr3]]; If[debug && expr3 =!= 0, Print["Error in the construction of the system for the c[i]!"]; Print["Aborting the computations!"]; CloseLog[]; Abort[] ]; ]; Return[eqlist]; ]; (* end constructeqlist *) (*****************************************************************************) (* mysimplify[]: simplifies the system of equations for the c[i] *) (*****************************************************************************) mysimplify[list_] := Module[ {newlist,lennewlist,k,fterm,simplelist = {}}, Print["Starting the simplification of the system."]; Print[" "]; newlist = Complement[list,{0}]; lennewlist = Length[newlist]; For[k = 1 ,k <= lennewlist,k++, fterm = Factor[newlist[[k]]]; If[Part[fterm,0] === Times && NumberQ[Part[fterm,1]], fterm = fterm/Part[fterm,1] ]; If[Intersection[Expand[simplelist],Expand[{fterm*(-1)}]] === {}, simplelist = Union[simplelist,{fterm}]]; ]; Return[simplelist]; ]; (* end mysimplify *) (*****************************************************************************) (* analyzer[]: determines the coefficients c[i] that must be zero, used in *) (* the search for compatibility conditions *) (*****************************************************************************) analyzer[list_,len_] := Module[ {i,partlist,analyzelist = {}}, Print[" "]; Print["Starting the analysis of the system."]; Print[" "]; For[i = 1,i <= len,i++, partlist = Part[list,i]; If[Length[partlist] === 1 && MemberQ[unknownlist,partlist], analyzelist = Union[analyzelist,{partlist}], If[Length[partlist] === 2 && MemberQ[unknownlist,partlist[[2]]] && (MemberQ[parameters,partlist[[1]]] || MemberQ[weightpars,partlist[[1]]]), analyzelist = Union[analyzelist,{partlist[[2]]}] ] ]; ]; Return[analyzelist]; ]; (* end analyzer *) (*****************************************************************************) (* picklhs[]: makes a system, setting LHS == 0 *) (*****************************************************************************) picklhs[list_,k_] := Part[list,k] == 0; (*****************************************************************************) (* coefmat[]: creates the coefficient matrix of the system for the c[i] *) (*****************************************************************************) coefmat[system_List,unknowns_List] := Module[ {i,mat = {},lensys}, If[system === {} || unknowns === {}, Print["Fatal Error"], lensys = Length[system]; For[i = 1,i <= lensys,i++, mat = Append[mat,Table[Coefficient[Expand[Part[Part[system,i],1]], Part[unknowns,k]],{k,1,Length[unknowns]}]]]]; Return[mat]; ]; (* end coefmat *) (*****************************************************************************) (* main[]: solves the system for the c[i], computes the compatibility *) (* conditions (if applicable) *) (*****************************************************************************) main[utlist_] := Module[ {rules = {},trules,ressubst,formj,tformj,dtformj,comcond,k,comcondfactab, myeqlist,maineqlist,lengthpar,syscondlist1={},syscond,inputval,inputrule,sol1, lengthsol1,nrules1,par,newmaineqlist,solc,nrules2,nsolc,funclist,j,seqlist, lenseqlist,analyzelist,inputpart,inputlist,i,iii,complexanalysis,comcondfac}, rhseq[i_] := D[u[i][x,t],t] /. Part[utlist,i]; rule[i_,n_] := D[D[u[i][x,t],{x,n}],{t,1}] -> D[rhseq[i],{x,n}]; For[i = 1,i <= noeqs,i++, trules[i] = Table[rule[i,n],{n,0,highestorderrho}]; rules = Union[rules,trules[i]]; ]; If[debug, Print["These are the replacement rules from the given system: ",rules]]; (* Computation of j *) rhot[x,t] = D[formrho[x,t],t]; ressubst = rhot[x,t] /. rules; ressubst = Expand[-ressubst]; If[debug,Print["The Euler operator will be applied to: ", pdeform[ressubst]]]; funclist = Table[u[i][x,t],{i,1,noeqs}]; Print[" "]; Print["Applying the Euler operator (the variational derivative)."]; formj = EulerD[ressubst,funclist,{x,t}]; Print[" "]; Print["Finished with variational differentiation."]; tformj = formj; (* begin If point 1 *) If[Length[Cases[tformj,x_0]] === noeqs, Print[" "]; Print["End of the computations."]; evaluate[formrho[x,t],rhot[x,t],{},rules], (* else point 1 *) dtformj = tformj; dtformj = Expand[dtformj]; Print[" "]; Print["Starting the construction of the system for the c[i]."]; Print[" "]; myeqlist = {}; For[j = 1,j <= noeqs,j ++, If[Part[dtformj,j] =!= 0, myeqlist = Union[myeqlist,constructeqlist[Part[dtformj,j]]]]; ]; seqlist = mysimplify[myeqlist]; lenseqlist = Length[seqlist]; maineqlist = Flatten[Table[picklhs[seqlist,i],{i,1,lenseqlist}]]; system = MapAll[Factor,maineqlist]; If[ (Not[debug] && Length[maineqlist]>50) && (parameters==={} && weightpars==={}), Print["Since there are no parameters, and the linear system"]; Print["for the c[i] has ",Length[maineqlist], " equations, which is more than 50,"]; Print["the system will not be shown."]; Print["To see the system at the end of the computations, type: system"]; Print[" "], Print["This is the system for the c[i]:"]; Print[" "]; Print[system]; Print[" "]; Print["List of unknown coefficients c[i]:"]; Print[" "]; Print[unknownlist]; Print[" "] ]; Print["Starting the solution of this system with ", Length[maineqlist], " equations."]; If[parameters =!= {}, (* then start for parameter test *) analyzelist = analyzer[seqlist,lenseqlist]; lengthpar = Length[parameters]; Print["The system for the coefficients c[i] has ", lengthpar," parameter(s)."]; complexanalysis=lengthpar+Length[maineqlist]+Length[unknownlist]; If[complexanalysis > 20, coefmatrix = coefmat[system,unknownlist]; Print["The system of PDEs, rank and form of rho, the system for"]; Print["c[i], its coefficient matrix, and the lists of unknowns and"]; Print["parameters, will all be saved in the file worklog.m for"]; Print["further analysis. For analysis use the Mathematica functions"]; Print["Reduce, Solve, etc. Load the file worklog.m with the command"]; Print["<0]]]] /. x:(c[a_]->0)->c[a], Complement[inputlist,First[#]& /@ solc]]]; Clear[nrules1,nrules2,solc,newmaineqlist]; ] (* for *), nrules1 = ToRules[sol1]; If[comcondfac=!=True, For[k = 1,k <= lengthpar,k++, par[k]=parameters[[k]] /. nrules1; Print[" "]; Print["For ",parameters[[k]]," = ", par[k]]; ] ]; newmaineqlist = maineqlist /. nrules1; solc = Flatten[Solve[newmaineqlist /. inputrule,unknownlist]]; solc = Union[solc,inputrule]; Print[" "]; Print["The solution of the system is "]; Print[" "]; Print[solc]; nrules2 = rules /. nrules1; evaluate[formrho[x,t],rhot[x,t],solc,nrules2]; inputlist = Intersection[inputlist, Union[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a], Complement[inputlist,First[#]& /@ solc]]]; Clear[nrules1,nrules2,solc,newmaineqlist]; ] (* if *), Print[" "]; Print["The system becomes inconsistent, or the compatibility conditions"]; Print["require that one or more of the parameters are zero."]; Print["Not acceptable!"]; ] ] (* closes While *), nsolc = Flatten[Solve[maineqlist,unknownlist]]; Print[" "]; Print["Solution of the system: "]; Print[" "]; Print[nsolc]; evaluate[formrho[x,t],rhot[x,t],nsolc,rules]; Clear[solc1,solc2,solc3,nsolc]; ]; ]; (* end point 1 *) ]; (* end main *) (* ------------------------------------------------------------------------- *) (* set-up for collecting the data in a log file, transcript of computations *) logfile=""; OpenLog[filename_String] := (logfile = OpenWrite[filename]; If[logfile === $Failed, Return[]]; AppendTo[$Echo,logfile]; AppendTo[$Output,logfile];); CloseLog[] := ($Echo = Complement[$Echo,{logfile}]; $Output = Complement[$Output,{logfile}]; Close[logfile];); (* End of the auxiliary functions and procedures *) (* Start of the executable code lines *) menu; OpenLog[myfile]; commentinter[]; If[formrho[x,t] === 0, formrho[x,t] = {}]; If[Not[ListQ[formrho[x,t]]], formrho[x,t] = {formrho[x,t]}]; Print[" "]; Print["Working with the data file for the ",name,"."]; Print[" "]; For[i = 1,i <= noeqs,i++, eq[i][x,t] = Expand[eq[i][x,t]]; Print[" "]; Print["Equation ",i," of the system with ",noeqs," equation(s):"]; Print[" "]; Print[" ",pdeform[eq[i][x,t]]," = 0"]; ]; eq[x,t] = Sum[eq[i][x,t],{i,1,noeqs}]; highestordereq = highestorder[eq[x,t],False,False,False,False,True,False]; Print[" "]; Print["Highest order of derivative in the system is ",highestordereq,"."]; If[formrho[x,t] === {}, (* start point rhogiven *) eqlist = Flatten[Table[eq[i][x,t],{i,1,noeqs}]]; scalerules = scaling[eqlist]; counternegweight = 0; listfreeweights = {}; varscalelist = {}; allchoices = {}; Print[" "]; Print["Program determines the weights of the variables (and parameters)."]; For[i = 1,i <= noeqs,i++, weightu[i] = weightu[i] /. scalerules; If[i === 1, Print[" "]; Print["For the given system:"]; Print[" "]]; Print["* weight of u[",i,"] is ",weightu[i],"."]; varscalelist = Union[varscalelist,{{u[i][x,t],weightu[i]}}]; If[Not[NumberQ[weightu[i]]], If[Length[weightu[i]] === 1, listfreeweights = Union[listfreeweights,{weightu[i]}]], If[weightu[i] < 0,counternegweight++]; ]; ]; lenweightpars = Length[weightpars]; For[i = 1,i <= lenweightpars,i++, weight[Part[weightpars,i]] = weight[Part[weightpars,i]] /. scalerules; Print["* weight of ",Part[weightpars,i]," is ", weight[Part[weightpars,i]],"."]; varscalelist = Union[varscalelist, {{Part[weightpars,i],weight[Part[weightpars,i]]}}]; If[Not[NumberQ[weight[Part[weightpars,i]]]], If[Length[weight[Part[weightpars,i]]] === 1, listfreeweights = Union[listfreeweights,{weight[Part[weightpars,i]]}]], If[weight[Part[weightpars,i]] < 0,counternegweight++]; ]; ]; weight[d/dt] = weight[d/dt] /. scalerules; Print["* weight of d/dt is ",weight[d/dt],"."]; tvarscalelist = Union[varscalelist,{{"d/dt",weight[d/dt]}}]; If[Not[NumberQ[weight[d/dt]]], If[Length[weight[d/dt]] === 1, listfreeweights = Union[listfreeweights,{weight[d/dt]}]], If[weight[d/dt] < 0, Print[" "]; Print["CAUTION! Weight of d/dt is negative. Computations"]; Print["will continue, but the program is not designed for this"]' Print["case. Results may not be completely reliable."]; ]; ]; counterfreeweight = Length[listfreeweights]; If[counternegweight =!= 0, (* start of If 1 *) Print[" "]; Print["One or more of the weights are negative."]; Print["Negative weights are not allowed."]; Print["Aborting the computations!"]; CloseLog[]; Abort[], If[counterfreeweight =!= 0, If[counterfreeweight > 1, Print[" "]; Print["Two or more of weights have freedom."]; Print["Enter your values for the free weights or type Abort[];"]; Print["Your values for the free weights can be entered by typing"]; Print["`weightu[label] = val' or `weight[variable] = val' and"]; Print["putting `;' between your entries."]; Input[": "], (* else of counterfreeweight > 1 *) Print[" "]; Print["Program will try to determine CHOICES for ", listfreeweights[[1]],":"]; Print[" "]; tempvarscalelist = Sort[Table[Part[tvarscalelist,k][[2]], {k,1,Length[tvarscalelist]}]]; tempvarscalelist = Complement[tempvarscalelist,Select[tempvarscalelist,NumberQ]]; lentempvarscalelist = Length[tempvarscalelist]; Do[ (* begin Do loop *) Print["* CHOICE ",k]; Print[" "]; Print["Solving the equation: ",Part[tempvarscalelist,k]," = 1."]; Print[" "]; attempt = Flatten[Solve[Part[tempvarscalelist,k] == 1]]; If[attempt =!= {}, solfreeweight = Part[listfreeweights,1] /. attempt; While[solfreeweight <= 0, solfreeweight = solfreeweight + 1; Print["Since the weight was zero or negative,"]; Print["it was incremented with 1."]; Print[" "] ]; Print[Part[listfreeweights,1]," = ",solfreeweight]; weightrule = Part[listfreeweights,1] -> solfreeweight; tempvarscalelistval = tempvarscalelist /. weightrule; If[MemberQ[NonNegative[tempvarscalelistval],False], Print["Choice is rejected!"], varscalelistval = varscalelist /. weightrule; tvarscalelistval = tvarscalelist /. weightrule; Print[" "]; Print["List of all the variables (and parameters) with their"]; Print["weights:"]; Print[" "]; Print[tvarscalelistval]; Print[" "]; testvarscalelist = Table[Part[varscalelistval,k][[2]], {k,1,Length[varscalelist]}]; If[Union[Positive[testvarscalelist]] === {False}, Print["Since all the weights of u[i] are zero"]; Print["this choice is rejected!"]; Print[" "], allchoices = Union[allchoices,{solfreeweight}]; ]; ]; ],{k,1,lentempvarscalelist} ]; (* end do *) Print[" "]; Print["Simple POSITIVE choices are considered."]; Print[" "]; intchoices = Select[allchoices,IntegerQ]; fracchoices = Complement[allchoices,intchoices]; If[intchoices =!= {}, (* Picking the minimum integer choice *) intchoices = Min[intchoices]; choicerule = Part[listfreeweights,1] -> intchoices, (* else *) If[fracchoices =!= {}, (* Picking the minimum fractional choice *) fracchoices = Min[fracchoices]; choicerule = Part[listfreeweights,1] -> fracchoices, (* else *) Print["Enter your choice for the value of ",Part[listfreeweights,1], "by typing its value. NO semi-colon at the end!"]; choicerule = Part[listfreeweights,1]-> Input[": "]; ]; ]; varscalelistchoice = varscalelist /. choicerule; tvarscalelistchoice = tvarscalelist /. choicerule; If[Length[allchoices] > 1, Print["List of all positive choices considered: "]; Print[" "]; Print[Union[Table[Part[listfreeweights,1] -> Part[allchoices,k], {k,1,Length[allchoices]}]]]; Print[" "]; Print["Some of the weights of u[i], but not all, could be zero."]; Print["In the data file you could enter your choices in the format"]; Print["`weightu[label] = value;' or `weight[variable] = value;'."]; Print[" "]; Print["Program continues with the choice: ", choicerule]; Print[" "]; Print["corresponding to the weights:"]; Print[" "]; Print[tvarscalelistchoice]; Print[" "]; ]; varscalelist = varscalelistchoice; ] ] ]; (* end if 1 *) Print[" "]; Print["The rank of rho should be an integer multiple of the lowest weight"]; Print["of the DEPENDENT variable(s). Fractional weights are allowed."]; Print[" "]; If[ formrho[x,t] === {}, If[Not[NumberQ[rhorank]], rhorank = Input["Enter the rank of rho: "]]; Print[" "] ]; If[formrho[x,t] === {}, Print["Computation of the density (and flux) of RANK = ",rhorank]; Print[" "]; ]; nodims = noeqs + Length[weightpars]; varscalelist = Reverse[Sort[varscalelist,OrderedQ[{#1[[2]],#2[[2]]}]&]]; constructformrho[nodims,varscalelist], (* else point for rhogiven *) formrho[x,t] = Part[formrho[x,t],1]; formrho[x,t] = Expand[formrho[x,t]]; If[Head[formrho[x,t]] === Times, lenformrho = 1,lenformrho = Length[formrho[x,t]]; ]; Print[" "]; Print["The program will only test the given form of the density rho."]; Print["No determination of scaling properties."]; Print["You have given this form for rho :"]; Print[" "]; Print[pdeform[formrho[x,t]]]; Print[" "]; ]; (* end point rhogiven *) highestorderrho=highestorder[formrho[x,t],False,False,False,False,True,False]; Print["The highest-order of the terms in rho is ", highestorderrho,"."]; eqlist = Flatten[Table[eq[i][x,t],{i,1,noeqs}]]; unknownlist = Table[c[i],{i,1,lenformrho}]; parameters = Union[parameters,weightpars]; weightpars = {}; Clear[i,k,lenformrho,counternegweight,counterfreeweight,scalerules, lenweightpars,varscalelist,tempvarscalelist,tvarscalelist, attempt,listfreeweights,allchoices,varscalelistval,temparscalelistval, tvarscalelistval,intchoices,fracchoices,choicerule,weightrule, nodims,highestordereq,weightu,lentempvarscalelist,solfreeweight, varscalelistchoice,tvarscalelistchoice,testvarscalelist,name,weight]; (* Solving for u[i]_t and calculation of the x-derivatives *) mainsolforlist = Table[D[u[i][x,t],t],{i,1,noeqs}]; If[debug, Print["*************************************************"]; Print["Trying to solve the equation(s) for ",mainsolforlist,","]; Print["respectively."]; Print[" "]; Print["*************************************************"] ]; utlist = {}; For[i = 1,i <= noeqs,i++, ut[i] = Solve[eq[i][x,t] == 0,mainsolforlist[[i]]][[1]]; utlist = Flatten[Append[utlist,ut[i]]]; ]; Clear[ut,mainsolforlist,i]; main[utlist]; Print[" "]; Print["*********************** SUMMARY ***********************"]; Print[" "]; Print["Total CPU time used in the current session is ",Round[TimeUsed[]], " seconds."]; Print[" "]; Print["To see the density type: rho[x,t], or use pdeform[rho[x,t]]."]; Print[" "]; Print["To see the flux type: j[x,t], or pdeform[j[x,t]]."]; Print[" "]; Print["To split the density in independent pieces type: split[rho[x,t]],"]; Print["or split[pdeform[rho[x,t]]]."]; CloseLog[]; (* ******************************* END ************************************* *)