(* File PDEHOPER.M: all pieces for continuous homotopy operator in one file! *)
(* WH 10/04/2003 *)
(* Based on code of 10/01/2003 *)

(* ************************* BEGIN of all ***************************** *)

(* ********************* BEGIN loadHomotopyFunctions.m ********** *)

(************************************************************************)
(* loadHomotopyFunctions                                                *)
(* Purpose: To load the common functions used for finding the conserved *)
(*          fluxes of systems of PDEs                                   *)
(* Authors: Ingo Kabirschke, Kara Namanny, Frances Martin               *)
(* Input: None                                                          *)
(* Output: None                                                         *)
(* Created: May 21, 2003                                                *)
(* Last Modified: May 21, 2003 @ 4:59 PM @ MLRC                         *)
(************************************************************************)

(*
<<eulerOperator.m
<<interiorProduct.m
<<orderOfExpression.m
<<lambdaReplace.m
<<homotopyOperator.m
*)

(* ********************* END loadHomotopyFunctions.m ********** *)

(* ********************* BEGIN eulerOperator.m ************************ *)

(* Time-stamp: <Fri May 30 11:50:43 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* eulerOperator[expression_, orderOfOperator_, dependVar_,             *)
(*               independVar_]                                          *)
(* Purpose: To calculate a specified order Euler operator of a given    *)
(*          expression with respect to the dependent variable           *)
(* Authors: Ingo Kabirschke, Kara Namanny, Frances Martin, Adam Ringler *)
(* Input: An expression: expression,                                    *)
(*        the order of the desired operator: orderOfOperator,           *)
(*        a dependent variable: dependVar,                              *)
(*        an independent variable: independVar                          *)
(* Output: A new expression that is the result of applying the specified*)
(*         order Euler Operator to the original expression.             *)
(* Created: May 19, 2003 at CSM                                         *)
(* Last Modified: May 30, 2003 at 11:30 AM at CSM                       *)
(************************************************************************)

(* Debug Flags for the Euler Operator *)
debugEulerOperatorInput = False;
debugEulerTerms = False;
debugEulerResult = False;

eulerOperator[expression_, orderOfOperator_, dependVar_, independVar_] :=
  Module[{i, orderOfExp, eulerResult},

    (* Debugging input to Euler Operator *)
    If[debugEulerOperatorInput, 
      Print["Applying the ", orderOfOperator, " order Euler Operator to "]; 
      Print[expression];
    ];

    (* Finding the highest order derivative in the expression. This    *)
    (* will be used to limit the sum in the Euler operator.            *)
    orderOfExp = orderOfExpression[expression, dependVar, independVar];
                                    
    (* Calculating specified euler operator (variational derivative)   *)
    (* using \cal{L}_u ^{(i)} on page 16 of Hereman's notes            *)
    eulerResult = Expand[
                    Sum[((-1)^(i - orderOfOperator)) * 
                        Binomial[i, orderOfOperator] * 
                        D[ D[expression, 
                             {Derivative[i, 0][dependVar][independVar, t], 1}],
                          {independVar, (i - orderOfOperator)}], 
                       {i, orderOfOperator, orderOfExp}]
                  ];

    (* Debugging terms: Prints terms, produced by applying the Euler   *)
    (* operator, in tabular form.                                      *)
    If[debugEulerTerms,
      Print["The terms from applying the Euler operator are: ",
            Expand[
              Table[((-1)^(i - orderOfOperator)) * 
                      Binomial[i, orderOfOperator] * 
                      D[ D[expression, 
                           {Derivative[i, 0][dependVar][independVar, t], 1}],
                        {independVar, (i - orderOfOperator)}], 
                    {i, orderOfOperator, orderOfExp}]
            ]
      ]
    ];
    
    (* Debugging result of applying the Euler Operator *)
    If[debugEulerResult,
      Print["The result of applying the ", orderOfOperator, " order
            Euler Operator to "];
      Print[expression];
      Print[" is: "];
      Print[eulerResult];
    ];
 
    (* Returning Euler Operator of specified order *)
    Return[eulerResult];
  ];

(* Print["Euler Operator loaded successfully."]; *)                 

(* ********************* END eulerOperator.m ************************ *)

(* ********************* BEGIN interiorProduct.m ************************ *)

(* Time-stamp: <Tue Jun 10 11:57:36 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* interiorProduct[expression_, dependVar_, independVar_]               *)
(* Purpose: To calculate the interior product of an expression          *)
(* Authors: Ingo Kabirschke, Kara Namanny, Adam Ringler                 *)
(* Input: An expression: expression,                                    *)
(*        a dependent variable: dependVar,                              *)
(*        an independent variable: independVar                          *)
(* Output: The interior product of the expression, using definitions    *)
(*         from Hereman's notes, page 17                                *)
(* Created: May 21, 2003 at CSM                                         *)
(* Last Modified: May 30, 2003 at 11:50 AM at CSM                       *)
(************************************************************************)

(* Debug Flags for the Interior Product *)
debugInteriorProductInput = False;
debugInteriorProductTerms = False;
debugInteriorProductResult = False;

interiorProduct[expression_, dependVar_, independVar_] := 
  Module[{i, interiorProductResult},

    (* Debugging input of interiorProduct *)
    If[debugInteriorProductInput, 
       Print["Interior Product acting on "];
       Print[expression];
       Print[" with dependent variable ", dependVar, 
             " with independent variable ", independVar];
    ];

    (* Calculates the interior product (little j) using the definition *)
    (* on page 17 of Hereman's notes                                   *)
    interiorProductResult = 
      Expand[
        Sum[D[dependVar[independVar, t] *
              eulerOperator[expression, i + 1, dependVar,
                            independVar], {independVar, i}],
          {i, 0, orderOfExpression[expression, dependVar, independVar] - 1}
        ]
      ];

    (* Debugging terms - prints terms of interior product in table form *)
    If[debugInteriorProductTerms, 
       Print["The terms of the interior product are ", 
         Expand[Table[D[dependVar[independVar, t] *
                        eulerOperator[expression, i + 1, dependVar,
                                      independVar], 
                        {independVar, i}],
                  {i, 0, orderOfExpression[expression, 
                                           dependVar, independVar] - 1}
                ] 
         ]
       ]
    ];

    (* Debugging result *)
    If[debugInteriorProductResult,
       Print["The Interior Product of "];
       Print[expression];
       Print[" with respect to ", dependVar, " is: "];
       Print[interiorProductResult];
    ];

    (* Returns interior product of given expression *)
    Return[interiorProductResult]
]; 

(* Print["Interior Product loaded successfully."]; *)

(* ********************* END interiorProduct.m ************************ *)

(* ********************* BEGIN orderOfExpression.m ******************** *)

(* Time-stamp: <Mon Sep 22 21:55:54 MDT 2003>  *)
(******************************************************************)
(* orderOfExpression[expression_, dependVar_, independVar_]       *)
(* Purpose: To determine the maximum derivative in expression of  *)
(*          dependVar with respect to independVar                 *)
(* Authors: Ryan Sayers, Frances Martin, Ingo Kabirschke,         *)
(*          Kara Namanny, Adam Ringler                            *)
(* Input: An expression: expression,                              *)
(*        a dependent variable: dependVar,                        *)
(*        and an independent variable: independVar                *)
(* Output: The maximum derivative in expression of dependVar      *)
(*         with respect to independVar                            *)
(* Created: May 15, 2003, 2:39 PM at CSM                          *)
(* Last Modified: Sep 22, 2003, 9:53 AM by Mike C.                *)
(******************************************************************)

(* Debug Flags for the Order of the Expression *)
debugOrderList = False;
debugOrderCasesList = False;
debugOrderInput = False;

orderOfExpression[expression_, dependVar_, independVar_] := 
  Module[{xDerivRule, orderCasesList, orderList, orderResult, order},

    (* Debugging input of orderOfExpression *)
    If[debugOrderInput, 
       Print["Input expression is "];
       Print[expression];
    ];

    (* Replaces derivatives of dependVar with their order         *)
    (* Replaces dependVar with 0, since it has 0 derivatives      *)
    (* Mike C.: *)
    (* Changed dummy variable "i" to "order."  "i" is too risky   *)
    xDerivRule = {Derivative[order_, 0][dependVar][independVar, t] -> order,
                  dependVar[independVar, t]-> 0};

    (* Cases will find every occurance of dependVar and its       *)
    (* derivatives in the expression and put them in a list       *)
    orderCasesList = Union[
                       Cases[expression, 
                             Derivative[order_, 0][dependVar][independVar, t], 
                         {0, Infinity}],
                       Cases[expression, 
                             dependVar[independVar, t], 
                         {0, Infinity}]
                     ];

    (* Debugging list of occurances of the dependent variable *)
    If[debugOrderCasesList, 
       Print["Cases are ", orderCasesList]
    ];

    (* Apply replacement rule *)
    orderList = orderCasesList /. xDerivRule;

    (* Debug list of orders *)
    If[debugOrderList, 
       Print["Order List is ", orderList]
    ];

    (* Find the maximum order in the list *)
    orderResult = Max[orderList];

    (* Return the highest order of the expression *)
    Return[orderResult];
  ];

(* Print["Order of Expression loaded successfully."]; *)

(* ********************* END orderOfExpression.m ******************** *)

(* ********************* BEGIN lambdaReplace.m ******************** *)

(* Time-stamp: <Fri May 30 11:23:55 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* lambdaReplace[expression_, dependVarHead_, independVar_]             *)
(* Purpose: To replace the dependent variable with lambda * dependVar   *)
(* Authors: Kara Namanny and Adam Ringler                               *)
(* Input: An expression: expression,                                    *) 
(*        a dependent variable head: dependVarHead,                     *)
(*          (For example, for u[1], u[2], and u[3] the head would be u) *)
(*        and an independent variable: independVar                      *)
(* Output: The expression with the dependent variable replaced          *)
(*         by lambda * dependent variable                               *)
(* Created: May 21, 2003 @ 3:41 PM @ CSM                                *)
(* Last Modified: May 30, 2003 @ 11:19 AM @ MLRC                        *)
(************************************************************************)

(* Debug Flags *)
debugLambdaReplaceInput = False;
debugLambdaReplaceResult = False;

lambdaReplace[expression_, dependVarHead_, independVar_] :=
  Module[{i, j, lambdaReplaceRule, lambdaReplaceResult},

    (* Debugging Lambda Replacement input *)
    If[debugLambdaReplaceInput, 
      Print["Replacing the dependent variables, ", 
        dependVarHead, "[i], in the expression, "];
      Print[expression]; 
      Print[", with lambda * ", dependVarHead, "[i]"];
    ];

    (* Defining replacement rules: multiplies every occurance of the *)
    (* dependent variable and its derivatives by lambda              *)
    lambdaReplaceRule = {
      dependVarHead[i_][independVar, t] -> 
        lambda * dependVarHead[i][independVar, t], 
      Derivative[i_, j_][dependVarHead[k_]][independVar, t] -> 
        lambda * Derivative[i, j][dependVarHead[k]][independVar, t]};

    (* Apply replacement rule to expression *)
    lambdaReplaceResult = Expand[expression /. lambdaReplaceRule];

    (* Debugging result of Lambda Replacement *)
    If[debugLambdaReplaceResult, 
      Print["The result of replacement is: "];
      Print[lambdaReplaceResult];
    ];

    (* Return expression with lambda replacement *)
    Return[lambdaReplaceResult];
  ];

(* Print["Lambda Replace loaded successfully."]; *)

(* ********************* END lambdaReplace.m ******************** *)

(* ********************* BEGIN continuousHomotopyOperator.m ************ *)

(* Time-stamp: <Tue Jul  1 03:34:48 MDT 2003>        *)
(************************************************************************)
(* continuousHomotopyOperator[expression_,dependVarHead_,numDependVar_, *)
(*                  independVar_]                                       *)
(* Purpose: To calculate the continuous Homotopoy Operator of an  *)
(*          expression                                                  *)
(* Authors: Ingo Kabirschke, Frances Martin                             *)
(* Input: An expression: expression,                                    *) 
(*        a dependent variable head: dependVarHead,                     *)
(*          (For example, for u[1], u[2], and u[3] the head would be u) *)
(*        the number of elements in the vector: numDependVar,           *) 
(*        an independent variable: independVar                          *)
(* Output: The continuous Homotopy Operator of expression               *)
(* Created: May 21, 2003 at CSM                                         *)
(* Last Modified: July 1, 2003, at 3:30 AM by Mike C.                   *)
(************************************************************************)

(* Debug Flags *)
debugHomotopyOpInput = True;
debugHomotopyOpResult = True;

continuousHomotopyOperator[expression_, dependVarHead_, numDependVar_, 
  independVar_] := 
  Module[{i, fullIntProduct, lambdaExpression, homotopyOpResult, integrand},

    (* Debugging input *)
    If[debugHomotopyOpInput,
       Print["Applying the continuous Homotopy Operator to "];
       Print[expression];
    ];

    (* Checking for an integrable expression (Zeroth Euler Operator = 0) *)
    For[i = 1, i <= numDependVar, i++,
      If[eulerOperator[expression, 0, dependVarHead[i], independVar] != 0,
         (* Then (if not integrable) *)
         Print["The expression input, "];
         Print[expression];
         Print[" is not integrable."];
         Print["Aborting calculations."];
         Return[False];
      ];
    ];

    (* Steps below taken from Hereman's Notes, pg 17 *)
    (*************************************************)

    (* Calculate the interior product (little j) of the given expression *)
    fullIntProduct = Sum[     
                       interiorProduct[expression, dependVarHead[i], 
                                       independVar],
                       {i, 1, numDependVar}
                     ];

    (* Replace dependent variable with lambda * dependent variable *)
    lambdaExpression = 
      lambdaReplace[fullIntProduct, dependVarHead, independVar];

    (* Divide by lambda *)
    integrand = Expand[lambdaExpression / lambda];

    (* Integrate with respect to lambda to find flux (big J) *)
    homotopyOpResult = Expand[Integrate[integrand, {lambda, 0, 1}]];

    (* 06/11/03 Mike C: *)
    homotopyOpResult = TrigExpand[homotopyOpResult];


    (* Debugging result *) 
    If[debugHomotopyOpResult,
       Print["The result of applying the continuous Homotopy Operator is "];
       Print[homotopyOpResult];
    ];

    (* Return final result of applying the PDE Homotopy Operator *)
    Return[homotopyOpResult];
];

(* Print["Homotopy Operator loaded successfully."]; *)

(* ********************* END continuousHomotopyOperator.m **************** *)

(* ********************* BEGIN pdeform.m ******************** *)

(*****************************************************************************)
(* pdeform[]: takes an expression and prints it with subscripted derivatives *)
(*            and removes the functional dependence (i.e. f[x,t] -> f)       *)
(*            Written by Tracy Otto & Tony Miller (CSM-1995).                *)
(*            Modified by U. Goktas.                                         *)
(*****************************************************************************)
pdeform[expres_] := expres /. {
        Derivative[n__][u[k_]][x__] :>
                SequenceForm[u, Subscript[k],Subscript[","],
                        Subscript[
                         SequenceForm @@ Flatten[Table[#[[1]], {#[[2]]}]& /@
                                        Transpose[{{x}, {n}}]]]],
        u[n_][x__] :> SequenceForm[u,Subscript[n]] }; (* end pdeform *)

(* ********************* END pdeform.m ******************** *)

Print["Continuous Homotopy Operator code loaded successfully."]; 

(* ********************************** end of all ************************ *)