(* Last updated: October 3, 2007 at 11:25 at Christchurch *) (* Added dynamic stripping as the equations are extracted by the *) (* function constructeqlist. *) (* File: c0828-07.m based on c0811-07.m and c0804-07.m *) (* Code now works under Mathematica versions 4.0, 5.0, 5.2 *) (* Bug fixes see: 08/04/2007 WH and 08/10/2007 WH *) (* Look for TO DO items at top and in code. *) (* TO DO: test under Mathematica version 6.0 *) (* TO DO: should check is a newly computed rho is linearly independent *) (* of all previously computed rhos -- with an algorithm similar to the one *) (* used in the construction of rho with Euler images. *) (* Largely based in file : c060717h.m but updated with pieces of c060718h.m *) (* Previously updated: 07/17/2006 by Hereman at office (Golden) at 15:00 *) (* IMPORTANT: *) (* Code provides the option of using the old algorithm or the new algorithm *) (* to investigate the compatibility conditions. *) (* Flags: UseMinimalAlgorithm = True for the old algorithm and *) (* UseMaximalAlgorithm = True for the new algorithm. *) (* The old algorithm (from Unal's Masters' Thesis) considers a minimal *) (* number of choices for the c[i]. *) (* The new algorithm uses the maximal number of choices for the nonzero *) (* c[i], but it produces reduncancy in the solutions, i.e. some conditions *) (* appear repeatedly. *) (* However, the algorithm comes normalized forms of the densities and *) (* reports whether or not densities are multiples of each other. *) (* It does not perform a complete independence test. *) (* 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 *) (* 02/23/2002 WH *) (*****************************************************************************) (* *) (* *** M A T H E M A T I C A P R O G R A M *** *) (* *) (* SYMBOLIC COMPUTATION OF CONSERVATION LAWS *) (* OF SYSTEMS OF NONLINEAR EVOLUTION EQUATIONS *) (* Output: CONSERVED DENSITIES, FLUXES AND *) (* COMPATIBILITY CONDITIONS *) (* *) (* program name: condens.m *) (* *) (* purpose: computation of the conserved densities and fluxes *) (* 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 coefficients *) (* *) (* u[i]_t = f(u[1],...,u[N],u[1]_{nx},...,u[N]_{nx}) *) (* with i=1,...,N; and n=1,2,... *) (* *) (* algorithm has two options: minimal-algorithm and maximal-algorithm *) (* *) (* output : density and correspoinding flux of desired rank (if it exists), *) (* and compatibility conditions on parameters, if applicable *) (* *) (* tested on : various workstations under Linux; PCs with Windows XP *) (* *) (* language : Mathematica v. 4.0, 5.0, 5.2 and 6.0 *) (* *) (* authors: Unal Goktas and Willy Hereman *) (* Department of Mathematical and Computer Sciences *) (* Colorado School of Mines *) (* Golden, CO 80401-1887, USA *) (* *) (* Version 4: October 3, 2007 (previous version: May 4, 1998) *) (* *) (* Copyright: 1998-2007 *) (* *) (*****************************************************************************) (* 08/04/2007 WH print statement added *) Print["Loading code condens2007.m of October 3, 2007."]; Clear["Global`*"]; (* 08/04/2007 WH and 07/10/2007 WH set flags *) UseMaximalAlgorithm = False; UseMinimalAlgorithm = True; (* 07/03/2006 WH and 08/28/2007 WH added debugflags *) debugstripper = False; debugmystripper = False; debugconstructeqlist = False; debugroutine4 = False; debugnewanalyze = False; debugsolve = False; debugnormal = False; debugevaluate = False; debugmysimplify2 = False; (*****************************************************************************) (* commentinter[]: prints a welcome message *) (*****************************************************************************) (* 08/29/2007 WH separated the comment statement in two pieces *) commentinter1[] := Block[{}, Print["*********************************************************"]; Print[" WELCOME TO THE MATHEMATICA PROGRAM "]; Print[" by UNAL GOKTAS and WILLY HEREMAN "]; Print[" FOR THE COMPUTATION OF CONSERVED DENSITIES, "]; Print[" FLUXES, AND COMPATIBILITY CONDITIONS "]; Print[" Last updated: October 3, 2007 "]; Print[" First released on May 4, 1998 "]; Print[" Copyright 1998-2007 "]; Print["*********************************************************"]; ]; (* end commentinter1 *) commentinter2[] := Block[{}, Print["*********************************************************"]; Print["The list of all parameters is ", allparameters, ". "]; Print["The software assumes that all these parameters are nonzero."]; Print["Cases where parameters are set zero must be run separately."]; Print["*********************************************************"] ]; (* end commentinter2 *) (*****************************************************************************) (* 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 or flux, it 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 expression has no free coefficients."], listofc=Union[Cases[rest,c[__],12]]; lenlistofc=Length[listofc]; Print["There is/are ", lenlistofc, " free coefficient(s)."]; Print["The free coefficients are: ", listofc, "."]; Print["Splitting the expression in independent pieces."]; Do[ ci=Part[listofc,i]; rhoi=Expand[Coefficient[rest,ci]]; rest=Expand[rest-ci*rhoi]; Print["Part of the expression with coefficient ", ci, ":"]; Print[rhoi],{i,1,lenlistofc}]; (* end do *) Print["Part of the expression without 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; (* 08/04/2007 WH 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 *) (* 08/28/2007 WH new version of mystripper *) (*********************************************************************** *) (* mystripper[]: cancels numerical factors, non-zero parameters, and *) (* any factor which 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[givenexprfac] === 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 = Product[Part[factorswithc[[kk]],1],{kk,1,Length[factorswithc]}]* Product[factorstokeep[[k]][[1]]^factorstokeep[[k]][[2]], {k,1,Length[factorstokeep]}]; If[debugmystripper, Print["At pt. MYS7, strippedequation: ", strippedequation] ]; factorscancelled = Factor[givenexprfac/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 *) (* begin taken out 00, old version mystripper (* 01/23/07 WH renamed variable givenexpr to givenexprfac *) (* 07/11/2006 WH mystripper *) (*********************************************************************** *) (* mystripper[]: cancels numerical factors, non-zero parameters, and *) (* any factor which 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 *) end taken out 00, old version 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], (* else *) 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]]; (* 08/28/2007 WH print statement added *) (* Print[newrhot[x,t]]; *) Print["Starting the computation of the flux."]; Print["Starting the integration of D_t(rho) with ", Length[newrhot[x,t]]," terms,"]; 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 *) (* 08/04/2007 WH 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 DO: should check is a newly computed rho is linearly independent *) (* of all previously computed rhos -- with an algorithm similar to the one *) (* used in the construction of rho with Euler images. *) 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 *) (* 08/04/2007 WH 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]]; (* 08/04/2007 WH 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]]]; split[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 of 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 *) (* 08/28/2007 WH added dynamic stripping in the constructeqlist routine *) (*****************************************************************************) (* constructeqlist[]: forms the system of equations for coefficients c[i]. *) (* Input: expression to be unravalled and listzeroc of c[i] that are zero. *) (* Output: double list of equations for c[i] and updated listzeroc. *) (*****************************************************************************) constructeqlist[expr_,listofzeroc_] := Module[ {expr1,rules1,rules2,eqlist = {},tempexpr,newtempexpr,expr2,templist, lentemplist,s,tempcoef,origtempcoef,parttemplist,expr3,listnewzeroc = {}, listzerocrules = {}}, rules1 = D[u[i_][x,t],{x,j_}] :> 0; rules2 = u[i_][x,t] :> 0; expr1 = Expand[expr]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-1, before setting c[i] -> 0, expr1:"]; Print[expr1] ]; (* 08/28/2007 WH produces rules of type c[i] -> 0 *) If[listofzeroc =!= {}, (* start if 0 *) listzerocrules = Union[Map[Rule[#,0]&,listofzeroc]]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-2, setting c[i] -> 0, listzerocrules:"]; Print[listzerocrules] ]; (* 08/28/2007 WH apply rules to expr1 *) expr1 = expr1 /. listzerocrules; expr1 = Expand[expr1]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-3, after setting c[i] -> 0, expr1:"]; Print[expr1] ], (* else *) listzerocrules = {} ]; (* end if 0 *) tempexpr = expr1 //. rules1; tempexpr = tempexpr //. rules2; If[tempexpr =!= 0, If[debugconstructeqlist, Print["At pt. CONSTREQLIST-4, before applying mystripper, tempexpr:"]; Print[tempexpr] ]; newtempexpr = mystripper[tempexpr]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-5, after applying mystripper, newtempexpr:"]; Print[newtempexpr] ]; eqlist = Union[eqlist,{newtempexpr}]; ]; (* removing the constant term from expr2 *) expr2 = expr1-tempexpr; If[expr2 =!= 0, (* start if -1 *) If[debugstripper, Print["At pt. CONSTREQLIST-6, calling stripper in constructeqlist!"] ]; templist = Union[stripper[expr2]]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-7, after applying stripper to expr2, templist:"]; Print[templist] ]; lentemplist = Length[templist]; (* initializing expr3 *) expr3 = expr2; (* initializing the listzeroc *) listzeroc = listofzeroc; For[s = 1,s <= lentemplist,s++, (* start of for 1 *) parttemplist = Part[templist,s]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-8a, working with parttemplist:"]; Print[parttemplist] ]; origtempcoef = Coefficient[Expand[expr2],parttemplist]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-8b, working with origtempcoef:"]; Print[origtempcoef] ]; origtempcoef = origtempcoef //. rules1; origtempcoef = origtempcoef //. rules2; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-8c, after rules1 and rules2, origtempcoef:"]; Print[origtempcoef] ]; (* 08/28/2007 WH apply c[i] -> 0 rules to tempcoef *) If[listzerocrules =!= {}, tempcoef = origtempcoef /. listzerocrules; tempcoef = Factor[tempcoef]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-8d, after setting c[i] -> 0, tempcoef:"]; Print[tempcoef] ], (* else *) tempcoef = origtempcoef ]; If[debug, expr3 = expr3 - (origtempcoef*parttemplist) ]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-9, before applying mystripper, tempcoef:"]; Print[tempcoef] ]; tempcoef=mystripper[tempcoef]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-10, after applying mystripper, tempcoef:"]; Print[tempcoef] ]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-10b, before appending, eqlist:"]; Print[eqlist] ]; eqlist = Append[eqlist,tempcoef]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-10c, after appending tempcoef, eqlist:"]; Print[eqlist] ]; eqlist = Union[eqlist]; (* 08/28/2007 WH pick out the c[i] that must be zero *) newlistzeroc = Cases[eqlist,_c,1]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-10d, after appending, newlistzeroc:"]; Print[newlistzeroc] ]; listzeroc = Join[listzeroc,newlistzeroc]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-10e, joining with newlistzeroc, listzeroc:"]; Print[listzeroc] ]; listzeroc = Union[listzeroc]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-11, pick c[i] that must be zero, listzeroc:"]; Print[listzeroc] ]; (* 08/28/2007 WH produces rules of type c[i] -> 0 *) If[listzeroc =!= {}, (* start if 2 *) listzerocrules = Union[Map[Rule[#,0]&,listzeroc]]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-12, setting c[i] -> 0, listzerocrules:"]; Print[listzerocrules] ]; (* 08/28/2007 WH apply rules to eqlist *) eqlist = eqlist /. listzerocrules; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-13, after setting c[i] -> 0, eqlist:"]; Print[eqlist] ] ]; (* end if 2 *) (* 08/28/2007 WH remove duplicates with Union and also zeros *) eqlist = Union[eqlist]; eqlist = Complement[eqlist,{0}]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-14, after Union and removal of zeros, eqlist:"]; Print[eqlist] ] ]; (* end of for 1 *) 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 -1 *) (* 08/28/2007 WH adding the c[i] that are zero to eqlist *) eqlist = Join[eqlist,listzeroc]; eqlist = Union[eqlist]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-15, after adding zero c[i], eqlist:"]; Print[eqlist] ]; (* 08/28/2007 WH showing listzeroc *) If[debugconstructeqlist, Print["At pt. CONSTREQLIST-16, listzeroc to be returned, listzeroc:"]; Print[listzeroc] ]; Return[{eqlist,listzeroc}] ]; (* end constructeqlist *) (* 08/10/2007 WH simplify2 strips off parameters and powers of parameters *) (* and weigthpars. Also removes duplicates and sorts according to length. *) (*****************************************************************************) (* mysimplify2[]: simplifies the system of equations for the c[i] *) (*****************************************************************************) mysimplify2[list_,parameters_] := Module[ {newlist,lennewlist,k,fterm,simplelist = {}, param}, Print["*********************************************************"]; Print["Starting the second simplification of the system."]; newlist = Complement[list,{0}]; lennewlist = Length[newlist]; For[k = 1 ,k <= lennewlist,k++, fterm = Factor[newlist[[k]]]; If[debugmysimplify2, Print["At pt. MYS2-1, working with fterm:"]; Print[fterm] ]; Do[ param = Part[parameters,p]; If[debugmysimplify2, Print["At pt. MYS2-2, working with param:"]; Print[param] ]; If[Part[fterm,0] === Times && (Cases[fterm,param] =!= {} || Cases[fterm,param^_] =!= {}), If[Cases[fterm,param] =!= {}, If[debugmysimplify2, Print["At pt. MYS2-2, before devision by param, fterm:"]; Print[fterm] ]; fterm = fterm/param; If[debugmysimplify2, Print["At pt. MYS2-3, after devision by param, fterm:"]; Print[fterm] ] ]; If[Cases[fterm,param^_] =!= {}, If[debugmysimplify2, Print["At pt. MYS2-4, before devision by power of param, fterm:"]; Print[fterm] ]; fterm = fterm/(Cases[fterm,param^_][[1]]); If[debugmysimplify2, Print["At pt. MYS2-5, after devision by power of param, fterm:"]; Print[fterm] ] ] ]; (* 08/28/2007 WH changed Union to Append *) If[Intersection[Expand[simplelist],Expand[{fterm*(-1)}]] === {}, simplelist = Append[simplelist,fterm]; simplelist = Union[simplelist] ], {p,1,Length[parameters]} ] (* end do loop *) ]; (* end for loop *) If[debugmysimplify2, Print["At pt. MYS2-6, before applying union to the list, simplelist:"]; Print[simplelist] ]; simplelist = Union[simplelist]; If[debugmysimplify2, Print["At pt. MYS2-7, after union, before sorting the list, simplelist:"]; Print[simplelist] ]; (* sort according to length *) simplelist = Sort[simplelist, (Length[#1] < Length[#2])&]; If[debugmysimplify2, Print["At pt. MYS2-8, after sorting the list, simplelist:"]; Print[simplelist] ]; Return[simplelist] ]; (* end mysimplify2 *) (* 08/10/2007 WH mysimplify strips off numerical factors *) (*****************************************************************************) (* mysimplify[]: simplifies the system of equations for the c[i] *) (*****************************************************************************) mysimplify[list_] := Module[ {newlist,lennewlist,k,fterm,simplelist = {}}, Print["Starting the first 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] ]; (* 08/28/2007 WH changed Union to Append *) If[Intersection[Expand[simplelist],Expand[{fterm*(-1)}]] === {}, simplelist = Append[simplelist,fterm]; simplelist = Union[simplelist] ]; ]; (* 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["*********************************************************"]; Print["Starting the repeated analysis and simplification of the system."]; currentsystem = Part[doublelist,1]; (* 08/28/2007 WH print statement added *) Print["Entering the new analyzer"]; 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 !!!!!!!!!! "]; ]; (* 08/28/2007 WH added print statement *) Print["Leaving the new analyzer"]; 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,dtformjpartj,comcond,k, comcondfactab,myeqlist,originalmaineqlist,maineqlist,maineqlistlhs, lengthpar,syscondlist1={},syscond,inputval,inputrule,sol1,lengthsol1, nrules1,par,newmaineqlist,solc,nrules2,nsolc,funclist,j,seqlist,showseqlist, tempseqlist,tempeqlist,allparameters,lenseqlist,analyzelist,analyzeliststart, analyzelistend,listzeroc,listzerocequations, listzerocrules,listnonzeroc,myruletrick,inputpart, inputlist,i,iii,complexanalysis,comcondfac,equationadjuster}, 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]] ]; (* 08/29/2007 WH changed the message *) 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]; (* 08/29/2007 WH changed the message *) If[debug, Print["The Euler operator will be applied to D_t(rho): ", 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]."]; (* 08/29/2007 WH dynamically apply the simplification of the expressions *) (* from which the system for the c[i] is extracted. *) (* Pass on listzeroc = {} to constructeqlist on first run, then update *) (* listzeroc with all c[i] that are zero, and pass the updated list on *) (* in the next run to the new expression *) (* Input: constructeqlist takes two arguments: the first is the *) (* expression to operate on, the second is listzeroc *) (* Output: constructeqlist returns doublelist {myeqlist,listzeroc} *) myeqlist = {}; listzeroc = {}; For[j = 1,j <= noeqs,j ++, (* begin for 3 *) dtformjpartj = Part[dtformj,j]; If[dtformjpartj =!= 0, (* begin if -3 *) doublelistj = constructeqlist[dtformjpartj,listzeroc]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-X1, returning, doublelistj:"]; Print[doublelistj] ]; myeqlist = Join[myeqlist,Part[doublelistj,1]]; myeqlist = Union[myeqlist]; (* sort according to length *) myeqlist = Sort[myeqlist, (Length[#1] < Length[#2])&]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-X2, updated and resorted myeqlist:"]; Print[myeqlist] ]; listzeroc = Join[listzeroc, Part[doublelistj,2]]; listzeroc = Union[listzeroc]; If[debugconstructeqlist, Print["At pt. CONSTREQLIST-X3, updated listzeroc:"]; Print[listzeroc] ] ] (* end if -3 *) ]; (* end for 3 *) (* 08/10/2007 WH added a second simplification, mysimplify2, to remove *) (* all nonzero parameters and powers of such parameters that appear *) (* as factors. *) (* Sums of parameters are not removed! *) tempseqlist = mysimplify[myeqlist]; Print["*********************************************************"]; Print["This is the system of the c[i] (after the first simplification):"]; tempeqlist = Table[Part[tempseqlist,k] == 0, {k,1,Length[tempseqlist]} ]; Print[tempeqlist]; (* 08/28/2007 WH added print statement *) Print["The simplified system has ",Length[tempeqlist]," equations."]; (* 08/28/2007 WH changed Union to Join *) allparameters = Join[parameters,weightpars]; If[allparameters =!= {}, seqlist = mysimplify2[tempseqlist, allparameters]; (* 08/28/2007 WH added print statement *) Print["This is the system for the c[i] (after second simplification):"]; (* Print[seqlist]; *) showseqlist = Union[Map[Equal[#,0]&,seqlist]]; Print[showseqlist]; (* 08/28/2007 WH added print statement *) Print["The simplified system has ",Length[showseqlist]," equations."], seqlist = tempseqlist ]; lenseqlist = Length[seqlist]; originalmaineqlist = Flatten[Table[picklhs[seqlist,i],{i,1,lenseqlist}]]; system = MapAll[Factor,originalmaineqlist]; (* 08/29/2007 WH begin taken out 00 (* since system is no longer available at the end of the computation *) 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[system], " 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] (after second simplification):"]; Print[system]; (* 08/28/2007 WH added print statement *) Print["The simplified system has ",Length[system]," equations."]; Print["List of unknown coefficients c[i]:"]; Print[unknownlist] ]; (* end if 0 *) end of taken out 00 *) (* 02/23/2002 WH: 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 simplified 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] ]; equationadjuster = {lhs__ == rhs__ -> lhs - rhs == 0}; (* 08/08/2007 WH adjust the equation format of comcond *) (* 07/03/2006 WH move terms to left hand side of equation *) comcond = comcond /. equationadjuster; 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. SO3c, factored compatibility condition, comcondfac:"]; Print[comcondfac] ]; (* 08/08/2007 WH *) (* Reason: write a condition of the form lhs == rhs as lhs - rhs == 0 *) (* and factor the resulting condition *) Print["This is the compatibility condition (with duplicates):"]; Print[comcondfac]; If[debugsolve, Print["At pt. SO4, before removing duplicates, 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 (without duplicates):"]; Print[comcondfac]; (* 06/10/2004 LogicalExpand added to make it work under version Math 5.0 *) If[debugsolve, Print["At pt. SO6, entering Reduce, with parameters: "]; Print[parameters] ]; (* 08/04/2007 WH *) (* TO DO: use mystripper on each piece of comcondfac to strip off *) (* the parameters and their powers that are nonzero! *) (* 08/04/07 WH have the code work with comcondfac instead of comcond *) (* begin taken out 3 If[debugsolve, Print["At pt. SO7, entering Reduce, comcond && syscond: "]; Print[comcond && syscond] ]; end taken out 3 *) (* 06/18/06 WH entering with comcondfac instead of comcond !!! *) If[debugsolve, Print["At pt. SO7, entering Reduce, with comcondfac && syscond: "]; Print[comcondfac && syscond] ]; (* ############### *) If[debugsolve, Print["At pt. SO7bis, entering Reduce, comcondfac: "]; Print[comcondfac]; Print["Code now uses a slightly different solver dependent on the "]; Print["version of Mathematica used. Reason: Reduce works differently "]; Print["under Mathematica v. 4.0 and v. 5.0 !!!"] ]; (* Bug fix 08/04/2007 WH different version to make code work under *) (* versions of Mathematica earlier than v. 5.0 *) (* For Mathematica version 5.0 and higher *) (* 08/04/2007 WH ---> IMPORTANT: work with comcondfac instead of comcond *) If[$VersionNumber >= 5., flip[equ_] := equ /. { Equal[a_,b_] -> Equal[b,a] }; sol1 = LogicalExpand[Reduce[comcondfac && syscond, parameters]]; sol2 = LogicalExpand[Reduce[comcondfac && 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[comcondfac && syscond, parameters]; sol2 = Reduce[comcondfac && 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 ]; 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] ]; (* 08/28/2007 WH changed Union to Join *) solc = Join[solc,inputrule]; solc = Union[solc]; If[debugsolve, Print["At pt. SO17a, after union with inputrule, solc: "]; Print[solc] ]; (* 08/28/2007 WH changed Union to Join *) solc = Join[solc,listzerocrules]; solc = Union[solc]; 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[UseMinimalAlgorithm, (* if UseMinimalAlgorithm 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 UseMaximalAlgorithm (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 UseMinimalAlgorithm (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] ]; (* 08/28/2007 WH changed Union to Join *) solc = Join[solc,inputrule]; solc = Union[solc]; If[debugsolve, Print["At pt. OLDALGO or NEWALGO SO32a,"<> " after union with inputrule, solc: "]; Print[solc] ]; (* 08/28/2007 WH changed Union to Join *) solc = Join[solc,listzerocrules]; solc = Union[solc]; 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[UseMinimalAlgorithm, (* if UseMinimalAlgorithm 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 use-mimimal-algorithm (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 UseMinimalAlgorithm (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["*********************************************************"]; 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]; (* 08/28/2007 WH added list of parameters to commentinter *) (* use Join instead of Union *) allparameters = Join[parameters,weightpars]; commentinter1[]; (* 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]; commentinter2[]; 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["*********************************************************"]; 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}]; (* 08/28/2007 use Join instead of Union *) parameters = Join[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}]; (* 08/29/2007 WH changed the message *) If[debug, Print["*********************************************************"]; Print["Solving the given equation(s) for variables ",mainsolforlist,"."]; 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 pdeform[listallrhos]."]; Print["To see the last density type: rho[x,t], or pdeform[rho[x,t]]."]; Print["To split the last density in independent pieces type: "]; Print["split[rho[x,t]], or split[pdeform[rho[x,t]]]."]; Print["To see a list of all the fluxes type: listallfluxes,"]; Print["or pdeform[listallfluxes]."]; Print["To see the last flux type: j[x,t], or pdeform[j[x,t]]."]; Print["To split the last flux in independent pieces type: split[j[x,t]],"]; Print["or split[pdeform[j[x,t]]]."]; Print["********************** END SUMMARY **********************"]; CloseLog[]; (* 08/04/2007 WH print statement added *) Print["Code condens2007.m of October 3, 2007 successfully loaded."]; (* ******************************* END ************************************* *)