(* Last updated: August 4, 2007 at 00:55 at home in Christchurch *) (* File: c0804-07.m *) (* Code now works under Mathematica versions 4 and 5 *) (* Bug fixes see: WH 08/04/2007 *) (* Based in file : c060717h.m *) (* Previously updated: 07/17/2006 by Hereman at office (Golden) at 15:00 *) (* Provided the option of using the old algorithm or the new algorithm *) (* to investigate the compatibility conditions. *) (* See useoldalgorithm and usenewalgorithm *) (* Rewrote the normal function to normalize the density. *) (* Add newanalyzer to dynamically clean up the system. *) (* New method of selecting the c[i] for the case were there are parameters *) (* in the system. *) (* The code now determines the listzeroc of the c[i] that must be zero, *) (* computes listnonzeroc as the complement relative to the list unknowns, *) (* and set these c[i] equal to one (one at the time). *) (* Added a function to check is the next rho is independent from the *) (* previously computed rho *) (* Keep track of all the independent rhos in listallrhos. *) (* Similarly, keep track of independent fluxes in a listallfluxes *) (* Changes marked by 07/03/2006 WH, also look for word debugflags *) (* Purpose: Put print statements in the routines that analyze and solve *) (* the system for the constants c[i] and find the compatibility conditions *) (* ********************************************************************** *) (* Based on condens.m from 06/10/2003 *) (* 06/10/2003 Added a LogicalExpand to Reduce in solver for the c[i] *) (* Replaced worklog.m by workcond.m *) (* Always saves the algebraic system etc in workcond.m now *) (* 02/23/2002 Hereman *) (*****************************************************************************) (* *) (* *** 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 : PCs with Windows XP and workstations with Linux *) (* *) (* language : Mathematica v. 5.0, backward compatible with v. 4.0 *) (* *) (* authors: Unal Goktas and Willy Hereman *) (* Department of Mathematical and Computer Sciences *) (* Colorado School of Mines *) (* Golden, CO 80401-1887, USA *) (* *) (* Last updated: August 4, 2007 *) (* Based on version 3.2: May 4, 1998 *) (* *) (* Copyright 1998-2007 *) (* *) (*****************************************************************************) (* WH 08/04/2007 print statement added *) Print["Loading package condens.m of August 4, 2007."]; Clear["Global`*"]; (* WH 08/04/2007 set debug flags *) (* 07/03/2006 WH debugflags *) debugstripper = False; debugmystripper = False; debugroutine4 = False; debugnewanalyze = False; debugsolve = False; debugnormal = False; debugevaluate = False; usenewalgorithm = True; useoldalgorithm = False; (* ------------------------------------------------------------------------- *) (*****************************************************************************) (* 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[" Last updated: August 4, 2007 "]; Print[" Based on version 3.2 released on May 4, 1998 "]; Print[" Copyright 1998-2007 "]; 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[" *** 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["The density has no free coefficients."], listofc=Union[Cases[rest,c[__],12]]; lenlistofc=Length[listofc]; 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["Part of the density with coefficient ", ci, ":"]; Print[rhoi],{i,1,lenlistofc}]; (* end do *) Print["Part of the density that has no free coefficient:"]; Print[rest] ] (* end if *) ]; (* end module *) (* new normal function *) (*****************************************************************************) (* 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[ {expandrho,firsttermrho,firstpartmylist,normalizationcoef, normalizationcoefcipart}, expandrho = Expand[rho[x,t]]; If[debugnormal, Print["At pt. NORMAL01 in normal, expanded form of rho, expandrho: "]; Print[expandrho] ]; If[Head[expandrho] === Plus, (* then 1 *) firsttermrho = Part[expandrho,1]; If[debugnormal, Print["At pt. NORMAL02 in normal, rho is a sum, firsttermrho: "]; Print[firsttermrho] ], (* else *) firsttermrho = expandrho; If[debugnormal, Print["At pt. NORMAL03 in normal, rho is not a sum, firsttermrho: "]; Print[firsttermrho] ] ]; (* end if 1 *) If[Head[firsttermrho] === Times, (* then 2 *) firstpartmylist=Part[mylist,1]; If[debugnormal, Print["At pt. NORMAL04 in normal, Part[mylist,1], firstpartmylist:"]; Print[firstpartmylist] ]; normalizationcoef = firsttermrho/firstpartmylist; (* WH 08/04/2007 taking the part of normalizationcoef that is independent *) (* of c[i] *) If[debugnormal, Print["At pt. NORMAL05 in normal, normalizationcoef:"]; Print[normalizationcoef] ]; normalizationcoefcipart = Cases[normalizationcoef,c[_]]; If[debugnormal, Print["At pt. NORMAL06 in normal, normalizationcoefcipart:"]; Print[normalizationcoefcipart] ]; If[ (normalizationcoefcipart =!= {} && Length[normalizationcoefcipart] >= 1), (* then *) normalizationcoef = Coefficient[normalizationcoef,normalizationcoefcipart[[1]]] ]; If[debugnormal, Print["At pt. NORMAL07 in normal, updated normalizationcoef:"]; Print[normalizationcoef] ]; If[debugnormal, Print["At pt. NORMAL08 in normal, divide firsttermrho by firstpartmylist,"<> " normalizationcoef:"]; Print[normalizationcoef] ]; rho[x,t] = Expand[Factor[rho[x,t]/normalizationcoef]]; If[debugnormal, Print["At pt. NORMAL09 in normal, after division by normalizationcoef,"<> " rho[x,t]: "]; Print[rho[x,t]] ]; j[x,t] = Expand[Factor[j[x,t]/normalizationcoef]]; If[debugnormal, Print["At pt. NORMAL10 in normal, after division by normalizationcoef,"<> " j[x,t]: "]; Print[j[x,t]] ]; If[Not[NumberQ[normalizationcoef]], Print["WARNING!"]; Print["Both rho[x,t] and j[x,t] were divided by: ",normalizationcoef,"."] ] (* no else 2, so leave rho and j unscaled *) ]; (* end if 2 *) rhojlist = {rho[x,t],j[x,t]}; If[debugnormal, Print["At pt. NORMAL11 in normal, rhojlist :"]; Print[rhojlist] ]; Return[rhojlist] ]; (* end normal *) (* old normal function *) (* begin taken out (*****************************************************************************) (* 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]]; If[debugnormal, Print["At pt. NORMAL01 in normal, factored form of rho, facexpr: "]; Print[facexpr] ]; denrho = Denominator[facexpr]; If[debugnormal, Print["At pt. NORMAL02 in normal, denominator of rho, denrho: "]; Print[denrho] ]; If[Not[NumberQ[denrho]], denrho = Factor[denrho]]; (* if 1 *) If[Not[FreeQ[Head[denrho],Times]], (* if 2 *) If[NumberQ[Part[denrho,1]] && Part[denrho,1] =!= 0, denrho = denrho/Part[denrho,1]; If[debugnormal, Print["At pt. NORMAL03 in normal, after division by part[denhro,1],"<> " denrho: "]; Print[denrho] ] ] (* end if 2 *) ]; (* end if 1 *) If[Not[NumberQ[denrho]], Print["WARNING!"]; Print["Both rho[x,t] and j[x,t] were multiplied by: ",denrho,"."]; rho[x,t] = Expand[Factor[rho[x,t]*denrho]]; If[debugnormal, Print["At pt. NORMAL04 in normal, after multiplication by denrho,"<> " rho[x,t]: "]; Print[rho[x,t]] ]; j[x,t] = Expand[Factor[j[x,t]*denrho]]]; If[debugnormal, Print["At pt. NORMAL05 in normal, after multiplication by denrho,"<> " j[x,t]: "]; Print[j[x,t]] ]; leadingcoef = Coefficient[rho[x,t],Part[mylist,1]]; If[debugnormal, Print["At pt. NORMAL06 in normal, leadingcoef :"]; Print[leadingcoef] ]; If[NumberQ[leadingcoef] && leadingcoef =!= 0, rho[x,t] = Expand[rho[x,t]/leadingcoef]; If[debugnormal, Print["At pt. NORMAL07 in normal, after division by leadingcoef,"<> " rho[x,t]: "]; Print[rho[x,t]] ]; j[x,t] = Expand[j[x,t]/leadingcoef]; If[debugnormal, Print["At pt. NORMAL08 in normal, after division by leadingcoef,"<> " j[x,t]: "]; Print[j[x,t]] ] ]; rhojlist = {rho[x,t],j[x,t]}; If[debugnormal, Print["At pt. NORMAL09 in normal, rhojlist :"]; Print[rhojlist] ]; Return[rhojlist] ]; (* end normal *) end taken out *) (* 01/23/07 WH renamed variable givenexpr to givenexprfac *) (* 07/11/2006 WH mystripper *) (*********************************************************************** *) (* mystripper[]: cancels the numerical factors and members of nzsymbols, *) (* and anything that has no C in each term of the argument *) (*********************************************************************** *) freeofc[x_]:=FreeQ[x,c]; notfreeofc[x_]:=Not[FreeQ[x,c]]; mystripper[givenexpr_] := Module[ {faclist,factorsfreeofc,factorswithc,factorstokeep,factorstocancel, factorscancelled,givenexprfac}, givenexprfac = Factor[givenexpr]; If[debugmystripper, Print["At pt. MYS1, factored form givenexprfac: ", givenexprfac] ]; If[Head[givenexpr] === Times, (* then *) faclist = FactorList[givenexprfac]; If[debugmystripper, Print["At pt. MYS2, faclist: ", faclist] ]; factorsfreeofc = Select[faclist,freeofc[#]&]; If[debugmystripper, Print["At pt. MYS3, factorsfreeofc: ", factorsfreeofc] ]; factorswithc = Complement[faclist, factorsfreeofc]; If[debugmystripper, Print["At pt. MYS4, factorswithc: ", factorswithc] ]; factorstocancel = Select[factorsfreeofc,(Not[Head[Part[#,1]] === Plus])&]; If[debugmystripper, Print["At pt. MYS5, factorstocancel: ", factorstocancel] ]; factorstokeep = Select[factorsfreeofc,(Head[Part[#,1]] === Plus)&]; If[debugmystripper, Print["At pt. MYS6, factorstokeep: ", factorstokeep] ]; strippedequation = Part[Flatten[factorswithc],1]*Product[ factorstokeep[[k]][[1]]^factorstokeep[[k]][[2]], {k,1,Length[factorstokeep]}]; If[debugmystripper, Print["At pt. MYS7, strippedequation: ", strippedequation] ]; factorscancelled = Factor[givenexpr/strippedequation]; If[debugmystripper, Print["At pt. MYS8, factorscancelled: ", factorscancelled] ], (* else head is not times *) strippedequation = givenexpr ]; (* end if for head testing *) Return[strippedequation] ]; (* end of mystripper *) (*****************************************************************************) (* 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[debugstripper, Print["At pt. STR1, expand givenexpr: ", expr1] ]; If[Head[expr1] === Plus, lenexpr1 = Length[expr1], lenexpr1 = 1 ]; (* start for 1 *) For[k = 1,k <= lenexpr1,k++, If[lenexpr1 == 1, expr2 = expr1, expr2 = Part[expr1,k] ]; expr2 = expr2 /. c[i_] :> 1; If[debugstripper, Print["At pt. STR2, for k=",k,", setting c's to one, expr2: "]; Print[expr2] ]; expr3 = FactorList[expr2]; If[debugstripper, Print["At pt. STR3, for k=",k,", factorlist of expr2, expr3: "]; Print[expr3] ]; lenexpr3 = Length[expr3]; If[debugstripper, Print["At pt. STR4, for k=",k,", lenexpr3: ", lenexpr3] ]; (* start for 2 *) For[s = 1,s <= lenexpr3,s++, expr3s = Part[expr3,s]; If[debugstripper, Print["At pt. STR5, for k=",k," and part s=",s," of expr3, expr3s:"]; Print[expr3s] ]; If[FreeQ[expr3s,u], If[FreeQ[weightpars,Part[expr3s,1]], expr2 = expr2/Part[expr3s,1]^Part[expr3s,2]; If[debugstripper, Print["At pt. STR6, for k=",k,", after division, expr2: ", expr2] ] ], s = lenexpr3+1; ]; ]; (* end for 2 *) list = Append[list,expr2]; If[debugstripper, Print["At pt. STR7, for k=",k,", updated list: "]; Print[list] ] ]; (* end for 1 *) 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 which *) ]; (* 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]; If[debugstripper, Print["At pt. SCALING100, calling stripper in scaling function!"] ]; 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["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["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 *) (*****************************************************************************) (* 07/03/2006 WH added print statements to function subroutine4 *) (* 07/15/2006 WH checked for presence of c[i] first, to avoid redundant work *) subroutine4[givenexpr_] := Module[ {expr1,nexpr1,freerule={},part1,lennexpr1,i}, (* 07/03/2006 WH replaced expand by factor since if was factor in next line *) (* expr1 = Expand[givenexpr]; *) If[debugroutine4, Print["At pt. SUB01, expresssion for rho, givenexpr: "]; Print[givenexpr] ]; expr1 = Factor[givenexpr]; If[debugroutine4, Print["At pt. SUB02, factored from of givenexpr (for rho), expr1: "]; Print[expr1] ]; If[ Cases[expr1,_c,Infinity] === {}, (* then *) freerule={}; If[debugroutine4, Print["At pt. SUB03, no c[i] were present, hence, freerule: "]; Print[freerule] ], (* else continue with computation of freerule *) nexpr1 = Numerator[expr1]; If[debugroutine4, Print["At pt. SUB04, numerator of expr1, nexpr1: "]; Print[nexpr1] ]; lennexpr1 = length[nexpr1]; If[debugroutine4, Print["At pt. SUB05, length of nexpr1, lennexpr1: "]; Print[lennexpr1] ]; For[i = 1,i<= lennexpr1,i++, part1 = Part[nexpr1,i]; If[debugroutine4, Print["At pt. SUB06, first part of nexpr1, part1: "]; Print[part1] ]; If[MemberQ[unknownlist,part1], freerule = Union[freerule,{part1 -> 1}] ]; If[debugroutine4, Print["At pt. SUB07, freerule: "]; Print[freerule] ] ]; (* end for loop *) If[freerule =!= {}, Print["Caution! There is a common (free) factor in the density."]; Print["Setting ",freerule,"."]; ]; If[debugroutine4, Print["At pt. SUB08, final selected, freerule: "]; Print[freerule] ] ]; (* end if *) 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["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}], If[debugstripper, Print["At pt. CONSTRFORMRHO200, calling stripper in constructformrho!"] ]; difflist1 = stripper[D[temppar,{x,inttest-1}]]; lendifflist1 = Length[difflist1]; For[j = 1, j <= lendifflist1,j++, If[debugstripper, Print["At pt. CONSTRFORMRHO300, calling stripper in constructformrho!"] ]; 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["For RANK = ",rhorank,", this is the form of density rho[x,t]: "]; Print[pdeform[Expand[formrho[x,t]]]]; If[formrho[x,t] === 0, Print["The only density is the trivial density rho[x,t] = 0."]; 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,printflag, leadingtermofrho}, 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]]; (* start if 100 *) If[FreeQ[j[x,t],Integrate] && rho[x,t] =!= 0, (* then if 100 *) 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, (* if factortest not zero, then *) Print["Automatic verification of the conservation law FAILED!"], (* else factortest was zero *) If[debugstripper, Print["At pt. EVALUATE400, calling stripper in evaluate function!"] ]; leadingtermofrho = stripper[Part[Expand[rho[x,t]],1]]; (* first application of normal function *) (* WH 08/04/2007 normal function may return c[i] in denominators *) rhojlist = normal[leadingtermofrho]; rho[x,t] = Part[rhojlist,1]; rho[x,t] = Factor[rho[x,t]]; rho[x,t] = Expand[rho[x,t]]; (* 07/14/2006 WH capturing the rho in list *) (* Extra test to verify if the new rho is a multiple of first *) (* computed rho *) (* TO BE DONE *) (* Extra test to verify if the new rho is independent from all previously *) (* computed rhos *) If[listallrhos =!= {}, (* if x *) ratiorhos = Factor[rho[x,t]/listallrhos[[1]]]; If[debugevaluate, Print["At pt. EVAL1, ratiorhos = "]; Print[ratiorhos] ], (* else x *) ratiorhos = u ]; (* end if x *) If[rho[x,t] =!= 0 && Not[FreeQ[ratiorhos,u]], (* then *) (* WH 08/04/2007 to avoid duplicates, added a Union *) listallrhos = Append[listallrhos,rho[x,t]]; listallrhos = Union[listallrhos]; printflag = True, (* else *) printflag = False ]; If[debugevaluate, Print["At pt. EVAL2, listallrhos = "]; Print[listallrhos] ]; j[x,t] = Part[rhojlist,2]; j[x,t] = Factor[j[x,t]]; j[x,t] = Expand[j[x,t]]; (* 07/14/2006 WH capturing the fluxes in list *) If[j[x,t] =!= 0 && Not[FreeQ[ratiorhos,u]], listallfluxes = Append[listallfluxes,j[x,t]]; (* WH 08/04/2007 to avoid duplicates, added a Union *) listallfluxes = Union[listallfluxes]; printflag = True, (* else *) printflag = False ]; If[debugevaluate, Print["At pt. EVAL3, listallfluxes = "]; Print[listallfluxes] ]; If[printflag === True, (* if printflag is true, then *) Print["*******************************************************"]; Print["The normalized density rho[x,t] is:"]; Print[pdeform[rho[x,t]]]; split[pdeform[rho[x,t]]]; If[Length[j[x,t]] <= 20 || debug, Print["*******************************************************"]; Print["The corresponding flux j[x,t] is:"]; Print[pdeform[j[x,t]]], (* else length *) 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["Result of explicit verification: (rho_t + J_x) = ",factortest]; Print["*******************************************************"], (* else when printflag is false *) Print["The rho[x,t] and j[x,t] are linearly dependent on the ones"<> " computed above!"] ], (* end if printflag *) (* else if 100 *) Print["*******************************************************"]; Print["There is NO DENSITY of this form!"]; Print["*******************************************************"] ] (* end if factortest *) ] (* end if 100 *) ]; (* 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, If[debugstripper, Print["At pt. CONSTREQLIST500, calling stripper in constructeqlist!"] ]; 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[] ]; ]; (* end if *) 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."]; 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}]]; ]; (* end for *) Return[simplelist] ]; (* end mysimplify *) (*****************************************************************************) (* newanalyzer[]: determines the coefficients c[i] that must be zero, *) (* Dynamically simplifies system and makes list of c[i] that must be zero *) (* used in the search for compatibility conditions *) (*****************************************************************************) (* analyzelist will be list of two pieces, *) (* first piece: list of simplified lhs of equations *) (* second piece: list of the c[i] that must be zero *) newanalyzer[doublelist_List] := Module[{currentsystem,currentzeroc,strippedsystem,updatedzeroc,analyzelist, updatedzerocrules,updatedsystem}, Print["Starting the repeated analysis and simplification of the system."]; currentsystem = Part[doublelist,1]; If[debugnewanalyze, Print["At pt. ANN01 in newanalyzer, currentsystem: "]; Print[currentsystem] ]; currentzeroc = Part[doublelist,2]; If[debugnewanalyze, Print["At pt. ANN02 in newanalyzer, currentzeroc: "]; Print[currentzeroc] ]; strippedsystem = Union[Map[mystripper,currentsystem]]; If[debugnewanalyze, Print["At pt. ANN03 in newanalyzer, after mapping mystripper,"<> " strippedsystem: "]; Print[strippedsystem] ]; updatedzeroc = Union[Join[currentzeroc,Cases[strippedsystem,_c,1]]]; If[debugnewanalyze, Print["At pt. ANN04 in newanalyzer, updatedzeroc: "]; Print[updatedzeroc] ]; updatedzerocrules = Union[Map[Rule[#,0]&, updatedzeroc]]; If[debugnewanalyze, Print["At pt. ANN05 in newanalyzer, updatedzerocrules: "]; Print[updatedzerocrules] ]; updatedsystem = Complement[strippedsystem /. updatedzerocrules,{0}]; If[debugnewanalyze, Print["At pt. ANN06 in newanalyzer, after application of"<> " updatedzerocrules, updatedsystem: "]; Print[updatedsystem] ]; analyzelist = {updatedsystem, updatedzeroc}; If[debugnewanalyze, Print["At pt. ANN07 in newanalyzer, analyzelist: "]; Print[analyzelist] ]; If[debugnewanalyze, Print["At pt. ANN08 in newanalyzer, "]; Print[" !!!!!!!!! END OF ONE RUN OF NEWANALYZER !!!!!!!!!! "]; ]; Return[analyzelist] ]; (* end newanalyzer *) (*****************************************************************************) (* 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,originalmaineqlist,maineqlist,maineqlistlhs,lengthpar,syscondlist1={}, syscond,inputval,inputrule,sol1,lengthsol1,nrules1,par,newmaineqlist, solc,nrules2,nsolc,funclist,j,seqlist,lenseqlist,analyzelist,analyzeliststart, analyzelistend,listzeroc,listzerocrules,listnonzeroc,myruletrick,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["Applying the Euler operator (the variational derivative)."]; formj = EulerD[ressubst,funclist,{x,t}]; Print["Finished with variational differentiation."]; tformj = formj; (* begin if point 1 *) If[Length[Cases[tformj,x_0]] === noeqs, (* then point 1 *) Print["End of the computations."]; evaluate[formrho[x,t],rhot[x,t],{},rules], (* else point 1 *) dtformj = tformj; dtformj = Expand[dtformj]; Print["Starting the construction of the system for the c[i]."]; 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]; originalmaineqlist = Flatten[Table[picklhs[seqlist,i],{i,1,lenseqlist}]]; system = MapAll[Factor,originalmaineqlist]; If[ (* start if 0 *) (Not[debug] && Length[originalmaineqlist]>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["This is the system for the c[i]:"]; Print[system]; Print["List of unknown coefficients c[i]:"]; Print[unknownlist] ]; (* end if 0 *) (* 02/23/2002 Hereman: workcond.m instead of worklog.m *) (* complexanalysis > 20 replaced by complexanalysis > 0 to force it *) (* to always save the information in a file workcond.m *) lengthpar = Length[parameters]; Print["The unsimplified system for the coefficients c[i] has ", lengthpar," parameter(s)."]; complexanalysis=lengthpar+Length[originalmaineqlist]+Length[unknownlist]; If[complexanalysis > 0, 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 workcond.m for"]; Print["further analysis. For analysis use the Mathematica functions"]; Print["Reduce, Solve, etc. Load the file workcond.m with the command"]; Print["< " run separately!"]; Print["The program will set the nonzero coefficients equal to 1,"<> " one at the time."]; (* remove all c[i] that must be zero *) inputlist = Complement[unknownlist,listzeroc]; If[debugsolve, Print["At pt. SO1, c[i] to be set to one (one at the time), inputlist : "]; Print[inputlist] ]; If[inputlist === {}, (* then 6 no else 6 *) Print["All of the coefficients c[i] in the density have to vanish."]; Print["The search for compatibility conditions ends."]; ]; While[inputlist =!= {}, inputpart = Part[inputlist,1]; inputval = inputpart == 1; inputlist = Complement[inputlist,{inputpart}]; inputrule = ToRules[inputval]; Print["-------------------------------------------------------"]; Print["* Setting ",inputrule[[1]],":"]; Print["Computation of the compatibility conditions."]; If[debugsolve, Print["At pt. SO2a, entering Eliminate after applying inputrule,"<> " maineqlist: "]; Print[maineqlist /. inputrule] ]; If[debugsolve, Print["At pt. SO2b, entering Eliminate, listnonzeroc : "]; Print[listnonzeroc] ]; comcond = Eliminate[maineqlist /. inputrule, listnonzeroc]; If[debugsolve, Print["At pt. SO3a, leaving Eliminate, comcond : "]; Print[comcond] ]; (* 07/03/2006 WH move terms to left hand side of equation, then factor! *) myruletrick = {lhs__ == rhs__ -> lhs - rhs == 0}; comcond = comcond /. myruletrick; If[debugsolve, Print["At pt. SO3b, moving all terms to left hand side, comcond : "]; Print[comcond] ]; comcondfac = MapAll[Factor,comcond]; If[debugsolve, Print["At pt. SO4, after mapall factor on comcond, comcondfac : "]; Print[comcondfac] ]; (* Remove possible duplicates *) If[Head[comcondfac] === Or, comcondfactab = Table[Part[comcondfac,k], {k,1,Length[comcondfac]}]; comcondfac = Union[comcondfactab]; If[debugsolve, Print["At pt. SO5, after removing possible duplicates, comcondfac : "]; Print[comcondfac] ] ]; (* end if head *) Print["This is the compatibility condition:"]; Print[comcondfac]; (* 06/10/2004 LogicalExpand added to make it work under version Math 5.0 *) If[debugsolve, Print["At pt. SO6, entering Reduce, parameters: "]; Print[parameters] ]; If[debugsolve, Print["At pt. SO7, entering Reduce, comcond && syscond: "]; Print[comcond && syscond] ]; If[debugsolve, Print["At pt. SO7bis, entering Reduce, comcond: "]; Print[comcond] ]; (* HERE HERE HERE *) (* Bug fix WH 08/04/2007 different version to make code work under *) (* earlier versions of Mathematica *) (* For Mathematica version 5.0 and higher *) If[$VersionNumber >= 5., flip[equ_] := equ /. { Equal[a_,b_] -> Equal[b,a] }; sol1 = LogicalExpand[Reduce[comcond && syscond, parameters]]; sol2 = LogicalExpand[Reduce[comcond && syscond, Reverse[parameters]]]; sol2 = Map[flip, sol2]; rootfindersol2 = Cases[sol2, Power[_,1/2], Infinity]; If[rootfindersol2 === {}, finalsol = sol2]; rootfindersol1 = Cases[sol1, Power[_,1/2], Infinity]; If[rootfindersol1 === {}, finalsol = sol1]; sol1 = finalsol ]; (* For Mathematica versions lower than 5.0 *) If[$VersionNumber < 5., flip[equ_] := equ /. { Equal[a_,b_] -> Equal[b,a] }; sol1 = Reduce[comcond && syscond, parameters]; sol2 = Reduce[comcond && syscond, Reverse[parameters]]; sol2 = Map[flip, sol2]; rootfindersol2 = Cases[sol2, Power[_,1/2], Infinity]; If[rootfindersol2 === {}, finalsol = sol2]; rootfindersol1 = Cases[sol1, Power[_,1/2], Infinity]; If[rootfindersol1 === {}, finalsol = sol1]; sol1 = finalsol ]; (* END HERE HERE HERE *) If[debugsolve, Print["At pt. SO8, leaving Reduce, sol1: "]; Print[sol1] ]; If[comcondfac === True, Print["The compatibility condition is satisfied without constraints"]; Print["on the parameters."] ]; If[sol1 =!= False, (* then 4 *) If[!FreeQ[sol1,Or], (* then 5 case where sol1 has or in it *) lengthsol1 = Length[sol1]; (* for iii *) For[iii = 1,iii <= lengthsol1,iii++, If[debugsolve, Print["At pt. SO9, taking Part[sol1,iii] : "]; Print[Part[sol1,iii]] ]; nrules1 = ToRules[Part[sol1,iii]]; If[debugsolve, Print["At pt. SO10, nrules1 : "]; Print[nrules1] ]; If[comcondfac=!=True, For[k=1,k <= lengthpar,k++, If[debugsolve, Print["At pt. SO11, before applying nrules1, parameters[[k]]: "]; Print[parameters[[k]]] ]; par[k]=parameters[[k]] /. nrules1; If[debugsolve, Print["At pt. SO12, after applying nrules1, par[k]: "]; Print[par[k]] ]; Print["For ",parameters[[k]]," = ",par[k] ]; ] ]; If[debugsolve, Print["At pt. SO13a, before applying nrules1, maineqlist: "]; Print[maineqlist] ]; newmaineqlist = maineqlist /. nrules1; If[debugsolve, Print["At pt. SO13b, after applying nrules1, newmaineqlist: "]; Print[newmaineqlist] ]; newmaineqlist = newmaineqlist /. inputrule; If[debugsolve, Print["At pt. SO14, after applying inputrule, going into solver"<> " newmaineqlist: "]; Print[newmaineqlist] ]; If[debugsolve, Print["At pt. SO15, unknowns going into solver, listnonzeroc: "]; Print[listnonzeroc] ]; solc = Flatten[Solve[newmaineqlist, listnonzeroc]]; If[debugsolve, Print["At pt. SO16, after solving for c[i], solc: "]; Print[solc] ]; solc = Union[solc,inputrule]; If[debugsolve, Print["At pt. SO17a, after union with inputrule, solc: "]; Print[solc] ]; solc = Union[solc,listzerocrules]; If[debugsolve, Print["At pt. SO17b, after union with listzerocrules, solc: "]; Print[solc] ]; Print["The solution of the system is:"]; Print[solc]; If[debugsolve, Print["At pt. SO18, rules: "]; Print[rules] ]; If[debugsolve, Print["At pt. SO19, nrules1: "]; Print[nrules1] ]; nrules2 = rules /. nrules1; If[debugsolve, Print["At pt. SO20, rules /. nrules1 gives nrules2: "]; Print[nrules2] ]; evaluate[formrho[x,t],rhot[x,t],solc,nrules2]; If[useoldalgorithm, (* if useoldalgorithm then case 1 *) (* begin taken out 1 *) (* case with parameters, and sol1 has or in it, *) (* begin old routine for compatibility analysis *) If[debugsolve, Print["At pt. OLDALGO SO21a, before INTERSECTION,"<> " original inputlist: "]; Print[inputlist] ]; If[debugsolve, Print["At pt. OLDALGO SO21b, solc: "]; Print[solc] ]; If[debugsolve, Print["At pt. OLDALGO SO21c, Position[solc,_->0]: "]; Print[Position[solc,_->0]] ]; If[debugsolve, Print["At pt. OLDALGO SO21d, Flatten[Position[solc,_->0]]: "]; Print[Flatten[Position[solc,_->0]]] ]; If[debugsolve, Print["At pt. OLDALGO SO21e,"<> " solc[[ Flatten[Position[solc,_->0]] ]]:"]; Print[solc[[Flatten[Position[solc,_->0]]]]] ]; If[debugsolve, Print["At pt. OLDALGO SO21f, after rule"<> " x:(c[a_]->0)->c[a] gives:"]; Print[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a] ] ]; If[debugsolve, Print["At pt. OLDALGO SO21g, First[#]& /@ solc :"]; Print[First[#]& /@ solc] ]; If[debugsolve, Print["At pt. OLDALGO SO21h,"<> " Complement[inputlist,First[#]& /@ solc] :"]; Print[Complement[inputlist,First[#]& /@ solc]] ]; If[debugsolve, Print["At pt. OLDALGO SO21i, union of modified solc and"<> " complement :"]; Print[Union[ solc[[ Flatten[Position[solc,_->0]] ]] /. x:(c[a_]->0)->c[a], Complement[inputlist,First[#]& /@ solc] ] ]]; inputlist = Intersection[ inputlist, Union[ solc[[ Flatten[Position[solc,_->0]] ]] /. x:(c[a_]->0)->c[a], Complement[inputlist,First[#]& /@ solc] ] ]; (* end intersection *) If[debugsolve, Print["At pt. OLDALGO SO22, after INTERSECTION,"<> " updated inputlist: "]; Print[inputlist] ], (* else for usenewalgorithm case 1 *) (* case with parameters, and sol1 has or in it, *) (* end old routine for compatibility analysis *) (* end taken out 1 *) (* case with parameters, and sol1 has or in it, *) (* begin new routine for compatibility analysis *) If[inputlist =!= {}, (* then *) inputlist = Drop[inputlist,1]; If[debugsolve, Print["At pt. NEWALGO SO23, case with parameters, after"<> " removing the first element, updated inputlist: "]; Print[inputlist] ], (* else *) inputlist = {} ] (* end if *) ]; (* end if useoldalgorithm case 1 *) (* case with parameters, and sol1 has or in it, *) (* end new routine for compatibility analyis *) Clear[nrules1,nrules2,solc,newmaineqlist]; ], (* for iii *) (* else 5 where sol1 has no or in it *) nrules1 = ToRules[sol1]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO24, nrules1: "]; Print[nrules1] ]; If[comcondfac=!=True, (* then 6 *) For[k = 1,k <= lengthpar,k++, If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO25,"<> " before applying nrules1, parameters[[k]]: "]; Print[parameters[[k]]] ]; par[k]=parameters[[k]] /. nrules1; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO26,"<> " after applying nrules1, par[k]: "]; Print[par[k]] ]; Print["For ",parameters[[k]]," = ", par[k]]; ] ]; (* end if 6 *) If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO27,"<> " before applying nrules1, maineqlist: "]; Print[maineqlist] ]; newmaineqlist = maineqlist /. nrules1; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO28,"<> " after applying nrules1, newmaineqlist: "]; Print[newmaineqlist] ]; newmaineqlist = newmaineqlist /. inputrule; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO29,"<> " after applying inputrule, going into solver"<> " newmaineqlist: "]; Print[newmaineqlist] ]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO30,"<> " unknowns going into solver, listnonzeroc: "]; Print[listnonzeroc] ]; solc = Flatten[Solve[newmaineqlist, listnonzeroc]]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO31,"<> " after solving for c[i], solc: "]; Print[solc] ]; solc = Union[solc,inputrule]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO32a,"<> " after union with inputrule, solc: "]; Print[solc] ]; solc = Union[solc,listzerocrules]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO32b,"<> " after union with listzerocrules, solc: "]; Print[solc] ]; Print["The solution of the system is "]; Print[solc]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO33,"<> " system replacement substitution rules: "]; Print[rules] ]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO34, nrules1: "]; Print[nrules1] ]; nrules2 = rules /. nrules1; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO35,"<> " rules /. nrules1 gives nrules2: "]; Print[nrules2] ]; evaluate[formrho[x,t],rhot[x,t],solc,nrules2]; If[useoldalgorithm, (* if useoldalgorithm then case 2 *) (* begin taken out 2 *) (* case with parameters, and sol1 has no or in it, *) (* begin old routine for compatibility analysis *) If[debugsolve, Print["At pt. OLDALGO SO36a, before INTERSECTION,"<> " original inputlist: "]; Print[inputlist] ]; If[debugsolve, Print["At pt. OLDALGO SO36b, solc: "]; Print[solc] ]; If[debugsolve, Print["At pt. OLDALGO SO36c, Position[solc,_->0]: "]; Print[Position[solc,_->0]] ]; If[debugsolve, Print["At pt. OLDALGO SO36d, Flatten[Position[solc,_->0]]: "]; Print[Flatten[Position[solc,_->0]]] ]; If[debugsolve, Print["At pt. OLDALGO SO36e,"<> " solc[[ Flatten[Position[solc,_->0]] ]]:"]; Print[solc[[Flatten[Position[solc,_->0]]]]] ]; If[debugsolve, Print["At pt. OLDALGO SO36f, after rule"<> " x:(c[a_]->0)->c[a] gives:"]; Print[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a] ] ]; If[debugsolve, Print["At pt. OLDALGO SO36g, First[#]& /@ solc :"]; Print["All c[i] that were fixed after solving:"]; Print[First[#]& /@ solc] ]; If[debugsolve, Print["At pt. OLDALGO SO36h,"<> " Complement[inputlist,First[#]& /@ solc] :"]; Print["All c[i] that remained undetermined after solving:"]; Print[Complement[inputlist,First[#]& /@ solc]] ]; If[debugsolve, Print["At pt. OLDALGO SO36i, union of modified solc and"<> " complement :"]; Print["All c[i] that are undetermined or zero:"]; Print[Union[ solc[[ Flatten[Position[solc,_->0]] ]] /. x:(c[a_]->0)->c[a], Complement[inputlist,First[#]& /@ solc] ] ] ]; inputlist = Intersection[inputlist, Union[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a], Complement[inputlist,First[#]& /@ solc]] ]; If[debugsolve, Print["At pt. OLDALGO SO37, after INTERSECION,"<> " updated inputlist: "]; Print["All c[i] that were zero after solving, plus all "]; Print["undetermined c[i], minus c[i] that must always be zero:"]; Print[inputlist]; Print["Investigate what happens when these c[i] are set to one"] ], (* else if usenewalgorithm case 2 *) (* case with parameters, and sol1 had no or in it, *) (* end old routine for compatibility analysis *) (* end taken out 2 *) (* case with parameters, and sol1 has no or in it, *) (* begin new routine for compatibility analysis *) If[inputlist =!= {}, (* then *) inputlist = Drop[inputlist,1]; If[debugsolve, Print["At pt. NEWALGO SO38, case without parameters,"<> " after removing the first element, updated inputlist: "]; Print[inputlist] ], (* else *) inputlist = {} ] ]; (* end if useoldalgorithm case 2 *) (* case with parameters, and sol1 had no or in it, *) (* end new routine for compatibility analysis *) Clear[nrules1,nrules2,solc,newmaineqlist]; ], (* end if 5 case where sol1 had or in it *) (* else 4 sol1 was false *) Print["The system becomes inconsistent, or the compatibility conditions"]; Print["require that one or more of the parameters are zero."]; Print["Not acceptable!"]; ] (* end if 4 *) ], (* closes While *) (* else 8 start case without parameters *) maineqlist = originalmaineqlist; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO39,"<> " before solving (case without parameters), maineqlist: "]; Print[maineqlist] ]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO40,"<> " before solving (case without parameters), unknownlist: "]; Print[unknownlist] ]; nsolc = Flatten[Solve[maineqlist,unknownlist]]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO41, after solving"<> " (case without parameters), nsolc: "]; Print[nsolc] ]; Print["Solution of the system: "]; Print[nsolc]; evaluate[formrho[x,t],rhot[x,t],nsolc,rules]; Clear[solc1,solc2,solc3,nsolc] ] (* end 8 end if parameter parameter test *) ] (* end point 1 *) ]; (* end module 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[]; (* HERE INITIALLY *) (* 07/14/2006 WH setting up empty lists for listallrhos, listallfluxes *) listallrhos = {}; listallfluxes = {}; If[formrho[x,t] === 0, formrho[x,t] = {}]; If[Not[ListQ[formrho[x,t]]], formrho[x,t] = {formrho[x,t]}]; Print["Working with the data file for the ",name,"."]; For[i = 1,i <= noeqs,i++, eq[i][x,t] = Expand[eq[i][x,t]]; Print["Equation ",i," of the system with ",noeqs," equation(s):"]; 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["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["Program determines the weights of the variables (and parameters)."]; For[i = 1,i <= noeqs,i++, weightu[i] = weightu[i] /. scalerules; If[i === 1, Print["For the given system:"] ]; 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["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["One or more of the weights are negative."]; Print["Negative weights are not allowed."]; Print["Aborting the computations!"]; CloseLog[]; Abort[], If[counterfreeweight =!= 0, (* start if 2 *) If[counterfreeweight > 1, (* start if 3 *) 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["Program will try to determine CHOICES for ", listfreeweights[[1]],":"]; 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["Solving the equation: ",Part[tempvarscalelist,k]," = 1."]; 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[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["List of all the variables (and parameters) with their"]; Print["weights:"]; Print[tvarscalelistval]; 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!"], (* else *) allchoices = Union[allchoices,{solfreeweight}]; ]; ]; ],{k,1,lentempvarscalelist} ]; (* end do loop *) Print["Simple POSITIVE choices are considered."]; 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[Union[Table[Part[listfreeweights,1] -> Part[allchoices,k], {k,1,Length[allchoices]}]]]; 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["Program continues with the choice: ", choicerule]; Print["corresponding to the weights:"]; Print[tvarscalelistchoice] ]; varscalelist = varscalelistchoice ] (* end if 3 *) ] (* end if 2 *) ]; (* end if 1 *) Print["The rank of rho should be an integer multiple of the lowest weight"]; Print["of the DEPENDENT variable(s). Fractional weights are allowed."]; If[ formrho[x,t] === {}, If[Not[NumberQ[rhorank]], rhorank = Input["Enter the rank of rho: "]] ]; If[formrho[x,t] === {}, Print["Computation of the density (and flux) of RANK = ",rhorank]; ]; 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["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[pdeform[formrho[x,t]]] ]; (* 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["*************************************************"] ]; 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]; (* 07/14/2006 WH added several statements to the summary *) Print["********************* BEGIN SUMMARY **********************"]; Print["Total CPU time used in the current session is ",Round[TimeUsed[]], " seconds."]; If[Length[listallrhos] === 1, (* then *) Print["The software computed 1 density with its corresponding flux."], (* else *) Print["The software computed ",Length[listallrhos]," different densities"<> " with their corresponding fluxes."] ]; Print["To see a list of all the densities type: listallrhos," ]; Print["or subscriptform[listallrhos]."]; Print["To see the last density type: rho[x,t], or subscriptform[rho[x,t]]."]; Print["To split the last density in independent pieces type: "]; Print["splitrho[rho[x,t]], or splitrho[subscriptform[rho[x,t]]]."]; Print["To see a list of all the fluxes type: listallfluxes,"]; Print["or subscriptform[listallfluxes]."]; Print["To see the last flux type: j[x,t], or subscriptform[j[x,t]]."]; Print["To split the last flux in independent pieces type: splitflux[j[x,t]],"]; Print["or splitflux[subscriptform[j[x,t]]]."]; Print["************************ END SUMMARY ***************************"]; CloseLog[]; (* WH 08/04/2007 print statement added *) Print["Package condens.m of August 4, 2007 successfully loaded."]; (* ******************************* END ************************************* *)