The Homotopy Operator: From Integration by Parts to the Computation of Conservation Laws of Nonlinear PDEs in Multiple Dimensions

Douglas Poole and Willy Hereman

Department of Mathematical and Computer Sciences
Colorado School of Mines
Golden, CO 80401

Homotopy Code

(*****************************************************************************)
(*             ***  M A T H E M A T I C A   P R O G R A M  ***               *)
(*          HOMOTOPY METHOD FOR INVERSION OF DERIVATIVES IN 1-D              *)
(*                   AND DIVERGENCES IN 2-D AND 3-D                          *)
(*                                                                           *)
(*  Program name:HomotopyOperator                                            *)
(*  Purpose: To apply the homotopy method so that a 1-dimensional (1-D)      *)
(*           differential expression can be integrated,or so that a          *)
(*           known divergence in 2-D or 3-D can be inverted.                 *)
(*  Input to homotopyOperator:                                               *)
(*           Any differential expression with one,two or three independent   *)
(*           variables.Independent variables can appear explicitly.          *)
(*           A list of independent variables.                                *)
(*           A list of parameters if any exist in the expression.            *)
(*  Output:  Integrated form of a one-dimensional expression,or an inverted  *)
 (*           divergence with most (if not all) curl terms removed.           *)
(*                                                                           *)
(*  Language: Mathematica 5.0,5.1,5.2,and 6.0                                *)
(*                                                                           *)
(*  Authors: Douglas Poole and Willy Hereman                                 *)
(*  Department of Mathematical and Computer Sciences                         *)
(*  Colorado School of Mines                                                 *)
(*  Golden,CO 80401-1887,U.S.A.                                              *)
(*                                                                           *)
(*  Version 1.1:   24 September,2008                                         *)
(*****************************************************************************)

(*########## Function:variationalDerivativeMultiD ##########*)
(************************************************************************)
(* variationalDerivativeMultiD[expression,dependentVar,indVarList]      *)
(* Purpose: To find the variational derivative of the given expression  *)
(*          in order to determine if the expression in integrable.      *)
(* Input:   An expression                                               *)
(*          A number for the dependent variable                         *)
(* Output:  A partial derivative of the expression with respect         *)
(*          to terms over the jet space                                 *)
(*          The order of the expression with respect to the             *)
(*          independent variable                                        *)
(* Created: 23 May, 2006 by DP                                          *)
(* Last Modified: 10 March, 2008, 15:00 by DP at CSM                    *)
(************************************************************************)variationalDerivativeMultiD[expression_,dependentVar_,indVarList_]:=Module[{lengthIndVars,derivativeList,orderList,higherOrderList,dummyderiv,depVarList,partialWRTpartialList,combinedList,variationalDerivative,termsWithDerList,i,j,indVar1=Part[indVarList,1],indVar2,indVar3,variationalDerivativeList},
(* Assign a standard form to the independent variables. *)
lengthIndVars=Length[indVarList];
If[lengthIndVars>1,indVar2=Part[indVarList,2]];
If[lengthIndVars>2,indVar3=Part[indVarList,3]];
(* Create a list of patial derivatives existing in the expression. *)termsWithDerList=Cases[expression,Derivative[__][u[dependentVar]][Apply[Sequence,indVarList]],{0,Infinity}];
(* Check for terms that have no derivatives. *)
depVarList=Cases[expression,u[dependentVar][Apply[Sequence,indVarList]],{0,Infinity}];
(* Put depVarList and derivative lists together. *)
combinedList=Union[depVarList,termsWithDerList];
(* Create a list of the orders of all derivatives that exist in the expression. *)
higherOrderList=Map[Head,Map[Head,termsWithDerList]]/.Derivative->List;
(* Include order if terms with no derivatives exist. *)orderList=If[depVarList=!={},Union[{Table[0,{i,1,lengthIndVars}]},higherOrderList],Union[higherOrderList],Union[higherOrderList]];
(* Create a list of the partial derivatives of the expression with  respect to the partial derivative list created above. *)partialWRTpartialList=Map[D[expression,#]&,combinedList];
derivativeList=Table[{Part[{indVar1,indVar2,indVar3},i],Part[orderList,j,i]},{j,1,Length[orderList]},{i,1,lengthIndVars}];
(* Compute a list of terms that will add together to form the variational derivative by taking a series of derivatives the previous list. *)
variationalDerivativeList=Map[(-1)^Apply[Plus,#]&,orderList]*ReleaseHold[Thread[dummyderiv[partialWRTpartialList,Map[Hold[Apply[Sequence,#]]&,derivativeList]]]]/.dummyderiv->D;
(* Add the terms in the variational derivative list. *)variationalDerivative=If[Cases[variationalDerivativeList,_?(Head[#]===Sin||Head[#]===Cos||Head[#]===Tan||Head[#]===Sinh||Head[#]===Cosh||Head[#]===Tanh&),{0,Infinity}]=!={},TrigExpand[Apply[Plus,variationalDerivativeList]],Expand[Apply[Plus,variationalDerivativeList]],variationalDerivativeList];
variationalDerivative=If[variationalDerivative=!=0,Factor[variationalDerivative],variationalDerivative,variationalDerivative];
Clear[lengthIndVars,derivativeList,orderList,higherOrderList,dummyderiv,partialWRTpartialList,combinedList,variationalDerivativeList,depVarList,termsWithDerList,i,j,indVar1,indVar2,indVar3];
Return[variationalDerivative]];
(* end Module variationalDerivativeMultiD *)

(*########## Function:homotopyStripper ##########*)
(**********************************************************************)
(* homotopyStripper[givenexpr,indVarsList,paramswithweight]           *)
(* Purpose: Cancels the numerical factors of each term of the         *)
(*          argument and puts these into a list after cancellation    *)
(* Input:   An expression                                             *)
(*          A list of independent variables                           *)
(*          A list of parameters is optional                          *)
(* Output:  A list of the terms of the expression with all constant   *)
(*          coefficients removed.                                     *)
(* Last Modified: 10 March, 2007, 12:31 by DP at CSM                  *)
(**********************************************************************)homotopyStripper[givenexpr_,indVarsList_,paramswithweight_:{}]:=Module[{expnList,expnFnList,coefficientList,singleNoCoefficientList,strippedTermsList,expandExpr,n1,rules1,rules2,rules3},(* If input expression is zero,deal with it separately. *)
If[givenexpr===0,strippedTermsList={0};
Return[strippedTermsList]];
(* Input expression must be in expanded form. *)
expandExpr=Expand[givenexpr];
(* Create a list of terms in the given expression. *)
expnList=If[Head[expandExpr]===Plus,Apply[List,expandExpr],{expandExpr},{expandExpr}];
(* Separate the coefficients and parameters from the terms. *)
expnFnList=Apply[Power,Map[FactorList,expnList],{2,2}];
(* Create a list of coefficients-not including any weighted parameters. Two rules are developed to remove all unwanted parts. *)rules1=_[n1_,___]/;MemberQ[indVarsList,n1]->{};
rules2=Map[Rule[#,{}]&,paramswithweight];
rules3=Map[Rule[#,{}]&,indVarsList];
coefficientList=expnFnList/.rules1;
coefficientList=coefficientList/.rules2;
coefficientList=coefficientList/.rules3;
(* The {} is used first to avoid getting zeros in denominators when rational terms are given in the problem. *)
coefficientList=coefficientList/.{}->1;
(* Combine the coefficients for each term. *)
singleNoCoefficientList=Apply[Times,coefficientList,{1,1}];
(* Form a list of terms without coefficients,but including weighted parameters. *)strippedTermsList=If[Length[expnFnList]>1,Thread[Times[expnList,1/singleNoCoefficientList]],expnList/singleNoCoefficientList];
strippedTermsList=Factor[strippedTermsList];
Clear[expnList,expnFnList,coefficientList,singleNoCoefficientList,expandExpr,n1,rules1,rules2,rules3];
Return[strippedTermsList];];
(*end Module homotopyStripper*)

(*########## Function:applyConstantsToVector ##########*)
(************************************************************************)
(* applyConstantsToVector[vector,lengthIndVar]                          *)
(* Purpose: To apply constants k[i] through the components of a vector  *)
(*          so that the indicies run from 1 to n through all            *)
(*          components.                                                 *)
(* Input:   An expression in rational form                              *)
(*          The number of independent variables                         *)
(* Output:  The numerator of the expression                             *)
(* Created: 3 September,2007 by DP                                      *)
(* Last Modified: 10 March, 2008, 14:56 by DP at CSM                    *)
(************************************************************************)applyConstantsToVector[vector_,lengthIndVar_]:=Module[{i1=1,lengthComp,homotopyResultList={},unknownConstants,homotopyResultPart,termsWithConstants,expnWithConstants},
(* Find the length of each component. *)
lengthComp[1]=If[Head[Part[vector,1]]===Plus,Length[Part[vector,1]],1];
lengthComp[2]=If[lengthIndVar>1,If[Head[Part[vector,2]]===Plus,Length[Part[vector,2]],1],0];
lengthComp[3]=If[lengthIndVar===3,If[Head[Part[vector,3]]===Plus,Length[Part[vector,3]],1],0];
(* Write all terms into a single list. *)
While[i1≤lengthIndVar,homotopyResultPart=If[Head[Part[vector,i1]]===Plus,
Apply[List,Part[vector,i1]],Part[vector,i1]];
homotopyResultList=Flatten[Append[homotopyResultList,homotopyResultPart]];
i1++];(*end While*)
(* Set the list of unknown constants and match the constants with the terms in the homotopy term list and with the ratings for each term. *)unknownConstants=Table[k[i],{i,1,Length[homotopyResultList]}];
(* Attach constants to the terms produced by the homotopy operator. *)
termsWithConstants=Thread[Times[homotopyResultList,unknownConstants]];
(* Remake the homotopy components from the list of homotopy terms with constants so that a divergence can be found.*)expnWithConstants[1]=Apply[Plus,Take[termsWithConstants,lengthComp[1]]];
expnWithConstants[2]=If[lengthIndVar>1,Apply[Plus,Take[termsWithConstants,{lengthComp[1]+1,lengthComp[1]+lengthComp[2]}]],{},{}];
expnWithConstants[3]=If[lengthIndVar===3,Apply[Plus,Take[termsWithConstants,{lengthComp[1]+lengthComp[2]+1,Length[termsWithConstants]}]],{},{}];
Clear[i1,lengthComp,homotopyResultList,homotopyResultPart,termsWithConstants];
Return[{expnWithConstants[1],expnWithConstants[2],expnWithConstants[3],unknownConstants}]];
(* end Module applyConstantsToVector *)

(*########## Function:expressionVarReassign ##########*)
(************************************************************************)
(* expressionVarReassign[expression,indepVariableList]                  *)
(* Purpose: To reassign independent and dependent variables to fit      *)
(*          the format used the all homotopy operator functions.        *)
(* Input:   An expression                                               *)
(*          A list of independent variables                             *)
(* Output:  An expression in the correct format to run in all           *)
(*          functions called by the homotopy operator.                  *)
(* Created: 17 March, 2006 by DP                                        *)
(* Last Modified: 10 March, 2008, 14:43 by DP at CSM                    *)
(************************************************************************)expressionVarReassign[expression_,indepVariableList_]:=Module[{i,depVarList,expressionTerms,n0,n1,n2,n3,n4,integrVarList,newexpression,revisedExpression,depVarAssignRule,positionIntegVars,depVarAssignRuleInvert,lengthIndVarList=Length[indepVariableList],indepCheckList,indList,freeIndVars,newNumDepVars,reverseRule},
(* Modify input expression and variable list to match format used in the homotopy operator functions. *)
(* Input Check 1: Are the independent variables entered into a list in the correct form? *)If[Flatten[indepVariableList]=!=indepVariableList,Print["Independent variables must be entered as a single List, {a,b,c}."];
Abort[]];(* end If Flatten[indepVariableList] *)
(* Input Check 2: Are there more than three independent variables given for integration? *)If[lengthIndVarList>3,Print["This program can integrate expressions with respect to one independent variable, or invert a divergence on two or three independent variables."<>" You have given ",lengthIndVarList,""<>" independent variables, which is beyond the capabilities of this program."];
Abort[]];(* end If lengthIndVarList *)
(* Make lists of dependent variables and independent variables by taking them from the input expression.At the same time check the input for consistency to avoid obscure error messages in later functions. Separate the terms in the given expression,then isolate functions. *)expressionTerms=Flatten[Expand[expression]/.Plus->List/.Times->List/.Power->List//._[n1_[n2__]]->n1[n2]];
expressionTerms=If[Head[expressionTerms]=!=List,{expressionTerms},expressionTerms,expressionTerms];
expressionTerms=Flatten[Map[FactorList,expressionTerms]];
expressionTerms=Union[DeleteCases[expressionTerms,n1_/;NumericQ[n1]||Head[n1]===Symbol]];
(* From expressionTerms,make a list of independent variables.*)indepCheckList=Union[expressionTerms/.{n2_[n1__]/;n2=!=List->List[n1],List[List[n3__]]->List[n3]}];
(* Input Check 3: Are the independent variables put into the functions in a consistent manner? *)If[Length[indepCheckList]>1,Print["The independent variables have not been placed in functions in a consistent manner.  That is, one function may show u[x,y] while"<>" another shows u[y,x] or u[y,z]."];
Abort[]];(* end If Length[indepCheckList] *)
indepCheckList=Flatten[indepCheckList];
(* Input Check 4: Do the independent variables given for integration match the independent variables in the expression? *)If[Not[Part[Union[Map[MemberQ[indepCheckList,#]&,indepVariableList]],1]],Print["Independent variables given for integration do not match the independent variables given in the differential function."];
Abort[]];(* end If Not Union Map *)
freeIndVars=Length[indepCheckList]-lengthIndVarList;
(* From expressionTerms,make a list of dependent variables.*)
indList=Apply[Sequence,indepCheckList];
depVarList=Union[expressionTerms/.n0_[n1_[indList]]/;n0=!=List->n1/.{Derivative[__][n2_][indList]->n2,Derivative[__][n3_]->n3,n4_[indList]->n4}];
(* Input Check 5: Are there any variables appearing in both the independent variable list and the dependent variable list at the same time? *)If[Intersection[indepCheckList,depVarList]=!={},Print["The variable(s), ",Intersection[indepCheckList,depVarList]," occur as both independent and dependent. The program does not know how to interpret this."];
Abort[]];(* end If Complement *)(* Input Check 6: Check that the order of the variables given for integration matches the order they are found in the independent variable list. *)positionIntegVars=Flatten[Map[Position[indepCheckList,#]&,indepVariableList],1];
If[Sort[positionIntegVars]=!=positionIntegVars,integrVarList=Extract[indepCheckList,Sort[positionIntegVars]],
(* else *)
integrVarList=indepVariableList];(* end If Sort *)
(* Build a list of rules reassigning the dependent variables. The expression with reassigned variables will be passed to other functions in the process of running the homotopy operator.*)
depVarAssignRule=Table[Part[depVarList,i]->u[i],{i,1,Length[depVarList]}];
depVarAssignRule=Apply[Rule,DeleteCases[Apply[List,depVarAssignRule,1],{_,_}?(Part[#,1]===Part[#,2]&)],1];
depVarAssignRuleInvert=Table[Reverse[Part[depVarAssignRule,i]],{i,1,Length[depVarAssignRule]}];
newexpression=If[depVarAssignRule=!={},expression/.depVarAssignRule,expression,expression];
(* Rewrite the expression,eliminating any free independent variables reassigning terms with derivatives on the free independent variables,so they do not interfere with the homotopy operator. *)If[freeIndVars>0,{revisedExpression,newNumDepVars,reverseRule}=changeIndependentParameters[newexpression,Length[depVarList],indepCheckList,integrVarList],(*else*)
revisedExpression=newexpression;
newNumDepVars=Length[depVarList];
reverseRule={};];(* end If freeIndVars>0 *)
Clear[i,depVarList,expressionTerms,n0,n1,n2,n3,n4,newexpression,depVarAssignRule,positionIntegVars,lengthIndVarList,indList,freeIndVars];
Return[{revisedExpression,depVarAssignRuleInvert,newNumDepVars,indepCheckList,integrVarList,reverseRule}]];
(* end Module expressionVarReassign *)

(*########## Function:changeIndependentParameters ##########*)
(************************************************************************)
(* changeIndependentParameters[expression,depVarList,indepVarList,      *)
(*      integVarList]                                                   *)
(* Purpose: To rewrite the given expression in terms of the variable(s) *)
(*          required in the integration process.New dependent           *)
(*          variables will replace derivatives on variables not         *)
(*          involved in the integration process.                        *)
(*          e.g.If D[u[1][x,y],x,y] is to be integrated with respect    *)
(*          to x only,the expression will become D[u[2][x],x] with      *)
(*          u[2][x]=D[u[1][x,y],y]                                      *)
(* Input:   An expression                                               *)
(*          A list of dependent variables                               *)
(*          A list of all independent variables found in the given      *)
(*               expression                                             *)
(*          A list of variables involved in integration                 *)
(* Output:  A revised expression.                                       *)
(*          A revised number of dependent variables                     *)
(*          Rules for rewriting the result in the same form as given    *)
(* Created: 10 July, 2007 by DP                                         *)
(* Last Modified: 10 March, 2008, 14:49 by DP at CSM                    *)
(************************************************************************)changeIndependentParameters[expression_,noDepVars_,indepVarList_,integrVarList_]:=Module[{listOfDiffTerms,lengthAllIndVars=Length[indepVarList],positionFreeVars,patternForFreeVars,patternForAllIndVars,newexpression,patternForIntegrVars,patternFreeVarsResult,patternIntegrVarsResult,patternForFreeDerivs,pattern,set,dd1,dd2,m1,m2,j,noOfDerivsListFreeVars,noOfDerivsListIntegrVars,conditions,termsWDerivsOnFreeIndVars,termsWDerivsSeparated,termsWFreeIndVarDers,ruleForFreeIndVarTerms,termsWDerivsOnIntegrVars,ruleToChangeExpn,freeIndVars,numNewDepVars,ruleToReverseChange},
(* Make a list of differential terms in the expression. *)
listOfDiffTerms=Union[Cases[expression,_[___,Part[integrVarList,1],___],{0,Infinity}]];
(* Make a list of free independent variables in the expression. *)
freeIndVars=Complement[indepVarList,integrVarList];
(* Find the positions of the free independent variables in the full independent variable list,then rewrite the list of free variables in the order in which they were found. *)
positionFreeVars=Flatten[Map[Position[indepVarList,#]&,freeIndVars],1];
freeIndVars=Extract[indepVarList,Sort[positionFreeVars]];
(* Create pattern maps for rules so that free variables can be removed and terms containing derivative on free variables can be renamed and kept track of. *)patternForFreeVars=Table[If[MemberQ[positionFreeVars,{j}],pattern[ToExpression["n"<>ToString[j]],Blank[]],Blank[]],{j,1,lengthAllIndVars}];
patternForAllIndVars=Table[pattern[ToExpression["n"<>ToString[j]],Blank[]],{j,1,lengthAllIndVars}];
patternForAllIndVars=patternForAllIndVars/.pattern->Pattern;
patternForIntegrVars=Table[If[!MemberQ[positionFreeVars,{j}],pattern[ToExpression["n"<>ToString[j]],Blank[]],Blank[]],{j,1,lengthAllIndVars}];
patternFreeVarsResult=DeleteCases[patternForFreeVars/.pattern->Sequence,Verbatim[Blank[]]];
patternIntegrVarsResult=DeleteCases[patternForIntegrVars/.pattern->Sequence,Verbatim[Blank[]]];
noOfDerivsListFreeVars=Thread[List[freeIndVars,patternFreeVarsResult]];
noOfDerivsListIntegrVars=Thread[List[integrVarList,patternIntegrVarsResult]];
patternForFreeDerivs=patternForFreeVars/.pattern->Pattern;
patternForIntegrVars=patternForIntegrVars/.pattern->Pattern;
conditions=Apply[Or,Map[Greater[#,0]&,patternFreeVarsResult]];
termsWDerivsOnFreeIndVars=Union[Cases[listOfDiffTerms,Derivative[Apply[Sequence,patternForFreeDerivs]][u[_]][__]/;Evaluate[conditions]]];
(* Separate the derivatives on free independent variables from the part with derivatives on integrating variables. *)termsWDerivsSeparated=termsWDerivsOnFreeIndVars/.Derivative[Apply[Sequence,patternForAllIndVars]][u[m1_]][m2__]->dd1[dd2[u[m1][m2],Apply[Sequence,noOfDerivsListFreeVars]],Apply[Sequence,noOfDerivsListIntegrVars]]/.dd2->D;
termsWFreeIndVarDers=Union[Map[Part[#,1]&,termsWDerivsSeparated]];
(* Reassign the parts with derivatives on free variables to new dependent variables. *)
numNewDepVars=Length[termsWFreeIndVarDers];
ruleForFreeIndVarTerms=Table[set[u[j][Apply[Sequence,indepVarList]],Part[termsWFreeIndVarDers,j-noDepVars]],{j,noDepVars+1,numNewDepVars+noDepVars}];
ruleForFreeIndVarTerms=ruleForFreeIndVarTerms/.set[m1_,m2_]->Rule[m2,m1];
termsWDerivsOnIntegrVars=termsWDerivsSeparated/.ruleForFreeIndVarTerms/.dd1->D;
ruleToChangeExpn=Thread[Rule[termsWDerivsOnFreeIndVars,termsWDerivsOnIntegrVars]];
newexpression=expression/.ruleToChangeExpn;
newexpression=newexpression/.{u[m1_][Apply[Sequence,indepVarList]]->u[m1][Apply[Sequence,integrVarList]],Derivative[Apply[Sequence,patternForIntegrVars]][u[m1_]][Apply[Sequence,indepVarList]]->Derivative[Apply[Sequence,patternIntegrVarsResult]][u[m1]][Apply[Sequence,integrVarList]]};
ruleToReverseChange=If[ruleForFreeIndVarTerms=!={},ruleForFreeIndVarTerms/.Rule[m1_,m2_]->Rule[m2,m1],{{}}];
Clear[listOfDiffTerms,lengthAllIndVars,positionFreeVars,patternForFreeVars,patternForAllIndVars,patternForIntegrVars,patternFreeVarsResult,patternIntegrVarsResult,patternForFreeDerivs,pattern,set,dd1,dd2,m1,m2,j,noOfDerivsListFreeVars,noOfDerivsListIntegrVars,conditions,termsWDerivsOnFreeIndVars,termsWDerivsSeparated,termsWFreeIndVarDers,ruleForFreeIndVarTerms,termsWDerivsOnIntegrVars,ruleToChangeExpn,freeIndVars];
Return[{newexpression,numNewDepVars+noDepVars,ruleToReverseChange}]];
(* end Module changeIndependentParameters *)

(*########## Function:returnAuxIndependentVars ##########*)
(************************************************************************)
(* returnAuxIndependentVars[homotopyResult,fullIndVarList,              *)
(*     integrateVarList,reverseRules]                                   *)
(* Purpose: To replace renamed dependent variables with the derivative  *)
(*          on auxillary variables removed in                           *)
(*          removeFreeIndpendentVariables,and reinsert the auxillary    *)
(*          variables back into the expression,returning it to its      *)
(*          original form.                                              *)
(* Input:   The result returned by the homotopy operator                *)
(*          A list of all independent variables original expression     *)
(*          A list of rules for reversing the changes made in the       *)
(*              function removeFreeIndpendentVariables                  *)
(* Output:  A revised expression.                                       *)
(* Created: 12 July, 2007 by DP                                         *)
(* Last Modified: 10 March, 2008, 13:04 by DP at CSM                    *)
(************************************************************************)returnAuxIndependentVars[homotopyResult_,fullIndVarList_,integrateVarList_,reverseRules_]:=Module[{positionIntegVars,patternForIntegVars,initialPattern,finalPattern,newExpression,addDepVarNos,auxDepTerms,auxDepTermsRev,ruleAuxToOldForm,pattern,dd1,dd2,j,m1,m2,m3,thread,list},
(* Put all auxillary independent variables back into the homotopy result. *)
positionIntegVars=Flatten[Map[Position[fullIndVarList,#]&,integrateVarList],1];
patternForIntegVars=Table[If[MemberQ[positionIntegVars,{j}],pattern[ToExpression["n"<>ToString[j]],Blank[]],0],{j,1,Length[fullIndVarList]}];
initialPattern=DeleteCases[patternForIntegVars,0]/.pattern->Pattern;
finalPattern=DeleteCases[patternForIntegVars/.pattern->Sequence,Verbatim[Blank[]]];
newExpression=homotopyResult/.{Derivative[Apply[Sequence,initialPattern]][u[m1_]][__]->Derivative[Apply[Sequence,finalPattern]][u[m1]][Apply[Sequence,fullIndVarList]],u[m2_][__]->u[m2][Apply[Sequence,fullIndVarList]]};
(* Apply reverse rules to rename terms with derivatives on auxillary variables. *)If[reverseRules=!={{}},addDepVarNos=Map[Part[#,1]&,reverseRules]/.u[m1_][__]->m1;
auxDepTerms=Cases[newExpression,Derivative[__][u[m1_]][__]/;MemberQ[addDepVarNos,m1],{0,Infinity}];
(* Restructure the auxillary dependent variable terms so that their derivatives are included in the original form.*)auxDepTermsRev=auxDepTerms/.Derivative[m1__][m2_][m3__]->dd2[m2[m3],thread[list[{m3},{m1}]]]/.thread->Thread/.list->List;
auxDepTermsRev=auxDepTermsRev/.reverseRules;
auxDepTermsRev=ReleaseHold[auxDepTermsRev/.dd2[m1__,m2_]->dd1[m1,Hold[Apply[Sequence,m2]]]]/.dd1->D;
auxDepTerms=auxDepTerms/.reverseRules;
ruleAuxToOldForm=Thread[Rule[auxDepTerms,auxDepTermsRev]];
newExpression=newExpression/.ruleAuxToOldForm/.reverseRules];(* end If reverseRules *)Clear[positionIntegVars,patternForIntegVars,initialPattern,finalPattern,addDepVarNos,auxDepTerms,auxDepTermsRev,ruleAuxToOldForm,pattern,dd1,dd2,j,m1,m2,m3,thread,list];
Return[newExpression]];
(* end Module returnAuxIndependentVars *)

(*########## Function:integrabilityTest ##########*)
(************************************************************************)
(* integrabilityTest1D[expression,dependentVar,listOfIndepVars,         *)
(*     (optional) parameters]                                           *)
(* Purpose: To determine if the expression is integrable by applying    *)
(*          the zero-Euler operator.If the expression is not            *)
(*          integrable and parameters have been given,then the          *)
(*          function will try to identify values for the parameters     *)
(*          to make the expression integrable.                          *)
(* Input:   An expression                                               *)
(*          A number for the dependent variable                         *)
(*          A list of independent variables                             *)
(*          A list of parameters (optional)                             *)
(* Output:  A partial derivative of the expression with respect to      *)
(*              terms over the jet space                                *)
(*          The order of the expression with respect to the             *)
(*              independent variable                                    *)
(* Created: 22 March, 2006 by DP                                        *)
(* Last Modified: 10 March, 2008 at 14:46 by DP                         *)
(************************************************************************)integrabilityTest[expression_,noOfDepVar_,listOfIndepVars_,parameters_:{}]:=Module[{j,variationalDerivativeResult,testBasis,parameterSolutions,parameterList,parameterRule,intelligibleAnswer1=0,lthiv,intelligibleAnswer2=0,answer1,answer2,secondTest,lps,lpm,lpa},
(* Check parameter list for proper format. *)
parameterList=If[Length[parameters]===0&&parameters=!={},{parameters},parameters,parameters];
(* Find the zero-Euler operator for the expression for each component of the dependent variable. *)variationalDerivativeResult=Table[variationalDerivativeMultiD[expression,j,listOfIndepVars],{j,1,noOfDepVar}];
testBasis=Table[0,{j,1,noOfDepVar}];
If[variationalDerivativeResult===testBasis,If[globalHomotopyVerbose,lthiv=Length[listOfIndepVars];
printOUT["The zero-Euler operator applied to the given expresion has produced zero for all dependent variables.  The expression is a "<>If[lthiv===1,"total derivative.","divergence."]<>"  The "<>If[lthiv===1,"total derivative","divergence"]<>" will be "<>If[lthiv===1,"integrated","inverted"]<>" with respect to "<>If[lthiv===1,Apply[ToString,listOfIndepVars],StringReplace[ToString[Insert[listOfIndepVars,"and",lthiv]],{"{"->"","}"->"",", and,"->" and"}]]<>"."];
printOUT["------------------------------------------"];];(* end globalHomotopyVerbose *)
Return[{"Integrable",{}}]];(* end If variationalDerivativeResult *)If[variationalDerivativeResult=!=testBasis&&parameterList==={},printOUT["The zero-Euler operator applied to the expression has not  produced zero for all dependent variables.  The expression is"<>" not integrable. The integration process will not be activated."];
Clear[j,variationalDerivativeResult,testBasis,parameterSolutions,parameterList,parameterRule,intelligibleAnswer1,printIT,intelligibleAnswer2,answer1,answer2,lps,lpm,lpa];
Return[{"Not Integrable",{}}]];(* end If variationalDerivativeResult *)If[variationalDerivativeResult=!=testBasis&&parameterList=!={},parameterSolutions=Flatten[Table[parameterTest[Part[variationalDerivativeResult,j],j,parameterList],{j,1,noOfDepVar}]];
lps=If[Length[parameterSolutions]>1,"values","a value"];
lpm=If[Length[parameterList]>1,"parameters","parameter"];
lpa=If[Length[parameterList]>1,"are","is"];
Print["The zero-Euler operator applied to the expression has not produced zero for all dependent variables.  The expression may"<>" become integrable if the "<>lpm<>" "<>lpa<>" set to the value(s) given in the following list"];
Print[parameterSolutions];
While[intelligibleAnswer1===0,answer1=ToString[Input["See the information about parameters on the main screen.  The program suggests "<>lps<>" for the "<>lpm<>" given in the list, "<>ToString[parameterList]<>", which may allow the integration process to continue.  Do you  want to set "<>lps<>" for the "<>lpm<>" and continue with the problem?  (Yes or No)"]];
If[MemberQ[{"No","N","no","n","NO"},answer1],intelligibleAnswer1=1;
Print["The expression is not integrable. The integration process"," will not be activated."];
Clear[j,variationalDerivativeResult,testBasis,parameterSolutions,parameterList,parameterRule,intelligibleAnswer1,printIT,intelligibleAnswer2,answer1,answer2,lps,lpm,lpa];
Return[{"Not Integrable",{}}]];(* end If answer1===No *)
If[MemberQ[{"Yes","Y","yes","y","YES"},answer1],intelligibleAnswer1=1;
answer2=Input["Enter a value or a list of values for the"<>" parameter list,   "<>ToString[parameterList]<>".         "<>"\nIf your parameter list is {alpha, beta, gamma} then enter a list of values you wish the parameters to have, e.g., {-2,"<>" Sqrt[3], gamma}.  If you do not wish to change a parameter or have a free parameter as suggested on the screen, re-enter the"<>" name of the parameter in the list."];
While[intelligibleAnswer2===0,If[Length[answer2]=!=Length[parameterList]||Head[answer2]=!=List,answer2=Input["Please re-enter your list.  The list must have "<>ToString[Length[parameterList]]<>" item(s) contained in braces, {}."],intelligibleAnswer2=1] (* end If Length *)] (* end While intelligibleAnswer2 *)] (* end If answer1===Yes *)];(* end While intelligilbeAnswer1 *)printOUT["The values given are  ",TableForm[{parameterList,"=",answer2},TableDirections->{Row,Column}]];
printOUT["After a recheck for integrability, the integration process will begin."];
parameterRule=Thread[Rule[parameterList,answer2]];
secondTest=integrabilityTest[expression/.parameterRule,noOfDepVar,listOfIndepVars]];(*end If variationalDerivativeResult*)Part[secondTest,2]=parameterRule;
Clear[j,variationalDerivativeResult,testBasis,parameterSolutions,parameterList,parameterRule,intelligibleAnswer1,intelligibleAnswer2,answer1,answer2,lps,lpm,lpa];
Return[secondTest]] 
(* end Module integrabilityTest *)

(*########## Function:parameterTest ##########*)
(************************************************************************)
(* parameterTest[variationalDerivativeResult,dependentVar,parameters]   *)
(* Purpose: To try to identify values for the parameters that will      *)
(*          make the expression integrable.                             *)
(* Input:   An expression                                               *)
(*          A number for the dependent variable                         *)
(*          A list of parameters                                        *)
(* Output:  Suggested values for parameters,or a statement that no      *)
(*              such parameter exists.                                  *)
(* Created: 22 March, 2006 by DP                                        *)
(* Last Modified: 10 March, 2008, 15:09 by DP at CSM                    *)
(************************************************************************)

parameterTest[variationalDerivative_,dependentVar_,parameterList_]:=Module[{factoredResult,factoredResultList,parameterSolutions1,parameterSolutions2={},parameterSolutions3={},parameterSolutions,paramindex=1,i,resultSolutions,listOfParameterTerms,listOfDerivativeTerms,possibleNewFactors,variationalDerivativeResult},(* Create a list of terms by applying Factor to the variational derivative result. *)
variationalDerivativeResult=Expand[variationalDerivative];
factoredResult=Factor[variationalDerivativeResult];
factoredResultList=Apply[List,factoredResult];
(* To further try to isolate values for parameters, set each list member equal to zero and solve.  A value for a parameter may result. *)resultSolutions=DeleteCases[Map[Resolve[#==0,Map[#&,parameterList]]&,factoredResultList],False];
(* Strip the parameter solutions off of the previous result. *)
parameterSolutions1=Flatten[Map[Cases[resultSolutions,Equal[#,___],Infinity]&,parameterList]];
parameterSolutions1=DeleteCases[parameterSolutions1,Equal[_,___?(MemberQ[#,Derivative[___][u[_]][___],Infinity]===True||MemberQ[#,u[_][___],Infinity]===True&)],Infinity];
(* Break down and isolate parameters in expression terms,then match terms in the result from the variational derivative to create additional factors. *)While[paramindex≤Length[parameterList],listOfParameterTerms=Cases[variationalDerivativeResult,_?(MemberQ[#,Part[parameterList,paramindex],Infinity]===True&)];
listOfDerivativeTerms=Union[Flatten[Table[Flatten[{Cases[listOfParameterTerms[[i]],Derivative[___][u[_]][___],Infinity],Cases[listOfParameterTerms[[i]],u[_][___],Infinity]}],{i,1,Length[listOfParameterTerms]}]]];
possibleNewFactors=Union[Table[Cases[variationalDerivativeResult,_?(MemberQ[#,listOfDerivativeTerms[[i]],Infinity]===True&)],{i,1,Length[listOfDerivativeTerms]}]];
parameterSolutions2=Append[parameterSolutions2,Solve[Apply[Plus,possibleNewFactors,1]==0,Part[parameterList,paramindex]]/.Rule->Equal];
parameterSolutions2=DeleteCases[Flatten[parameterSolutions2],Equal[_,___?(MemberQ[#,Derivative[___][u[_]][___],Infinity]===True||MemberQ[#,u[_][___],Infinity]===True&)],Infinity];
(* Directly solve for parameters using the result from the variational derivative. *)parameterSolutions3=Append[parameterSolutions3,Solve[variationalDerivativeResult==0,Part[parameterList,paramindex]]/.Rule->Equal];
parameterSolutions3=DeleteCases[Flatten[parameterSolutions3],Equal[_,___?(MemberQ[#1,Derivative[___][u[_]][___],Infinity]===True||MemberQ[#,u[_][___],Infinity]===True&)],Infinity];
paramindex++];(* end While *)
parameterSolutions=Union[parameterSolutions1,parameterSolutions2,parameterSolutions3];
Clear[factoredResult,factoredResultList,parameterSolutions1,parameterSolutions2,parameterSolutions3,paramindex,i,resultSolutions,listOfParameterTerms,listOfDerivativeTerms,possibleNewFactors,variationalDerivativeResult];
Return[parameterSolutions]];
(* end Module parameterTest *)
(*########## Function:degreeOperator ##########*)
(************************************************************************)
(* degreeOperator[expression,dependentVar,IndepVars]                    *)
(* Purpose: To determine the overall degree of a differential function  *)
(* Input:   An expression                                               *)
(*          The number of dependent variables                           *)
(*          A list of independent variables                             *)
(* Output:  A revised expression with information about its degree      *)
(*              contained in its coefficient                            *)
(* Created: 11 June, 2007 by DP                                         *)
(* Last Modified: 24 March, 2008, 17:28 by DP at CSM                    *)
(************************************************************************)degreeOperator[term_,numDependentVar_Integer,IndepVars_List]:=Module[{numIndVar,termsWithDerList,depVarList,combinedList,partialWRTpartialList,counter=1,finalDegree={}},numIndVar=Length[IndepVars];
(* Add results for all dependent variables. *)
While[counter≤numDependentVar,
(* Create a list of patial derivatives existing in the expression. *)
termsWithDerList=Cases[term,Derivative[__][u[counter]][__],{0,Infinity}];
(* Check for terms that have no derivatives. *)
depVarList=Cases[term,u[counter][__],{0,Infinity}];
(* Put depVarList and derivative lists together. *)
combinedList=Union[depVarList,termsWithDerList];
(* Create a list of the partial derivatives of the expression with respect to the partial derivative list created above. *)partialWRTpartialList=Map[D[term,#]&,combinedList];
(* Multiply the partial term to the partial derivative just found. *)
partialWRTpartialList=Thread[Times[partialWRTpartialList,combinedList]];
partialWRTpartialList=Expand[partialWRTpartialList];
finalDegree=Flatten[Append[finalDegree,partialWRTpartialList]];
counter++];(* end While *)
finalDegree=Factor[Apply[Plus,finalDegree]];
Clear[numIndVar,termsWithDerList,depVarList,combinedList,partialWRTpartialList,counter];
Return[finalDegree]];
(* end Module degreeOperator *)

degreeOperExpression[expression_,numDepVars_Integer,indVarsList_List]:=Module[{expandedExpression,degreeOperResult},
expandedExpression=Expand[expression];
expandedExpression=If[Head[expandedExpression]===Plus,Apply[List,expandedExpression],{expression},{expression}];
degreeOperResult=Map[degreeOperator[#,numDepVars,indVarsList]&,expandedExpression];
Return[{expandedExpression,degreeOperResult}]];
(* end Module degreeOperExpression *)

(*########## Function:rationalExpressionCheck ##########*)
(************************************************************************)
(* rationalExpressionCheck[expressionGiven,indepVars,numDepVars]        *)
(* Purpose: To check an expression for terms that fall in the kernel    *)
(*          of the degree operator,then to separate those terms from    *)
(*          the rest of the expression.The separate pieces will be      *)
(*          run through the homotopy operator independently             *)
(* Input:   An expression                                               *)
(*          A list of terms in the denominator                          *)
(* Output:  The combined result from the homotopy operator              *)
(* Created: 27 March, 2008 by DP                                        *)
(* Last Modified: 1 April, 2008, 18:03 by DP at CSM                     *)
(************************************************************************)rationalExpressionCheck[expressionGiven_,indepVars_,numDepVars_]:=Module[{expnDenom,revisedExpn,degree,zeroDegreeTerms,shiftedTerms,afterHomotopy,degreeChk},
(* Check the expression for an overall degree of zero.If the degree is zero,the homotopy operator will not work until changes are made. *)expnDenom=MemberQ[expressionGiven,Power[n1__,n2_]/;MemberQ[n1,u[_][__]|Derivative[__][u[_]][__],{0,Infinity}]&&n2<0,{0,Infinity}];
If[expnDenom,{revisedExpn,degree}=degreeOperExpression[expressionGiven,numDepVars,indepVars];
(* If the expression has a term with a degree of zero,shift one of the terms in the denominator of the degree zero term. *)If[MemberQ[degree,0],zeroDegreeTerms=Extract[revisedExpn,Position[degree,0,{0,1}]];
shiftedTerms=shiftDenominatorTerm[Apply[Plus,zeroDegreeTerms],indepVars,numDepVars];
revisedExpn={Apply[Plus,Complement[revisedExpn,zeroDegreeTerms]],shiftedTerms};
(* Recheck the degree before moving on.*)
degreeChk=degreeOperExpression[Part[revisedExpn,2],numDepVars,indepVars];
(* If degree sends 0 to homotopyFinalIntegration, Mathematica's Integrate will be called for 1D functions. *)
degreeChk=If[MemberQ[Part[degreeChk,2],0],0,1];
revisedExpn=If[degreeChk===0,{expressionGiven},revisedExpn];
(* Run the homotopy operator independently on the terms not in the kernel of the degree operator, the on the terms in the kernel of the degree operator. *)
afterHomotopy=Map[homotopyFinalIntegration[#,numDepVars,indepVars,degreeChk]&,revisedExpn];
(* Undo any shifts applied to terms in the kernel of the degree operator. *)
afterHomotopy=afterHomotopy/.Log[shift]->0/.shift->0;
afterHomotopy=MapThread[Plus,afterHomotopy,1],
(* else *)
degreeChk=1;
afterHomotopy=homotopyFinalIntegration[expressionGiven,numDepVars,indepVars,degreeChk]],(* end If MemberQ[degree] *)
(* else *)
degreeChk=1;
afterHomotopy=homotopyFinalIntegration[expressionGiven,numDepVars,indepVars,degreeChk]];(* end If expnDenom *)
(* Clear all variables not being passed back. *)
Clear[expnDenom,revisedExpn,degree,zeroDegreeTerms,shiftedTerms,degreeChk];
Return[afterHomotopy]]
 (* end Module rationalExpressionCheck *)

(*########## Function:shiftDenominatorTerm ##########*)
(************************************************************************)
(* shiftDenominatorTerm[expressionGiven,indepVars,numDepVars]           *)
(* Purpose: To shift a term in the denominator of a rational expression *)
(*          away from the origin. This alleviates the homogenous        *)
(*          degree problem.                                             *)
(* Input:   An expression                                               *)
(*          A list of terms in the denominator                          *)
(* Output:  A list of degrees corresponding to each term in the         *)
(*              expression                                              *)
(* Created: 19 March, 2008 by DP                                        *)
(* Last Modified: 22 May, 2008, 10:29 by DP at CSM                      *)
(************************************************************************)

(* The function tuples replaces Mathematics's Tuples which is not available in versions before 5.1. *)
tuples[listOflists_List]:=Module[{end=Length[listOflists],tableOfParts,e,i1,j1,tupledList,part},tableOfParts={Table[part[listOflists,j1,e[j1]],{j1,1,end}],Apply[Sequence,Table[{e[i1],1,Length[Part[listOflists,i1]]},{i1,1,end}]]};
tupledList=Flatten[Apply[Table,tableOfParts]/.part->Part,end-1];
Return[tupledList]]

shiftDenominatorTerm[expressionGiven_,indepVars_,numDepVars_]:=Module[{expression,zero,denominatorList,denominatorTermList,orderList,commonTerms,varDerResult=1,termsToShift,shiftedExpression,pls,shiftedTerms},(* Initializations *)
Clear[shift];
zero=Table[0,{i1,1,numDepVars}];
expression=Expand[expressionGiven];
(* Get the denominators of the terms with a degree of zero. *)
denominatorList=Map[Denominator,Apply[List,expression]];
denominatorTermList=Union[Map[Cases[#,u[_][__]|Derivative[__][u[_]][__],{0,Infinity}]&,denominatorList]];
denominatorTermList=DeleteCases[denominatorTermList,{}];
(* Sort the terms in the denominator by order. *)
orderList=denominatorTermList/.{u[_][__]->0,Derivative[n1__][u[_]][__]->pls[n1]}/.pls->Plus;
denominatorTermList=Sort[Thread[List[orderList,denominatorTermList]]];
denominatorTermList=Map[Part[#,2]&,denominatorTermList];
commonTerms=Apply[Intersection,denominatorTermList];
denominatorTermList=DeleteCases[denominatorTermList,Apply[Alternatives,commonTerms],{0,3}];
While[commonTerms=!={}&&varDerResult=!=zero,termsToShift=First[commonTerms];
shiftedTerms=Rule[termsToShift,termsToShift+shift];
shiftedExpression=expression/.shiftedTerms;
varDerResult=Factor[Table[variationalDerivativeMultiD[shiftedExpression,i1,indepVars],{i1,1,numDepVars}]];
commonTerms=If[commonTerms=!={},Drop[commonTerms,1],commonTerms];];(* end While *)denominatorTermList=If[varDerResult===zero,denominatorTermList,tuples[denominatorTermList]];
While[Flatten[denominatorTermList]=!={}&&varDerResult=!=zero,termsToShift=First[denominatorTermList];
shiftedTerms=Map[Rule[#,#+shift]&,termsToShift];
shiftedExpression=expression/.shiftedTerms;
varDerResult=Factor[Expand[Table[variationalDerivativeMultiD[shiftedExpression,i1,indepVars],{i1,1,numDepVars}]]];
denominatorTermList=If[denominatorTermList=!={},Drop[denominatorTermList,1],denominatorTermList];];(*end While*)shiftedExpression=If[varDerResult=!=zero,expression,shiftedExpression];
(* Clear all variables not being passed back. *)
Clear[expression,zero,denominatorList,denominatorTermList,orderList,commonTerms,varDerResult,termsToShift,shiftedTerms];
Return[shiftedExpression]] 
(* end Module shiftDenominatorTerm *)

(*########## Function:homotopyFinalIntegration ##########*)
(************************************************************************)
(* homotopyFinalIntegration[expression,noOfDependentVar,indepVarList,   *)
(*     degCheck]                                                        *)
(* Purpose: To integrate an expression involving partial derivatives    *)
(*          in 1D and to invert the divergence operator in 2D and 3D.   *)
(* Input:   An expression                                               *)
(*          The number of dependent variables                           *)
(*          A list of independent variables                             *)
(* Output:  The inverted expression.                                    *)
(* Created: 6 March, 2006 by DP                                         *)
(* Last Modified: 9 April, 2008, 12:55 by DP at CSM                     *)
(************************************************************************)
homotopyFinalIntegration[expression_,noOfDependentVar_,indepVarList_,degCheck_Integer]:=Module[{numIndVars=Length[indepVarList],zero,lambdaPAR,expression2,indVar,uReplaceuLambdaRule,integrandList,sumOfDepComponents,lambdaSum,lambdaResult,homotopyAtOne,homotopyAtZero,interimResult,constantsk,divInterimResult,termsWithConstants,termsWOConstants,coefSystem,coefSystemSoln,homotopyResult,checkResult,dd,m1,m2,m3,m4,m5},
(* If the expression is zero,return 0 vector. *)
zero=Table[0,{i1,1,numIndVars}];
If[expression===0,Return[zero]];
(* Check to see if lambda is a parameter in the expression. If so,lambda to another parameter, lambdaPAR. *)expression2=If[MemberQ[expression,lambda,{0,Infinity}],expression/.lambda->lambdaPAR,expression,expression];
indVar[1]=Part[indepVarList,1];
If[numIndVars>1,indVar[2]=Part[indepVarList,2]];
If[numIndVars>2,indVar[3]=Part[indepVarList,3]];
(* The uReplaceuLambda rule takes u from an expression and replaces it with lambda*u. *)uReplaceuLambdaRule={u[m1_][m2__]->lambda*u[m1][m2],Derivative[m3__][u[m4_]][m5__]->lambda*Derivative[m3][u[m4]][m5]};
If[degCheck=!=0,
(* Use the integrand functions to determine the integrands for the homotopy operator. *)integrandList=Switch[numIndVars,1,Table[{integrandLeftDerList1D[expression2,j1,indepVarList]},{j1,1,noOfDependentVar}],2,Table[integrandLeftDerList2D[expression2,j1,indepVarList],{j1,1,noOfDependentVar}],3,Table[integrandLeftDerList3D[expression2,j1,indepVarList],{j1,1,noOfDependentVar}]];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>5.0*numIndVars,printOUT["Beginning to sum the expressions determined for each"<>" dependent variable."]];(* end If globalHomotopyVerbose *)
(* Add the terms produced by the integrand operator produced by each dependent variable. *)
sumOfDepComponents=Apply[Plus,integrandList];
If[globalHomotopyVerbose,printOUT["Proceding to the integration step."];
If[TimeUsed[]-timeKeeper>6.0*numIndVars,printOUT["All jet space coordinates are being replaced by "<>ToString[λ]<>"*coordinate."]] (*end If TimeUsed*)];(* end If globalHomotopyVerbose *)
(* Apply the lambda replace rule. *)
lambdaSum=sumOfDepComponents/.uReplaceuLambdaRule;
(* Divide all terms by lambda and integrate with respect to lambda. *)
lambdaSum=Expand[lambdaSum/lambda];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>6.0*numIndVars,printOUT["Beginning integration over "<>ToString[λ]<>"."]];(* end If globalHomotopyVerbose *)
lambdaResult=Map[Integrate[#,lambda]&,lambdaSum,{1}],
(* else *)
If[numIndVars===1,lambdaResult={Integrate[expression,indVar[1]]};
printOUT["Mathematica's Integrate function has been called because a zero degree term existed in the expression to be integrated."];
Return[lambdaResult],
(* else *)
printOUT["The homotopy operator is unable to invert the divergence on the given expression because at least one term in the"<>" expression has an overall degree of zero."];
Return[{}]];(* end If numIndVars *)];(* end If degCheck *)If[MemberQ[lambdaResult,_Integrate,{0,Infinity}],printOUT["The homotopy operator was unable to perform the auxillary integration required in the final step because the function is"<>" to complicated to integrate.  It is unable to provide a result."];
lambdaResult={}];(* end If MemberQ *)If[globalHomotopyVerbose,If[TimeUsed[]-timeKeeper>6.0*numIndVars,printOUT["Evaluating the result from "<>ToString[λ]<>" = 0"<>" to "<>ToString[λ]<>" = 1."];];(* end If TimeUsed *)
printOUT["------------------------------------------"]];(* end If globalHomotopyVerbose *)homotopyAtOne=If[MemberQ[lambdaResult,Sin[_]|Cos[_]|Tan[_]|Cot[_]|Sec[_]|Csc[_]|Sinh|Cosh[_]|Tanh[_]|Coth[_]|Sech[_]|Csc[_],{0,Infinity}]==={},Expand[TrigExpand[lambdaResult/.{lambda->1}]],If[MemberQ[lambdaResult,Log[_],{0,Infinity}],Expand[PowerExpand[lambdaResult/.{lambda->1}]],Expand[lambdaResult/.{lambda->1}]]];
(* Check the result after lambda->1 by comparing its divergence with the given expression. If they match,return and finish.*)checkResult=divergenceCheck[expression2,homotopyAtOne,indepVarList];
If[checkResult===0,
(* Clear all variables not being passed back. *)
Clear[numIndVars,zero,expression2,indVar,uReplaceuLambdaRule,lambdaSum,integrandList,sumOfDepComponents,lambdaResult,checkResult];
Return[homotopyAtOne/.lambdaPAR->lambda]];(* end If checkResult] *)
Off[Power::infy];
Off[∞::indet];
homotopyAtZero=If[MemberQ[lambdaResult,lambda,{0,Infinity}],Expand[lambdaResult/.lambda->0],zero];
(* Remove numerical constants. *)
homotopyAtZero=If[MemberQ[homotopyAtZero,u[_][__]|Derivative[__][u[_]][__],{0,Infinity}],homotopyAtZero,zero];
(* Mathematica may produce "constants" in terms of the dependent variables when it integrates the integrand expression over lambda. These extra terms must be removed before reverting from lambda space back to the independent variable space. First adjust the expression for any lambdas in the denominator. *)
homotopyAtZero=homotopyAtZero/.ComplexInfinity->0;
homotopyAtZero=DeleteCases[homotopyAtZero,Times[n1_,___]/;n1===Infinity||n1===-Infinity,{0,Infinity}];
homotopyAtZero=If[homotopyAtZero==={},zero,homotopyAtZero];
(* Second,check the expression with lambda -> 1 minus the expression with lambda -> 0. *)
interimResult=Expand[PowerExpand[Factor[homotopyAtOne-homotopyAtZero]]];
checkResult=divergenceCheck[expression2,interimResult,indepVarList];
If[checkResult===0,On[Power::infy];
On[∞::indet];
(* Clear all variables not being passed back. *)
Clear[numIndVars,zero,expression2,indVar,uReplaceuLambdaRule,lambdaSum,integrandList,sumOfDepComponents,lambdaResult,homotopyAtOne,homotopyAtZero,checkResult];
Return[interimResult/.lambdaPAR->lambda]];(* end If checkResult *)
(* Third,it may be necessary to evaluate terms produced when lambda -> 0. Not all terms in this case can removed, so comparisons must be made with the given expression. *)
If[homotopyAtZero=!=zero,homotopyAtZero=applyConstantsToVector[homotopyAtZero,numIndVars];
constantsk=Part[homotopyAtZero,4];
homotopyAtZero=If[numIndVars===1,Part[homotopyAtZero,1],Table[Part[homotopyAtZero,j],{j,1,numIndVars}]];
interimResult=Expand[homotopyAtOne-homotopyAtZero];
divInterimResult=If[numIndVars===1,D[Part[interimResult,1],indVar[1]],If[numIndVars===2,divergence[Flatten[{interimResult,0}],indVar[1],indVar[2],dd],divergence[interimResult,indVar[1],indVar[2],indVar[3]]]];
(* Compare the derivative/divergence of the result with the given expression. *)
divInterimResult=Factor[divInterimResult-expression];
divInterimResult=Expand[Numerator[divInterimResult]];
termsWithConstants=Apply[Plus,Cases[divInterimResult,Times[k[_],___],{0,Infinity}]];
termsWOConstants=-Complement[divInterimResult,termsWithConstants];
coefSystem=makeCoefficientEquations[termsWithConstants,constantsk,termsWOConstants,indepVarList,"caseTWO"];
coefSystemSoln=Flatten[Solve[coefSystem,constantsk]];
homotopyResult=interimResult/.coefSystemSoln/.k[_]->0;
homotopyResult=Factor[homotopyResult]];(* end If homotopyAtZero *)
On[Power::infy];
On[∞::indet];
(* Clear all variables not being passed back. *)
Clear[numIndVars,zero,expression2,indVar,uReplaceuLambdaRule,lambdaSum,integrandList,sumOfDepComponents,lambdaResult,homotopyAtOne,homotopyAtZero,interimResult,constantsk,divInterimResult,termsWithConstants,termsWOConstants,coefSystem,coefSystemSoln];
Return[homotopyResult/.lambdaPAR->lambda]];
(* end Module homotopyFinalIntegration *)

(*########## Function:integrandLeftDerList1D ##########*)
(************************************************************************)
(* integrandLeftDerList1D[expression,dependentVar,indepVarList]         *)
(* Purpose: To find the integrand for the 1-D homotopy operator         *)
(*          using the form where the derivatives of the dependent       *)
(*          varialble have been propagated outside the inner most sum   *)
(*          (after Kruskal et.al.) In this version,operators are        *)
(*          determined by manipulating lists                            *)
(* Input:   An expression                                               *)
(*          A number for the dependent variable                         *)
(*          A list of independent variables                             *)
(* Output:  The result of applying the higher Euler operators to the    *)
(*              expression.                                             *)
(* Created: 22 February, 2006 by DP                                     *)
(* Last Modified: 9 April, 2008, 12:40 by DP at CSM                     *)
(************************************************************************)integrandLeftDerList1D[expression_,dependentVar_,indepVarList_]:=Module[{derivativeList,orderList,i,maxi,partialWRTpartialList,innerSumList,dummyderiv,innerSum,outerSumList,integrandResult,k,indVar1=Part[indepVarList,1]},timeKeeper=TimeUsed[];
If[globalHomotopyVerbose,printOUT["Beginning the computation of the homotopy integrand for"<>" dependent variable "<>ToString[dependentVar]<>".  The algorithm"<>" first finds partial derivatives with respect to the jet space"<>" coordinates found in the expression, then performs a series of"<>" total derivatives on the partial derivatives.  The results are"<>" then combined into a sum."]];(* end If globalHomotopyVerbose *)(* Create a list of patial derivatives existing in the expression. *)derivativeList=Union[Cases[expression,Derivative[___][u[dependentVar]][indVar1],{0,Infinity}]];
(* Create a list of the orders of all derivatives that exist in the expression.*)
orderList=Flatten[Map[Head,Map[Head,derivativeList]]/.Derivative->List];
maxi=Max[orderList,0];
(* Create a list of the partial derivatives of the expression with respect to the partial derivative list created above.*)partialWRTpartialList=Map[D[expression,#]&,derivativeList];
(*Create a list of terms for the inner sum of the formula using the partial derivative list created in the last step.*)innerSumList=Table[Map[(-1)^(#-i-1)&,Select[orderList,#>i&]]*Thread[dummyderiv[Extract[partialWRTpartialList,Position[orderList,_?(#>i&)]],Map[{indVar1,#-i-1}&,Select[orderList,#>i&]]]]/.dummyderiv->D,{i,0,maxi-1}];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>2.0,printOUT["Beginning to sum all total derivatives."]];(* end If globalHomotopyVerbose *)(* Add the terms in the sum list. *)
innerSum=Expand[Apply[Plus,innerSumList,1]];
(* Multiply the sum by the derivatives of the dependent varialble with respect to x according to the outer sum. The result is a list of terms. *)outerSumList=Expand[Table[Part[innerSum,i+1]*Derivative[i][u[dependentVar]][indVar1],{i,0,maxi-1}]];
(* Add the terms from the list in the last step to get the result of applying the higher Euler operator. *)
integrandResult=Expand[Apply[Plus,outerSumList]];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>3.0,printOUT["All sums for dependent variable "<>ToString[dependentVar]<>" have been compiled."]];(* end If globalHomotopyVerbose *)
(* Clear all variables not being passed back. *)
Clear[derivativeList,orderList,i,k,maxi,partialWRTpartialList,innerSumList,dummyderiv,innerSum,outerSumList,indVar1];
Return[integrandResult]] 
(* end Module integrandLeftDerList1D *)

(*########## Function:integrandLeftDerList2D ##########*)
(************************************************************************)
(* integrandLeftDerList2D[expression,dependentVar,indepVarList]         *)
(* Purpose: To find the integrand for the 2-D homotopy operator         *)
(*          using the form where the derivatives of the dependent       *)
(*          varialble have been propagated outside the inner most sum   *)
(*          (after Kruskal et.al.) In this version,operators are        *)
(*          determined by manipulating lists                            *)
(* Input:   An expression                                               *)
(*          A number for the dependent variable                         *)
(*          A list of independent variables                             *)
(* Output:  The result of applying the higher Euler operators to the    *)
(*              expression.                                             *)
(* Created: 27 February, 2006 by DP                                     *)
(* Last Modified: 9 April, 2008, 12:42 by DP at CSM                     *)
(************************************************************************)
integrandLeftDerList2D[expression_,dependentVar_,indepVarList_]:=Module[{derivativeList,ix,iy,kx,ky,dummyderiv,orderList,maxix,maxiy,partialWRTpartialList,innerSumListx,innerSumListy,innerSumx,innerSumy,positionNonZeroTermsx,positionNonZeroTermsy,listWOzerosx,listWOzerosy,indVar1=Part[indepVarList,1],indVar2=Part[indepVarList,2],integrandResultx,integrandResulty,outerSumListx,outerSumListy},
timeKeeper=TimeUsed[];
If[globalHomotopyVerbose,printOUT["Beginning the computation of the homotopy integrand for"<>" dependent variable "<>ToString[dependentVar]<>".  The algorithm"<>" first finds partial derivatives with respect to the jet space"<>" coordinates found in the expression, then performs a series of"<>" total derivatives on the partial derivatives.  The results are"<>" then combined into a sum."]];(* end If globalHomotopyVerbose *)(* Create a list of patial derivatives existing in the expression. *)derivativeList=Union[Cases[expression,Derivative[___][u[dependentVar]][indVar1,indVar2],{0,Infinity}]];
(* Create lists of the orders of all derivatives that exist in the expression. *)
orderList=Map[Head,Map[Head,derivativeList]]/.Derivative->List;
(* Determine the highest order derivative for each variable. *)
maxix=Max[Map[Part[#,1]&,orderList],0];
maxiy=Max[Map[Part[#,2]&,orderList],0];
(* Create lists for each variable of the partial derivatives of the expression with respect to the partial derivative list created above. *)partialWRTpartialList=Map[D[expression,#]&,derivativeList];
(* Create lists of terms for the inner sum of the new formula using the partial derivative list created in the last step. *)innerSumListx=Table[Map[(-1)^(Part[#,1]+Part[#,2]-ix-iy-1)*Binomial[Part[#,1]+Part[#,2]-ix-iy-1,Part[#1,1]-ix-1]/Binomial[Part[#,1]+Part[#,2],Part[#,1]]&,Select[Select[orderList,Part[#,1]>ix&],Part[#,2]≥iy&]]*Thread[dummyderiv[Extract[partialWRTpartialList,Position[orderList,{_,_}?(Part[#,1]>ix&&Part[#,2]≥iy&)]],Map[{indVar1,Part[#,1]-ix-1}&,Select[Select[orderList,Part[#,1]>ix&],Part[#,2]≥iy&]],Map[{indVar2,Part[#,2]-iy}&,Select[Select[orderList,Part[#,1]>ix&],Part[#,2]≥iy&]]]]/.dummyderiv->D,{iy,0,maxiy},{ix,0,maxix-1}];
innerSumListy=Table[Map[(-1)^(Part[#,1]+Part[#,2]-ix-iy-1)*Binomial[Part[#,1]+Part[#,2]-ix-iy-1,Part[#,2]-iy-1]/Binomial[Part[#,1]+Part[#,2],Part[#,2]]&,Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]>iy&]]*Thread[dummyderiv[Extract[partialWRTpartialList,Position[orderList,{_,_}?(Part[#,1]≥ix&&Part[#,2]>iy&)]],Map[{indVar1,Part[#,1]-ix}&,Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]>iy&]],Map[{indVar2,Part[#,2]-iy-1}&,Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]>iy&]]]]/.dummyderiv->D,{iy,0,maxiy-1},{ix,0,maxix}];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>2.0,printOUT["Beginning to sum all total derivatives."]];(* end If globalHomotopyVerbose *)(* Add the terms in the sum lists to produce lists of inner sums corresponding to each outer index for ix and iy. *)innerSumx=Expand[Apply[Plus,innerSumListx,{2,2}]];
innerSumy=Expand[Apply[Plus,innerSumListy,{2,2}]];
(* Remove all zeros from the lists of sums above and keep track of the positions of the remaining terms. The positions give the outer index,ix and iy,for each term (in reverse order). *)
positionNonZeroTermsx=Position[innerSumx,_?(#=!=0&&#=!=List&),{2,2}];
positionNonZeroTermsy=Position[innerSumy,_?(#=!=0&&#=!=List&),{2,2}];
listWOzerosx=Flatten[DeleteCases[innerSumx,0,2]];
listWOzerosy=Flatten[DeleteCases[innerSumy,0,2]];
(* Multiply the sum by the partial derivative of the dependent variable with respect to the outer sum indices (ix and iy) and* constant determined by the outer sum indicies. *)outerSumListx=Thread[Times[listWOzerosx,Map[Binomial[Part[#,2]+Part[#,1]-2,Part[#,2]-1]*Derivative[Part[#,2]-1,Part[#,1]-1][u[dependentVar]][indVar1,indVar2]&,positionNonZeroTermsx]]];
outerSumListy=Thread[Times[listWOzerosy,Map[Binomial[Part[#,2]+Part[#,1]-2,Part[#,1]-1]*Derivative[Part[#,2]-1,Part[#,1]-1][u[dependentVar]][indVar1,indVar2]&,positionNonZeroTermsy]]];
(* Add the terms from the lists in the last step to get the result of applying the higher Euler operator for each variable. *)integrandResultx=Expand[Apply[Plus,outerSumListx]];
integrandResulty=Expand[Apply[Plus,outerSumListy]];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>3.0,printOUT["All sums for dependent variable "<>ToString[dependentVar]<>" have been compiled."]];(* end If globalHomotopyVerbose *)
(* Clear all variables not being passed back. *)
Clear[derivativeList,ix,iy,kx,ky,dummyderiv,orderList,maxix,maxiy,partialWRTpartialList,innerSumListx,innerSumListy,innerSumx,innerSumy,positionNonZeroTermsx,positionNonZeroTermsy,listWOzerosx,listWOzerosy,indVar1,indVar2,outerSumListx,outerSumListy];
Return[{integrandResultx,integrandResulty}]] 
(* end Module integrandLeftDerList2D *)

(*########## Function:integrandLeftDerList3D ##########*)
(************************************************************************)
(* integrandLeftDerList3D[expression,dependentVar,indepVarList]         *)
(* Purpose: To find the integrand for the 3-D homotopy operator         *)
(*          using the form where the derivatives of the dependent       *)
(*          varialble have been propagated outside the inner most sum   *)
(*          (after Kruskal et.al.) In this version,operators are        *)                   
(*          determined by manipulating lists                            *)
(* Input:   An expression                                               *)
(*          A number for the dependent variable                         *)
(*          A list of independent variables                             *)
(* Output:  The result of applying the higher Euler operators to the    *)
(*              expression.                                             *)
(* Created: 2 March, 2006 by DP                                         *)
(* Last Modified: 9 April, 2008, 12:43 by DP at CSM                     *)
(************************************************************************)

integrandLeftDerList3D[expression_,dependentVar_,indepVarList_]:=Module[{derivativeList,orderList,maxix,maxiy,maxiz,partialWRTpartialList,innerSumListx,innerSumListy,innerSumListz,ix,iy,iz,kx,ky,kz,innerSumx,innerSumy,innerSumz,positionNonZeroTermsx,positionNonZeroTermsy,positionNonZeroTermsz,listWOzerosx,listWOzerosy,listWOzerosz,dummyderiv,indVar1=Part[indepVarList,1],indVar2=Part[indepVarList,2],indVar3=Part[indepVarList,3],outerSumListx,outerSumListy,outerSumListz,integrandResultx,integrandResulty,integrandResultz},
timeKeeper=TimeUsed[];
If[globalHomotopyVerbose,printOUT["Beginning the computation of the homotopy integrand for"<>" dependent variable "<>ToString[dependentVar]<>".  The algorithm"<>" first finds partial derivatives with respect to the jet space"<>" coordinates found in the expression, then performs a series of"<>" total derivatives on the partial derivatives.  The results are"<>" then combined into a sum."]];(* end If globalHomotopyVerbos e*)(* Create a list of patial derivatives existing in the expression. *)derivativeList=Union[Cases[expression,Derivative[___][u[dependentVar]][indVar1,indVar2,indVar3],{0,Infinity}]];
(* Create lists of the orders of all derivatives that exist in the expression. *)
orderList=Map[Head,Map[Head,derivativeList]]/.Derivative->List;
(* Determine the highest order derivative for each variable. *)
maxix=Max[Map[Part[#,1]&,orderList],0];
maxiy=Max[Map[Part[#,2]&,orderList],0];
maxiz=Max[Map[Part[#,3]&,orderList],0];
(* Create lists for each variable of the partial derivatives of the expression with respect to the partial derivative list created above.*)partialWRTpartialList=Map[D[expression,#]&,derivativeList];
(* Create lists of terms for the inner sum of the higher Euler formula using the partial derivative list created in the last step.*)innerSumListx=Table[Map[(-1)^(Part[#,1]+Part[#,2]+Part[#,3]-ix-iy-iz-1)*Binomial[Part[#,1]+Part[#,2]+Part[#,3]-ix-iy-iz-1,Part[#,1]-ix-1]*Binomial[Part[#,2]+Part[#,3]-iy-iz,Part[#,2]-iy]/(Binomial[Part[#,1]+Part[#,2]+Part[#,3],Part[#,1]]*Binomial[Part[#,2]+Part[#,3],Part[#,3]])&,Select[Select[Select[orderList,Part[#,1]>ix&],Part[#,2]≥iy&],Part[#,3]≥iz&]]*Thread[dummyderiv[Extract[partialWRTpartialList,Position[orderList,{_,_,_}?(Part[#1,1]>ix&&Part[#1,2]≥iy&&Part[#,3]≥iz&)]],Map[{indVar1,Part[#,1]-ix-1}&,Select[Select[Select[orderList,Part[#,1]>ix&],Part[#,2]≥iy&],Part[#,3]≥iz&]],Map[{indVar2,Part[#,2]-iy}&,Select[Select[Select[orderList,Part[#,1]>ix&],Part[#,2]≥iy&],Part[#,3]≥iz&]],Map[{indVar3,Part[#,3]-iz}&,Select[Select[Select[orderList,Part[#,1]>ix&],Part[#,2]≥iy&],Part[#,3]≥iz&]]]]/.dummyderiv->D,{iz,0,maxiz},{iy,0,maxiy},{ix,0,maxix-1}];
innerSumListy=Table[Map[(-1)^(Part[#,1]+Part[#,2]+Part[#,3]-ix-iy-iz-1)*Binomial[Part[#,1]+Part[#,2]+Part[#,3]-ix-iy-iz-1,Part[#1,2]-iy-1]*Binomial[Part[#,1]+Part[#,3]-ix-iz,Part[#,3]-iz]/(Binomial[Part[#,1]+Part[#,2]+Part[#,3],Part[#,2]]*Binomial[Part[#,1]+Part[#,3],Part[#,3]])&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]>iy&],Part[#,3]≥iz&]]*Thread[dummyderiv[Extract[partialWRTpartialList,Position[orderList,{_,_,_}?(Part[#,1]≥ix&&Part[#,2]>iy&&Part[#,3]≥iz&)]],Map[{indVar1,Part[#,1]-ix}&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]>iy&],Part[#,3]≥iz&]],Map[{indVar2,Part[#,2]-iy-1}&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]>iy&],Part[#,3]≥iz&]],Map[{indVar3,Part[#,3]-iz}&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]>iy&],Part[#,3]≥iz&]]]]/.dummyderiv->D,{iz,0,maxiz},{iy,0,maxiy-1},{ix,0,maxix}];
innerSumListz=Table[Map[(-1)^(Part[#,1]+Part[#,2]+Part[#,3]-ix-iy-iz-1)*Binomial[Part[#1,1]+Part[#,2]+Part[#,3]-ix-iy-iz-1,Part[#,3]-iz-1]*Binomial[Part[#,1]+Part[#,2]-ix-iy,Part[#,1]-ix]/(Binomial[Part[#,1]+Part[#,2]+Part[#,3],Part[#,3]]*Binomial[Part[#,1]+Part[#,2],Part[#,1]])&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]≥iy&],Part[#,3]>iz&]]*Thread[dummyderiv[Extract[partialWRTpartialList,Position[orderList,{_,_,_}?(Part[#,1]≥ix&&Part[#,2]≥iy&&Part[#,3]>iz&)]],Map[{indVar1,Part[#,1]-ix}&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]≥iy&],Part[#,3]>iz&]],Map[{indVar2,Part[#,2]-iy}&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]≥iy&],Part[#,3]>iz&]],Map[{indVar3,Part[#,3]-iz-1}&,Select[Select[Select[orderList,Part[#,1]≥ix&],Part[#,2]≥iy&],Part[#,3]>iz&]]]]/.dummyderiv->D,{iz,0,maxiz-1},{iy,0,maxiy},{ix,0,maxix}];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>2.0,printOUT["Beginning to sum all total derivatives."]];(* end If globalHomotopyVerbose *)(* Add the terms in the sum lists to produce lists of inner sums corresponding to each outer index for ix and iy. *)innerSumx=Expand[Apply[Plus,innerSumListx,{3,3}]];
innerSumy=Expand[Apply[Plus,innerSumListy,{3,3}]];
innerSumz=Expand[Apply[Plus,innerSumListz,{3,3}]];
(* Remove all zeros from the lists of sums above and keep track of the positions of the remaining terms. The positions give the outer index,ix,iy,and iz,for each term (The positions in the list are given in the order z,y,x).*)positionNonZeroTermsx=Map[Part[#]-{1,1,1}&,Position[innerSumx,_?(#=!=0&&#=!=List&),{3,3}]];
positionNonZeroTermsy=Map[Part[#]-{1,1,1}&,Position[innerSumy,_?(#=!=0&&#=!=List&),{3,3}]];
positionNonZeroTermsz=Map[Part[#]-{1,1,1}&,Position[innerSumz,_?(#=!=0&&#=!=List&),{3,3}]];
listWOzerosx=Flatten[DeleteCases[innerSumx,0,3]];
listWOzerosy=Flatten[DeleteCases[innerSumy,0,3]];
listWOzerosz=Flatten[DeleteCases[innerSumz,0,3]];
(* Multiply the sum by the dependent variable and conduct derivatives with respect to x and y according to the outer sum.The result is a list of terms.*)outerSumListx=Thread[Times[Map[Binomial[Part[#,3]+Part[#,2]+Part[#,1],Part[#,3]]*Binomial[Part[#,2]+Part[#,1],Part[#,2]]&,positionNonZeroTermsx]*Map[Derivative[Part[#,3],Part[#,2],Part[#,1]][u[dependentVar]][indVar1,indVar2,indVar3]&,positionNonZeroTermsx],listWOzerosx]];
outerSumListy=Thread[Times[Map[Binomial[Part[#,3]+Part[#,2]+Part[#,1],Part[#,2]]*Binomial[Part[#,3]+Part[#,1],Part[#,1]]&,positionNonZeroTermsy]*Map[Derivative[Part[#,3],Part[#,2],Part[#,1]][u[dependentVar]][indVar1,indVar2,indVar3]&,positionNonZeroTermsy],listWOzerosy]];
outerSumListz=Thread[Times[Map[Binomial[Part[#,3]+Part[#,2]+Part[#,1],Part[#,1]]*Binomial[Part[#,3]+Part[#,2],Part[#,3]]&,positionNonZeroTermsz]*Map[Derivative[Part[#,3],Part[#,2],Part[#,1]][u[dependentVar]][indVar1,indVar2,indVar3]&,positionNonZeroTermsz],listWOzerosz]];
(* Add the terms from the lists in the last step to get the result of applying the higher Euler operator for each variable. *)integrandResultx=Expand[Apply[Plus,outerSumListx]];
integrandResulty=Expand[Apply[Plus,outerSumListy]];
integrandResultz=Expand[Apply[Plus,outerSumListz]];
If[globalHomotopyVerbose&&TimeUsed[]-timeKeeper>3.0,printOUT["All sums for dependent variable "<>ToString[dependentVar]<>" have been compiled."]];(* end If globalHomotopyVerbose *)
(* Clear all variables not being passed back. *)
Clear[derivativeList,orderList,maxix,maxiy,maxiz,partialWRTpartialList,innerSumListx,innerSumListy,innerSumListz,ix,iy,iz,kx,ky,kz,innerSumx,innerSumy,innerSumz,positionNonZeroTermsx,indVar2,indVar3,positionNonZeroTermsy,positionNonZeroTermsz,listWOzerosx,listWOzerosy,listWOzerosz,dummyderiv,indVar1,outerSumListx,outerSumListy,outerSumListz];
Return[{integrandResultx,integrandResulty,integrandResultz}]] 
(* end Module integrandLeftDerList3D *)

(**********************************************************************)
(*                                                                    *)
(*    BEGINNING THE REMOVAL OF CURL TERMS FROM THE HOMOTOPY RESULT    *)
(*                                                                    *)
(**********************************************************************)

(*########## Function:removeCurlTerms ##########*)
(**********************************************************************)
(* removeCurlTerms[resultFromHomotopy, independVarlist]               *)
(* Purpose: To reduce the number of terms in the final homotopy       *)
(*          operator result to least amount possible.This is a        *)
(*          modified form of extractSystemForCoefficient              *)
(* Input:   The result from the homotopy operator                     *)
(*          A list of independent variables                           *)
(* Output:  A simplified result                                       *)
(* Created: 20 September, 2006 by DP at home                          *)
(* Last Modified: 4 April, 2008, 15:44 by DP at CSM                   *)
(**********************************************************************)
removeCurlTerms[resultFromHomotopy_,allIndependVarList_]:=Module[{lengthIndVar,lengthComp,unknownConstants,independVarList,homotopyExpnWithConstants,divergenceHomotopyWithK,divCheckDiv,divergenceHomotopyNoK,system,zeroEquations,zeroConstants,j1=1,j2,nonzeroConstants,solveConstantList,systemSoln,finalHomotopy,keepHomotopy,divCheck,ivar,i,i1=1,homotopyResultPart,commonDenominator},
(* If no result is given,return immediately. *)
If[resultFromHomotopy==={},Return[resultFromHomotopy]];
Clear[k];
printOUT["The homotopy operator has finished inverting the divergence.  The program is beginning the removal of curl terms."];
(* Determine the number of independent variables. *)
independVarList=Part[allIndependVarList,1];
lengthIndVar=Length[resultFromHomotopy];
(* Assign independent variables to consistent terms. *)
ivar[1]=Part[independVarList,1];
If[lengthIndVar>1,ivar[2]=Part[independVarList,2]];
If[lengthIndVar>2,ivar[3]=Part[independVarList,3]];
(* Find the number of terms in each component of the homotopy result. *)lengthComp[1]=If[Head[Part[resultFromHomotopy,1]]===Plus,Length[Part[resultFromHomotopy,1]],1];
lengthComp[2]=If[Head[Part[resultFromHomotopy,2]]===Plus,Length[Part[resultFromHomotopy,2]],1];
lengthComp[3]=If[lengthIndVar===3,If[Head[Part[resultFromHomotopy,3]]===Plus,Length[Part[resultFromHomotopy,3]],1],0];
(* If only one term exists in each component,there is no need to go through the reduction process. *)If[lengthComp[1]≤1&&lengthComp[2]≤1&&lengthComp[3]≤1,Return[resultFromHomotopy]];(* end If lengthPart *)If[globalHomotopyVerbose,printOUT["The inverted divergence vector given by the homotopy operator"<>" has ",lengthIndVar," components.  The first component has"<>" ",lengthComp[1],If[lengthComp[1]===1," term"," terms"],If[lengthIndVar===3,", "," and "],"the second component has ",lengthComp[2],If[lengthComp[2]===1," term"," terms"],If[lengthIndVar===3,", and the third component has "<>ToString[lengthComp[3]]<>If[lengthComp[3]===1," term."," terms."],"."]]];(* end If globalHomotopyVerbose *)printOUT["If the number of terms in the components is large, say several hundred, the removal of curl terms may take a long time."];
printOUT["------------------------------------------"];
(* Check for rational expressions. *)denom=MemberQ[resultFromHomotopy,Power[n2__,n1_]/;n1<0/;MemberQ[n2,u[_][__],{0,Infinity}]||MemberQ[n2,Derivative[__][u[_]][__],{0,Infinity}],{0,Infinity}];
(* Distribute constants through the terms of the homotopy result. *){homotopyExpnWithConstants[1],homotopyExpnWithConstants[2],homotopyExpnWithConstants[3],unknownConstants}=applyConstantsToVector[resultFromHomotopy,lengthIndVar];
(* Take the divergence on the original homotopy result and the homotopy result with constants. *)If[lengthIndVar===2,divergenceHomotopyWithK=divergence[{homotopyExpnWithConstants[1],homotopyExpnWithConstants[2],0},ivar[1],ivar[2],ivar[3]];
divergenceHomotopyNoK=divergence[Append[resultFromHomotopy,0],ivar[1],ivar[2],ivar[3]]];(* end If lengthIndVar===2 *)If[lengthIndVar===3,divergenceHomotopyWithK=divergence[{homotopyExpnWithConstants[1],homotopyExpnWithConstants[2],homotopyExpnWithConstants[3]},ivar[1],ivar[2],ivar[3]];
divergenceHomotopyNoK=divergence[resultFromHomotopy,ivar[1],ivar[2],ivar[3]]];(* end If lengthIndVar===3 *)divergenceHomotopyNoK=Expand[Factor[divergenceHomotopyNoK]];
divCheckDiv=divergenceHomotopyNoK;
If[denom,commonDenominator=takeOnlyNumerator[MapAll[Apart,divergenceHomotopyWithK]];
divergenceHomotopyWithK=Map[Factor,Apply[List,Map[Times[#*commonDenominator]&,divergenceHomotopyWithK,{1}]]];
divergenceHomotopyWithK=Apply[Plus,divergenceHomotopyWithK];
divergenceHomotopyNoK=Map[Factor,Apply[List,Map[Times[#*commonDenominator]&,divergenceHomotopyNoK,{1}]]];
divergenceHomotopyNoK=Apply[Plus,divergenceHomotopyNoK]];(* end If denom *)
(* Make a system of equations on the unknown constants. The first case will try to construct a small system,assuming that most of the k[i] are zero.This may lead to an inconsistent system in some situations.If so,then case two will be applied. *)system=makeCoefficientEquations[divergenceHomotopyWithK,unknownConstants,divergenceHomotopyNoK,Flatten[allIndependVarList],"caseONE"];
If[Not[MemberQ[system,False,{0,Infinity}]],solveConstantList=Cases[system,k[_],{0,Infinity}];
zeroConstants=Complement[unknownConstants,solveConstantList];
zeroConstants=Map[Rule[#,0]&,zeroConstants];
(* Find two sets of solutions by reversing the order of the solve list. These will be compared for balance. *)systemSoln[1]=solveCoefficientSystem[system,solveConstantList,zeroConstants];
systemSoln[2]=solveCoefficientSystem[system,Reverse[solveConstantList],zeroConstants];
If[systemSoln[1]===systemSoln[2],j2=1,j2=2];
While[j1≤j2,
(* Replace constants in the homotopy operator with their solutions. *)finalHomotopy[j1]=Flatten[{homotopyExpnWithConstants[1],homotopyExpnWithConstants[2],homotopyExpnWithConstants[3]}]/.systemSoln[j1];
j1++];(* end While *)
If[j2===2,keepHomotopy=compareCurlFreeSolutions[finalHomotopy[1],finalHomotopy[2]],(*else*)keepHomotopy=finalHomotopy[1]],(* end If j2 *)
(* else *)
system=makeCoefficientEquations[divergenceHomotopyWithK,unknownConstants,divergenceHomotopyNoK,Flatten[allIndependVarList],"caseTWO"];
(* Identify all equations in the system where the right hand side is equal to zero.*)
zeroEquations=Cases[system,Equal[_,0],{0,2}];
(* Make a list of all constants,k[i],in zeroEquations. *)
zeroConstants=Union[Cases[zeroEquations,k[_],{0,Infinity}]];
(* Make a list of all constants that don't show up in zeroEquations *)
nonzeroConstants=Complement[unknownConstants,zeroConstants];
(* Solve the system,solving for the constants that don't appear in the zeroEquations first,then solving for all other constants.  Two different solve orders will be applied to give choices for the solution.*)solveConstantList=Flatten[Append[nonzeroConstants,zeroConstants]];
systemSoln[1]=solveCoefficientSystem[system,solveConstantList,{}];
(* Reverse the order of the constants from the nonzero equations and solve the system again. This gives the alternative. *)solveConstantList=Flatten[Append[Reverse[nonzeroConstants],zeroConstants]];
systemSoln[2]=solveCoefficientSystem[system,solveConstantList,{}];
If[systemSoln[1]===systemSoln[2],j2=1,j2=2];
While[j1≤j2,
(*Replace constants in the homotopy operator with their solutions.*)finalHomotopy[j1]=Flatten[{homotopyExpnWithConstants[1],homotopyExpnWithConstants[2],homotopyExpnWithConstants[3]}]/.systemSoln[j1];
j1++];(* end While *)
(* Compare simplified results. Report the result with the most alanced set of terms across the components. *)If[j2===2,keepHomotopy=compareCurlFreeSolutions[finalHomotopy[1],finalHomotopy[2]],
(* else *)
keepHomotopy=finalHomotopy[1]] (* end If j2 *)];(* end If MemberQ *)
(* Do a divergence check on the result after curl terms have been removed. If the test fails, return the result given by the homotopy operator. *)divCheck=divergenceCheck[divCheckDiv,keepHomotopy,independVarList];
If[divCheck=!=0||keepHomotopy==={},printOUT["The program was unable to remove curl terms from the vector provided by the homotopy operator.  The result with curl terms is begin returned."];
keepHomotopy=resultFromHomotopy];(* end divCheck *)
Clear[lengthIndVar,lengthComp,homotopyResultList,unknownConstants,termsWithConstants,homotopyExpnWithConstants,divergenceHomotopyWithK,divergenceHomotopyNoK,system,zeroEquations,zeroConstants,j1,j2,nonzeroConstants,solveConstantList,systemSoln,finalHomotopy,divCheck,commonDenominator,ivar,i,divCheckDiv];
Return[keepHomotopy];];
(* end Module removeCurlTerms *)

(*########## Function:makeCoefficientEquations ##########*)
(**********************************************************************)
(* makeCoefficientEquations[expressionWithC,unknownCoefList,          *)
(*     expressionNoC,indVarList,case]                                 *)
(* modified form of extractSystemForCoefficients                      *)
(* Purpose: To extract a linear system to solve for the coefficients  *)
(*          and match the coefficients to coefficients in an          *)
(*          equivalent system.                                        *)
(* Input:   The divergence of the homotopy expression after           *)
(*              coefficients,c[i],have been added                     *)
(*          The divergence of the homotopy expression with no         *)
(*              coefficients                                          *)
(*          List of unknown coefficients                              *)
(*          A list of independent variables                           *)
(*          A string telling whether the long or short method is      *)
(*              being used                                            *)
(* Output:  A system of coefficient equations                         *)
(* Created: 17 September, 2006 by DP at CSM                           *)
(* Last Modified: 10 April, 2008, 9:06 by DP at CSM                   *)
(**********************************************************************)makeCoefficientEquations[expressionWithK_,unknownCoefList_,expressionNoK_,indVarList_,case_String]:=Module[{expandedExpressionWithK,expandedExpressionNoK,termsToRemove,dco,commonTermsListWithK,commonTermsListNoK,coefficientEqns,zeroConstants,n1,n2,n3,indVarToZeroRules,excludeList},
expandedExpressionWithK=Expand[expressionWithK];
expandedExpressionNoK=Expand[expressionNoK];
(* Make a list of terms needed to identify common coefficients. The first case will form a system for the constants where the rhs is never zero. This may or may not form a consistent system. The second case sets up the whole coefficient system and will be used when the first case fails. *)Switch[case,"caseONE",termsToRemove=Reverse[Union[homotopyStripper[expandedExpressionNoK,indVarList]]],"caseTWO",termsToRemove=Reverse[Union[homotopyStripper[expandedExpressionWithK,indVarList]]]];(* end Switch *)excludeList=Flatten[{u[_][__],Derivative[__][u[_]][__],indVarList,Map[Power[#,_]&,indVarList]}];
(* Match all terms in the expression with k[i] according to the terms given in the previous list. *)rememberCommonTerms=commonTermsListWithK=Apply[Plus,Map[Cases[expandedExpressionWithK,Times[___?(Not[MemberQ[#,Apply[Alternatives,excludeList],{0,Infinity}]]&),#],{0,2}]&,termsToRemove],{1}];
(* Match all terms in the expression without c[i] according to the terms given in the previous list. *)expandedExpressionNoK=Map[Times[#,dco]&,expandedExpressionNoK];
commonTermsListNoK=Apply[Plus,Map[Cases[expandedExpressionNoK,Times[___?(Not[MemberQ[#,Apply[Alternatives,excludeList],{0,Infinity}]]&),#],{0,2}]&,termsToRemove],{1}];
commonTermsListWithK=Expand[Thread[Times[commonTermsListWithK,Map[Power[#,-1]&,termsToRemove]]]];
commonTermsListNoK=Expand[Thread[Times[commonTermsListNoK,Map[Power[#*dco,-1]&,termsToRemove]]]];
(* Set the corresponding expressions from each list equal to each other,forming coefficient equations. *)coefficientEqns=Thread[Equal[commonTermsListWithK,commonTermsListNoK]];
If[case==="caseONE",remainingTerms=Complement[Apply[List,expandedExpressionWithK],rememberCommonTerms/.Plus->Sequence];
zeroConstants=Union[Cases[remainingTerms,k[_],{0,Infinity}]];
zeroConstants=Map[Rule[#,0]&,zeroConstants];
coefficientEqns=coefficientEqns/.zeroConstants;
coefficientEqns=coefficientEqns/.Equal[0,_]->False];(* end If case *)Clear[expandedExpressionWithK,expandedExpressionNoK,termsToRemove,tim,commonTermsListWithK,commonTermsListNoK,zeroConstants,r,n1,n2,n3];
Return[coefficientEqns];];
(* end Module makeCoefficientEquations *)

(*########## Function:takeOnlyNumerator ##########*)
(************************************************************************)
(* takeOnlyNumerator[expression]                                        *)
(* Purpose: To rewrite a rational expression in terms of a common       *)
(*          denominator,then take the numerator.                        *)
(* Input:   An expression in rational form.                             *)
(* Output:  The numerator of the expression                             *)
(* Created: 31 August, 2007 by DP                                       *)
(* Last Modified: 10 March, 2008, 15:53 by DP at CSM                    *)
(************************************************************************)takeOnlyNumerator[expression_]:=Module[{workingExpression,coefs,lcm,denominatorTerms,coefsNonNum,firstTerm,firstTermList,highestDegreeTerm,commonDenom},workingExpression=ExpandAll[expression];
workingExpression=Apply[List,expression];
denominatorTerms=Map[Denominator,workingExpression];
denominatorTerms=Flatten[Map[FactorList,denominatorTerms],1];
denominatorTerms=Union[Apply[Power,denominatorTerms,{1}]];
(* Find the LCM of all numeric values. *)
coefs=Cases[denominatorTerms,n1_/;NumberQ[n1]];
coefs=If[coefs==={},{1},coefs];
lcm=Apply[LCM,coefs];
coefsNonNum=Complement[Cases[denominatorTerms,n1_/;NumericQ[n1]],coefs];
commonDenom=lcm*Apply[Times,coefsNonNum];
denominatorTerms=Complement[denominatorTerms,Union[coefs,coefsNonNum]];
(* Find the highest degree term for each type of differential terms found in the denominator. *)While[denominatorTerms=!={},firstTerm=Part[denominatorTerms,1]/.If[Head[Part[denominatorTerms,1]]===Power,Power[n1_,_]->n1,q->q];
firstTermList=Union[Cases[denominatorTerms,firstTerm],Cases[denominatorTerms,Power[firstTerm,___]]];
highestDegreeTerm=Part[Cases[firstTermList,Power[_,Max[Exponent[firstTermList,firstTerm]]]],1];
commonDenom=commonDenom*highestDegreeTerm;
denominatorTerms=Complement[denominatorTerms,firstTermList]];(* end While *)Clear[workingExpression,coefs,lcm,denominatorTerms,coefsNonNum,firstTerm,firstTermList,highestDegreeTerm];
Return[commonDenom]];
(* end Module takeOnlyNumerator *)

(*########## Function:solveCoefficientSystem ##########*)
(**********************************************************************)
(* solveCoefficientSystem[coefSystem,coefOrder,knownZeros]            *)
(* Purpose: To solve a linear system of constants. In this case, all  *)
(*          free constants in the solution will be set to zero        *)
(*          automatically. Zero solutions will be listed.             *)
(* Input:   The system to be solved                                   *)
(*          A list of constants to be solved for given in the order   *)
(*              that they should be solved for                        *)
(*          A list of solutions known to be zero                      *)
(* Output:  The solution to the system with every constant accounted  *)
(*              for.                                                  *)
(* Created: 21 May, 2007 by DP at CSM                                 *)
(* Last Modified: 10 March, 2008, 16:05 by DP at CSM                  *)
(**********************************************************************)solveCoefficientSystem[coefSystem_,coefOrder_List,knownZeros_List]:=Module[{systemSoln,undeterminedConstants,n1},
(* Solve the system given to the function. *)
systemSoln=Flatten[Solve[coefSystem,coefOrder]];
(* Collect all undetermined constants. *)
undeterminedConstants=Union[Cases[systemSoln/.Rule[_,n1_]->n1,k[_],{0,Infinity}]];
(* Formulate rules setting all undetermined constants to zero. *)
undeterminedConstants=Map[Rule[#,0]&,undeterminedConstants];
(* Apply the undetermined constants rules to the solution. *)
systemSoln=systemSoln/.undeterminedConstants;
(* Make a list showing the solutions to all k[i]. *)
systemSoln=Union[systemSoln,undeterminedConstants,knownZeros];
Clear[undeterminedConstants,n1];
Return[systemSoln]];
(* end Module solveCoefficientSystem *)

(*########## Function:compareCurlFreeSolutions ##########*)
(**********************************************************************)
(* compareCurlFreeSolutions[invertedDiv1,invertedDiv2]                *)
(* Purpose: To compare two possibilities for the final form for an    *)
(*          inverted divergence and recommend the one that has the    *)
(*          most balanced solution.                                   *)
(* Input:   A pair of inverted divergences                            *)
(* Output:  The inverted divergence with the most balanced form.      *)
(* Created: 21 May, 2007 by DP at CSM                                 *)
(* Last Modified: 10 March, 2008, 17:13 by DP at CSM                  *)
(**********************************************************************)compareCurlFreeSolutions[invertedDiv1_,invertedDiv2_]:=Module[{mylength,lengthCheck,findDiffInLength,lengthDiffs,i1,i2,keepHomotopy},
(* Compare simplified results. Report the result with the most balanced set of terms across the components. *)
mylength:=If[Head[#]===Plus,Length[#],1]&;
lengthCheck[1]=Map[mylength,invertedDiv1];
lengthCheck[2]=Map[mylength,invertedDiv2];
findDiffInLength:=Abs[#1-#2]&;
lengthDiffs[1]=Table[Part[Append[lengthCheck[1],Part[lengthCheck[1],1]],i1],{i2,1,Length[lengthCheck[1]]-1},{i1,i2,i2+1}];
lengthDiffs[2]=Table[Part[Append[lengthCheck[2],Part[lengthCheck[2],1]],i1],{i2,1,Length[lengthCheck[2]]-1},{i1,i2,i2+1}];
lengthDiffs[1]=Max[Apply[findDiffInLength,lengthDiffs[1],{1}]];
lengthDiffs[2]=Max[Apply[findDiffInLength,lengthDiffs[2],{1}]];
keepHomotopy=If[lengthDiffs[1]≤lengthDiffs[2],invertedDiv1,invertedDiv2];
Clear[mylength,lengthCheck,findDiffInLength,lengthDiffs,i1,i2];
Return[keepHomotopy]];
(* end Module compareCurlFreeSolutions *)

(*########## Function:divergence ##########*)divergence[expr_List,iv1_,iv2_:0,iv3_:0]:=D[Part[expr,1],iv1]+If[iv2===0,0,D[Part[expr,2],iv2]]+If[iv3===0,0,D[Part[expr,3],iv3]];

(*########## Function:divergenceCheck ##########*)
(************************************************************************)
(* divergenceCheck[expression,homotopyOpResult,indepVarList]            *)
(* Purpose: To take the derivative of divergence on the result          *)
(*          provided by the homotopy operator and check to see          *)
(*          if it is the same as the given expression.                  *)
(* Input:   The original given expression                               *)
(*          The calculated result from the homotopy operator            *)
(*          A list of independent variables                             *)
(* Output:  Yes or No.                                                  *)
(* Created: 30 August, 2007 by DP                                       *)
(* Last Modified: 20 March, 2008, 11:30 by DP at CSM                    *)
(************************************************************************)
divergenceCheck[expression_,homotopyOpResult_,indepVarList_]:=Module[{numIndVars=Length[indepVarList],newDiv,divCheck},newDiv=divergence[homotopyOpResult,Apply[Sequence,indepVarList]];
divCheck=Expand[expression-newDiv];
If[divCheck=!=0,divCheck=Factor[divCheck]];(* end If divCheck *)
If[divCheck=!=0,divCheck=Expand[TrigToExp[divCheck]]];(* end If divCheck *)
If[divCheck=!=0,divCheck=FullSimplify[divCheck,TimeConstraint->100]];(* end If divCheck *)
Clear[numIndVars,newDiv];
If[divCheck===0,Return[0],Return[1]]];
(* end Module divergenceCheck *)

(*########## Function:pdeform ##########*)
(**********************************************************************)
(* Purpose: Changes an expression: removes functional dependence, and *)
(*          prints derivatives as subscripts.                         *)
(* Authors: Tracy Otto, Tony Miller,and U.Goktas                      *)
(* Input:   An expression                                             *)
(* Output:  Same expression with all of the changes made              *)
(* Created: 1995 at CSM                                               *)
(*Last Modified: 19 December, 2006, 11:40 by DP at CSM                *)
(**********************************************************************)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 *)

(*########## Function:printSystemInPdeForm ##########*)
(**********************************************************************)
(* printSystemInPdeForm[]:                                            *)
(* Author:Mike Colagrosso                                             *)
(* Major Changes Made: 15 May, 2007, 12:03 by DP at CSM               *)
(**********************************************************************)
(* Major changes to function made since differential expressions and independent variables are passed by the function now. *)printSystemInPdeForm[diffFunctionList_,name_,note_,numDepVars_]:=Module[{inPDEform,eqnIter},Print[If[Length[diffFunctionList]>1,"System given for ",""],name,"."];
If[Head[note]===String,Print["NOTE: ",note]];
inPDEform=Map[pdeform,diffFunctionList];
For[eqnIter=1,eqnIter≤numDepVars,eqnIter++,Print["\nEquation ",eqnIter," of the system with ",numDepVars," equation",If[Length[diffFunctionList]>1,"s",""],":"];
Print["\n  ",Part[inPDEform,eqnIter]," = 0\n\n"];];(* end For *)]
 (* end Module *)

(*########## Function:homotopyDriver ##########*)
(************************************************************************)
(* homotopyDriver[expression, indepVariableList, (optional) parameters] *)
(* Purpose: To integrate a 1D expression or to invert the divergence    *)
(*          operator on a 2D or 3D expression by applying a homotopy    *)
(*          operator.                                                   *)
(* Input:   An expression                                               *)
(*          A list of independent variables                             *)
(*          A list of parameters (can be given as an option)            *)
(* Output:  The result of applying the higher Euler operators to the    *)
(*              expression.                                             *)
(* Created: 23 March, 2006 by DP                                        *)
(* Last Modified: 2 October,2008, 11:20 by DP at CSM                    *)
(************************************************************************)

homotopyDriver[expression_,indepVariableList_,parameters_:{}]:=Module[{expressionAltered,integrability,interimResult,indepVariableList2,integrableExpression,finalResult,lengthIndVarList,expressionDenom,finalCheck,ddd,j,newIndVarList},
paramrules={};
If[globalHomotopyVerbose,printOUT["The expression given to the program is"];
printOUT[expression];
printOUT["------------------------------------------"]];(* end globalHomotopyVerbose *)
(* Check to see that the independent variables are in a list. *)indepVariableList2=If[Head[indepVariableList]===Symbol,{indepVariableList},indepVariableList,indepVariableList];
lengthIndVarList=Length[indepVariableList2];
(* If the expression given to the homotopy operator is 0,Return 0. *)
If[expression===0,Return[Switch[lengthIndVarList,1,0,2,{0,0},3,{0,0,0}]]];
(* Reformat the given expression so that all dependent variables have the form u[_] and and any auxillary independent variables are renamed so that they appear to be removed from the expression. *)
expressionAltered=expressionVarReassign[expression,indepVariableList2];
(* Determine if the expression is integrable (1D) or if it is a divergence (2D and 3D). *)integrability=integrabilityTest[Part[expressionAltered,1],Part[expressionAltered,3],Part[expressionAltered,5],parameters];
(* If the expression is not integrable or not a divergence, stop the program. *)
If[Part[integrability,1]=!="Integrable",Abort[]];
(* If the expression is integrable or is a divergence,pass it to the homotopy operator to integrate it or invert the divergence. *)integrableExpression=If[Part[integrability,2]=!={},Part[expressionAltered,1]/.Part[integrability,2],Part[expressionAltered,1],Part[expressionAltered,1]];
interimResult=rationalExpressionCheck[integrableExpression,Part[expressionAltered,5],Part[expressionAltered,3]];
(* If homotopy operator did not work,return empty list. *)
If[interimResult==={},Return[{}]];(* end If interimResult *)
(* If the interim result is an inverted divergence,remove the curl terms. *)If[lengthIndVarList>1,newIndVarList={Part[expressionAltered,5],Complement[Part[expressionAltered,4],Part[expressionAltered,5]]};
finalResult=removeCurlTerms[ExpandAll[interimResult],newIndVarList],
(* else *)
finalResult=interimResult;];(* end If lengthIndVarList *)
(* Revert the expression back to the variables given in the problem. Part[expressionAltered, 2] is the inversion rule for the dependent variables. Part[expressionAltered, 6] returns terms with derivatives on auxillary variables. *)If[Part[expressionAltered,6]=!={},finalResult=returnAuxIndependentVars[finalResult,Part[expressionAltered,4],Part[expressionAltered,5],Part[expressionAltered,6]]];(* end If Part[expressionAltered, 6] *)
(* Return the final result in simplified factored form. *)
finalResult=MapAll[Apart,finalResult];
finalResult=finalResult/.Part[expressionAltered,2];
If[Part[integrability,2]=!={},Print["In the result that follows, parameters have been reassigned under the rules given by:"];
Print[Part[integrability,2]];
paramrules=Part[integrability,2]];(* end If Part[integrability, 2] *)
(* Check the final result by comparing its derivative/divergence to the expression given to the program.They should match. *)finalCheck=divergenceCheck[expression/.paramrules,finalResult,Part[expressionAltered,5]];
If[finalCheck===1,printOUT["The program has been unable to verify that the"<>If[lengthIndVarList===1," derivative"," divergence"]<>" of the result determined by the homotopy operator matches that of the given expression.  The result may not be accurate"<>" or the program may not be able to correctly simplify the expression calculated in the final check."],
(* else *)
printOUT["Check: Given Expression - "<>If[lengthIndVarList===1,"D[","Div["]<>"Homotopy Result] = ",finalCheck]];(* end If finalCheck *)finalResult=If[lengthIndVarList===1,Part[finalResult,1],finalResult,finalResult];
printio["Leaving homotopyDriver with final result:"];
printio[finalResult];
Clear[expressionAltered,integrability,interimResult,indepVariableList2,integrableExpression,finalCheck,lengthIndVarList,expressionDenom,ddd,j];
Return[finalResult]];
(* end Module homotopyDriver *)

(**********************************************************************)
(* The list of Options controls print and debug statements.           *)
(* The option is given as a rule in the opening function of the       *)
(* program.                                                           *)
(**********************************************************************)Options[HomotopyOperator]={Verbose->False,printHomotopyDriver->False,printStandAloneOutput->True};

(*########## Function:HomotopyOperator ##########*)
(**********************************************************************)
(* HomotopyOperator[pde,indepVariables,parameters,options]            *)
(*                                                                    *)
(* This function sets all options,then starts the program by calling  *)
(* the driver function.                                               *)
(**********************************************************************)
HomotopyOperator[pde_,indepVariables_,parameters___?(MemberQ[{List,Symbol},Head[#]]&),opts___?OptionQ]:=Module[{result},Off[Solve::svars];
(* Verbose -> True prints a running dialog showing the progression of a function through the homotopy operator. *)globalHomotopyVerbose=Verbose/.{opts}/.Options[HomotopyOperator];
(* printStandAloneOutput -> False suppresses all output. *)
printOUT=printStandAloneOutput/.{opts}/.Options[HomotopyOperator];
If[printOUT,printOUT=Print,Clear[printOUT],Clear[printOUT]];
result=homotopyDriver[pde,indepVariables,parameters];
If[printOUT===Print,Off[Solve::svars],On[Solve::svars]];
Return[result]] 
(* end Module homotopyOperator *)
Print["Code Homotopy.m successfully loaded."];
(* **************************end of all******************************* *)

Overview

Introduction to the homotopy operator

Introduction to conservation laws with focus on the calculation of the flux - A need for the homotopy operator arises

The homotopy operator and how it works

Examples using the homotopy operator

The Homotopy Operator

The homotopy operator is a tool from differential geometry.

The homotopy operator works on exact differential expressions involving unspecified functions of one, two, or three independent variables.

On expressions with one independent variable, the homotopy operator replaces integration by parts with a series of differentiations followed by a single integration over an auxillary variable.

On expressions with two or three independent variables, the homotopy operator inverts an expression known to be a divergence.

Conservation Laws

A conservation law for a partial differential equation (PDE) in (1 + 1) dimensions satisfies

(Subscript[D, t]) ρ + (Subscript[D, x]) J = 0    (on the PDE)

where ρ is the conserved density and J is the corresponding flux.

The Korteweg-de Vries (KdV) Equation

(Subscript[u, t]) + u (Subscript[u, x]) + (Subscript[u, xxx]) = 0

where u = u(x, t)

A well known, well studied equation, the KdV equation has an infinite number of conservation laws.  The first three are shown.

Conservation Laws for the KdV Equation

Conservation of momentum:

(Subscript[D, t]) ( u ) + (Subscript[D, x]) ( (1/2)(u^2) + (Subscript[u, xx]) ) = 0

Conservation of energy:

(Subscript[D, t]) ( (u^2) ) + (Subscript[D, x]) ( (2/3)(u^3) - (Subscript[u, x]^2) + 2u (Subscript[u, xx]) ) = 0

Boussinesq's moment of instability:

(Subscript[D, t]) ( (u^3) - 3(Subscript[u, x]^2) ) + (Subscript[D, x]) ( (3/4)(u^4) - 6u (Subscript[u, x]^2) + 3(u^2)(Subscript[u, xx]) +3(Subscript[u, xx]^2) - 6(Subscript[u, x])(Subscript[u, xxx])) = 0

An Algorithm for Conservation Laws

Our Focus: Find the Flux

The density ρ is calculated by building an expression based on symmetries found in the PDE. Our program, CondenMD finds the density symbollically

rho=u[x,t]^3 - 3 D[u[x,t],x]^2

Once the density is known, the flux, J, is calculated using the following algorithm based on the continuity equation.

1.  Find (Subscript[D, t]) ρ.

tDerivative=D[rho,t]

2.  Replace all t-derivatives in (Subscript[D, t]) ρ using the given PDE.
      For the KdV equation, (Subscript[u, t]) -> -(u (Subscript[u, x]) + (Subscript[u, xxx])).

modtDerivative=tDerivative/.Derivative[n_,1][u][x,t]->-(D[u[x,t]*D[u[x,t],x]+D[u[x,t],{x,3}],{x,n}])

3.  From the continuity equation, (Subscript[D, t]) ρ is the same as -(Subscript[D, x])J, so the flux can be obtained by integrating -(Subscript[D, t]) ρ where terms have only x-derivatives.

kdvflux=Integrate[-modtDerivative,x]

4.  The homotopy operator produces the same result.

kdvfluxHO=Expand[HomotopyOperator[-modtDerivative,x]]

The Sine-Gordon (sG) Equation

The sG equation in laboratory coordinates can be written as the system

(Subscript[u, t]) = v

(Subscript[v, t]) = (Subscript[u, xx]) + α sin u

where u = u(x, t),  v = v(x, t), and α is a constant parameter.

The sG equation has an infinite number of conservation laws.  Three are

(Subscript[D, t]) (2α cos u + (v^2) + (Subscript[u, x]^2) ) + (Subscript[D, x]) (2 (Subscript[u, x]) v) = 0

(Subscript[D, t]) (2 (Subscript[u, x]) v)+ (Subscript[D, x]) (-2α cos u + (v^2) + (Subscript[u, x]^2) )  = 0

(Subscript[D, t]) (6α (Subscript[u, x]) v cos u + (Subscript[u, x])(v^3) + (Subscript[u, x]^3) v - 8(Subscript[u, xx])(Subscript[v, x]))  
          + (Subscript[D, x]) ((3 (α^2) cos u - 3α (v^2) + α (Subscript[u, x]^2) )cos u -(1/4)((v^4) -(Subscript[u, x]^4) ) -(3/2)(Subscript[u, x]^2) (v^2) +4((Subscript[v, x]^2) + (Subscript[u, xx]^2)))= 0

A Flux for the Third sG Conservation Law

rho=-6 alpha Cos[u[x,t]] v[x,t] Derivative[1,0][u][x,t]-((v[x,t])^3) Derivative[1,0][u][x,t]-v[x,t] ((Derivative[1,0][u][x,t])^3)+8 Derivative[1,0][v][x,t] Derivative[2,0][u][x,t]

1.  Find (Subscript[D, t]) ρ.

tDerivative=D[rho,t]

2.  Replace all t-derivatives in (Subscript[D, t]) ρ using the given PDE.
      For the sG equation, (Subscript[u, t]) -> v,  (Subscript[v, t]) -> (Subscript[u, xx]) + sin u.

modtDerivative=tDerivative/.{Derivative[m_,1][u][x,t]->D[v[x,t],{x,m}],Derivative[n_,1][v][x,t]->D[D[u[x,t],{x,2}]+alpha*Sin[u[x,t]],{x,n}]}

3.  From the continuity equation, (Subscript[D, t]) ρ is the same as -(Subscript[D, x])J, so the flux can be obtained by integrating -(Subscript[D, t]) ρ with respect to x.

sGflux=Integrate[-modtDerivative,x]

4.  Replace the integration by parts with the homotopy operator.

sGfluxHO=HomotopyOperator[-modtDerivative,x]

Expand[D[sGfluxHO,x]-(-modtDerivative)]

Conservation Laws for PDEs in (2+1) Dimensions

A conservation law for a PDE in (2 + 1) dimensions satisfies

(Subscript[D, t]) ρ + ∇·J = 0    (on the PDE)

where ρ  is the conserved density,
             J = ((Subscript[J, 1]), (Subscript[J, 2])) is the corresponding flux, and
             ∇ = ((Subscript[D, x]), (Subscript[D, y]))

Zakharov-Kuznetsov (ZK) Equation

The generalized ZK equation,

(Subscript[(Subscript[u, t]) + α (u^n) (Subscript[u, x]) + β ((Subscript[u, xx]) + (Subscript[u, yy])), x]) = 0

The ZK equation has three conservation laws.  The first two are

(Subscript[D, t]) (u) + (Subscript[D, x]) ((α/(n+1))(u^(n+1)) + β (Subscript[u, xx])) + (Subscript[D, y]) (β (Subscript[u, xy])) = 0

(Subscript[D, t]) ((u^2)) + (Subscript[D, x]) (((2α)/(n+2))(u^(n+2)) - β ((Subscript[u, x]^2) - (Subscript[u, y]^2) ) + 2β (u (Subscript[u, xx])+u (Subscript[u, yy])) ) + (Subscript[D, y]) (-2β (Subscript[u, x])(Subscript[u, y])) = 0

A Flux for the Second ZK Conservation Law

rho=u[x,y,t]^2

1.  Find (Subscript[D, t]) ρ.

tDerivative=D[rho,t]

2.  Replace all t-derivatives in (Subscript[D, t]) ρ using the given PDE.
      For the ZK equation, (Subscript[u, t]) -> -(α (u^n)(Subscript[u, x]) + β (Subscript[((Subscript[u, xx]) + (Subscript[u, yy])), x])).

modtDerivative=tDerivative/.{Derivative[m_,p_,1][u][x,y,t]->-(D[alpha*u[x,y,t]^n*D[u[x,y,t],x]+beta*D[D[u[x,y,t],{x,2}]+D[u[x,y,t],{y,2}],x],{x,m},{y,p}])}

3.  Invert the divergence using the homotopy operator.

HomotopyOperator[-modtDerivative,{x,y}]

The Homotopy Operator

The homotopy operator does integration by parts on expressions in one independent variable involving unspecified functions.

If f is a derivative, then

F = (Subscript[D, x]^( -1))f = ∫f dx = (Subscript[ℋ, u(x)])f

The homotopy operator inverts divergences on expressions in two or three independent variables involoving unspecified functions.

If f is a divergence with respect to x and y, then

F = (Div^(-1))f = ( (Subscript[ℋ, u(x, y)]^((x)))f, (Subscript[ℋ, u(x, y)]^((y))) f)

On functions with one independent variable,

(Subscript[ℋ, u(x)])f = (∫_( 0)^1)(∑_(j=1)^N)((Subscript[I, (u^((j)))])f)[λu]((d λ)/λ)

N is the number of dependent variables

((Subscript[I, (u^((j)))])f)[λu] means that in ((Subscript[I, (u^((j)))])f), one replaces u⟶λ u, (Subscript[u, x])⟶λ (Subscript[u, x]), v⟶λ v, etc.

(Subscript[I, (u^((j)))])f = (∑_(k=1)^(Subscript[M, x]^((j))))((∑_(i=0)^(k-1))(((Subscript[u, ix]^((j)))(-(Subscript[D, x])))^(k-(i+1)))) ((∂f)/(∂(Subscript[u, kx]^((j)))))

Integration by Parts with the Homotopy Operator

f = 3(Subscript[u, x])(v^2)sin u - (Subscript[u, x]^3)sin u - 6v (Subscript[v, x])cos u +2(Subscript[u, x])(Subscript[u, xx])cos u + 8(Subscript[v, x])(Subscript[v, xx])
     where u = u(x) and v = v(x)

f=3*D[u[x],x]*v[x]^2*Sin[u[x]] - D[u[x],x]^3*Sin[u[x]] - 6*v[x]*D[v[x],x]*Cos[u[x]]+2*D[u[x],x]*D[u[x],{x,2}]*Cos[u[x]]+8*D[v[x],x]*D[v[x],{x,2}]

Compute the integrand :

(Subscript[I, u])f = u((∂f)/(∂(Subscript[u, x]))) + ((Subscript[u, x])I - u (Subscript[D, x]))((∂f)/(∂(Subscript[u, xx])))
         = 3u (v^2)sin u - u (Subscript[u, x]^2)sin u + 2(Subscript[u, x]^2)cos u

(Subscript[I, v])f = v((∂f)/(∂(Subscript[v, x])))+ ((Subscript[v, x])I - v (Subscript[D, x]))((∂f)/(∂(Subscript[v, xx])))
         = -6(v^2)cos u + 8(Subscript[v, x]^2)

Finish the integration :

F = (Subscript[ℋ, u(x)])f = (∫_0^1)((Subscript[I, u])f + (Subscript[I, v])f) [λu]((dλ)/λ)
    = (∫_0^1)(3(λ^2)u (v^2)sin (λ u) - (λ^2)u (Subscript[u, x]^2)sin (λ u) + 2λ (Subscript[u, x]^2)cos(λu)
-6λ (v^2)cos(λ u) + 8λ (Subscript[v, x]^2)) d λ
    = 4(Subscript[v, x]^2) -( 3(v^2) - (Subscript[u, x]^2) )cos u

Integrate[f,x]

hoResult=HomotopyOperator[f,x]

Expand[D[hoResult,x]-f]

Other Examples

(Subscript[f, 1]) = (Subscript[v, xxx])sin u + (Subscript[u, x]) (Subscript[v, xx])cos u- (Subscript[v, xxxx]),  u = u(x)  and  v = v(x)

f1=D[v[x],{x,3}]*Sin[u[x]]+D[u[x],x]*D[v[x],{x,2}]*Cos[u[x]]-D[v[x],{x,4}]

Integrate[f1,x]

HomotopyOperator[f1,x]

(Subscript[F, 1]) = (Subscript[v, xx]) sin u - (Subscript[v, xxx])

(Subscript[f, 2])  = (((Subscript[u, xx]^2) + (Subscript[u, x]) (Subscript[u, xxx]))/(2(Sqrt[(Subscript[u, x]) (Subscript[u, xx])])))

f2=(1/2)*(D[u[x],{x,2}]^2+D[u[x],x]*D[u[x],{x,3}])/Sqrt[D[u[x],x]*D[u[x],{x,2}]]

Integrate[f2,x]

HomotopyOperator[f2,x]

(Subscript[F, 2]) = (Sqrt[(Subscript[u, x]) (Subscript[u, xx])])

Inverting a Divergence with the Homotopy Operator

f = (Subscript[u, x]) (Subscript[v, y]) - (Subscript[u, xx]) (Subscript[v, y]) - (Subscript[u, y]) (Subscript[v, x]) + (Subscript[u, xy]) (Subscript[v, x])

By hand,

(OverTilde[F]) = (Div^(-1))f = (({{u (Subscript[v, y]) - (Subscript[u, x]) (Subscript[v, y])}, {-u (Subscript[v, x]) + (Subscript[u, x]) (Subscript[v, x])}}))

Using the homotopy formula,

F = (Div^(-1))f = (({{(1/2)u (Subscript[u, y ]) + (1/4)(Subscript[u, y]) (Subscript[v, x]) - (1/2)(Subscript[u, x]) (Subscript[v, y]) + (1/4)u (Subscript[u, xy]) - (1/2)(Subscript[u, y]) v + (1/2)(Subscript[u, xy]) v}, {-(1/2)u (Subscript[v, x]) - (1/4)u (Subscript[v, xx]) + (1/4)(Subscript[u, x]) (Subscript[v, x]) + (1/2)(Subscript[u, x]) v - (1/2)(Subscript[u, xx]) v}}))

If  K = (OverTilde[F]) - F,  then  K = Curl θ =  (({{(Subscript[D, y])θ}, {-(Subscript[D, x])θ}})),  with
          θ = (1/2)u v - (1/4)u (Subscript[v, x]) - (1/2)(Subscript[u, x]) v

Since Div (Curl K) = 0, the terms in K can be removed from F,  leaving (OverTilde[F]).

The multi-variable Homotopy Operator contains an algorithm for removing curl terms.

f=D[u[x,y],x]*D[v[x,y],y]-D[u[x,y],{x,2}]*D[v[x,y],y]-D[u[x,y],y]*D[v[x,y],x]+D[u[x,y],x,y]*D[v[x,y],x]

HomotopyOperator[f,{x,y}]

Other Capabilities of the Homotopy Operator

An example with three independent variables

g = ((Subscript[u, x])/v) - ((u (Subscript[v, x]))/(v^2)) + ((Subscript[v, y])/w) - ((v (Subscript[w, y]))/(w^2)) + ((Subscript[w, z])/u) - ((w (Subscript[u, z]))/(u^2))

g=((Derivative[1,0,0][u][x,y,z])/(v[x,y,z]))-(u[x,y,z] Derivative[1,0,0][v][x,y,z])/((v[x,y,z])^2)+((Derivative[0,1,0][v][x,y,z])/(w[x,y,z]))-(v[x,y,z] Derivative[0,1,0][w][x,y,z])/((w[x,y,z])^2)+((Derivative[0,0,1][w][x,y,z])/(u[x,y,z]))-(w[x,y,z] Derivative[0,0,1][u][x,y,z])/((u[x,y,z])^2)

HomotopyOperator[g,{x,y,z}]

G = (u/v) + (v/w) + (w/u)

Integration of a parameterized differential function

g = α (Subscript[u, x]) v + β u (Subscript[v, x])

g=alpha*D[u[x],x]*v[x]+beta*u[x]*D[v[x],x]

HomotopyOperator[g,x,{alpha,beta}]

Further Improvements

Reconfigure the multi-dimensional versions so that curl terms are avoided early in the calculations

Integrate non-exact expressions, for example

∫((Subscript[u, x]) v + u (Subscript[v, x]) +(u^2)(Subscript[u, xx])) dx = u v +∫(u^2) (Subscript[u, xx]) dx

References

W. Hereman, B. Deconinck,  and L. D. Poole, Continuous and discrete homotopy operators : a theoretical approach made concrete, (2006) Math. Comput. Simul  74: 352 - 360.

W. Hereman, M. Colagrosso, R. Sayers, A. Ringler, B. Deconinck, M. Nivala, and M. S. Hickman, Continuous and discrete homotopy operators and the computation of conservation laws, in Differential Equations with Symbolic Computation, D. Wang, and Z. Zheng, Eds., (2005) Birkhäuser Verlag, Basel.

P. J. Olver, Applications of Lie Groups to Differential Equations, 2nd ed., (1993) Grad.Texts in Math. 107, Springer, New York.