(* File PDE3HOPE.M: all pieces for 3D continuous homotopy operator in a file!*)
(* WH 10/16/2003 at home in Boulder, at 19:20 *)
(* Based on code of 10/01/2003 *)

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

(* *************** BEGIN loadMultiDHomotopyFunctions.m ********** *)

(* Time-stamp: <Thu Jun  5 14:29:17 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* loadMultiDHomotopyOperator                                           *)
(* Purpose: Loads all necessary functions for the multidimensional      *)
(*          homotopy operator                                           *)
(* Authors: Ingo Kabirschke, Frances Martin                             *)
(* Input: None                                                          *)
(* Output: None                                                         *)
(* Created: June 3, 2003 @ CSM                                          *)
(* Last Modified: June 5 @ 2:29 PM @ CSM                                *)
(************************************************************************)

(* 
(* Load functions *)
<< orderMultiD.m;
<< partialDArray.m;
<< zerothEulerMultiD.m;
<< eulerMultiD.m;
<< eulerArray.m;
<< interiorProductMD.m;
<< lambdaReplaceMD.m;
<<homotopyOperatorMD.m;
*)

(* *************** END loadMultiDHomotopyFunctions.m ********** *)

(* *************** BEGIN orderMultiD.m ********** *)

(* Time-stamp: <Wed Jun 11 12:48:21 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* orderMultiD[expression, dependVarHead,                               *)
(*             independVarList_List, independVar]                       *)
(* Purpose: To find the maximum number of derivatives of dependVarHead  *)
(*          with respect to independVar                                 *)
(* Authors: Ryan Sayers, Ingo Kabirschke, Kara Namanny                  *)
(* Input: An expression: expression,                                    *)
(*        a dependent variable head: dependVarHead,                     *)
(*         (Example: for u[1], u[2], u[3], the head would be u)         *)
(*        a list of indpendent variables: independVarList,              *)
(*        the independent variable: independVar                         *)   
(* Output: The maximum derivative of dependVarHead with respect to      *)
(*         independVar                                                  *)
(* Created: June 02, 2003,                                              *)
(* Last Modified: June 06, 2003 at 2:39 PM at CSM                       *)
(************************************************************************)

(* Debug Flags for orderMultiD*)
debugOrderMDInput = False;
debugDerivCases = False;
debugDerivList = False;
debugOrderMDResult = False;

orderMultiD[expression_, dependVarHead_, independVarList_List, 
            independVar_] := 
Module[
  {i, derivCasesList, derivRule, derivList, orderMDresult, xIndex},

  (* Debug input to orderMultiD *)
  If[debugOrderMDInput, 
    Print["Input to orderMultiD is "];
    Print[expression];
  ];

  (* Checks that the independent variable given is in the *)
  (*   independent variable list, if not it returns False *)
  (* If the variable is in the list, the xIndex is set to *)
  (*   the position of the variable in the list           *)
  If[Position[independVarList, independVar] =!= {},
     xIndex = Position[independVarList, independVar][[1, 1]],
     (* Else *)
     Print["The independent variable ", independVar, 
           " given to orderMultiD is not in the variable list: "];
     Print[independVarList];
     Return[False];
  ];

  (* Extract the cases of dependVar and its derivatives *)
  (* in the expression, and put them in a list          *)
  derivCasesList = 
    Union[
      Cases[expression, 
            Derivative[listOfDerivs_List, 0][dependVarHead[i_]]
                      [independVarList, t], 
           {0, Infinity}
      ],
      Cases[expression, dependVarHead[i_][independVarList, t], 
           {0, Infinity}
      ]
    ];

  (* Debug list of cases of dependVar and its derivatives *)
  If[debugDerivCases, 
    Print["derivCasesList is "];
    Print[derivCasesList]
  ];

  (* Turn off the error message encountered since listOfDerivs *)
  (* hasn't been defined *)
  Off[Part::partd];

  (* Define replacement rule:                                       *)
  (* Replace derivatives with their order with respect to the given *)
  (*   independent variable;                                        *)
  (* Replace dependent variable with 0 if it has no derivatives     *)
  derivRule = {Derivative[listOfDerivs_List, 0][dependVarHead[i_]]
                         [independVarList, t] -> listOfDerivs[[xIndex]],
               dependVarHead[i_][independVarList, t] -> 0};

  (* Turn error message back on *)
  On[Part::partd];

  (* Apply replacement rule *)
  (* Create list of derivative orders *)
  derivList = derivCasesList /. derivRule;

  (* Debug list of derivatives *)
  If[debugDerivList, 
    Print["The list of derivative orders for the ", independVar,
          " variable is "];
    Print[derivList];
  ];

  (* Find maximum derivative order *)
  orderMDresult = Max[derivList];

  (* Debug result of orderMultiD *)
  If[debugOrderMDResult, 
    Print["The maximum order with respect to ", xIndex, " is "];
    Print[orderMDresult];
  ];

  (* Return the order of the expression with respect to the given *)
  (* independent variable *)
  Return[orderMDresult];
];

(* Print["orderMultiD loaded successfully."]; *)

(* *************** END orderMultiD.m ********** *)

(* *************** BEGIN partialDArray.m ********** *)

(* Time-stamp: <Mon Jun  9 13:57:42 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* partialDArray[expression, dependVar,                                 *)
(*               independVarList_List, expressionOrderList_List]        *)
(* Purpose: To produce an array of partial derivatives of expression    *)
(*          with respect to derivatives of dependVar with respect to    *)
(*          independVarList variables                                   *)
(* Authors: Ryan Sayers, Frances Martin and Adam Ringler                *)
(* Input: An expression: expression,                                    *)
(*        a dependent variable: dependVar,                              *)
(*        a list of independent variables: independVarList_List,        *)
(*        upper bounds for the array: expressionOrderList_List          *)
(* Output: A multi-dimensional array of partial derivatives             *)
(*         ranging from 0 to expressionOrderList                        *)
(* Created: June 03, 2003, at CSM                                       *)
(* Last Modified: June 6, 2003 at 2:16 PM at CSM                        *)
(************************************************************************)

(* Debug Flags for the partialDArray *)
debugPartialDInput = False;
debugPartialDResult = False;

partialDArray[expression_, dependVar_, 
              independVarList_List, expressionOrderList_List] := 
Module[
  {zeroList, partialDResultArray},

  (* Debug input to partialDArray *)
  If[debugPartialDInput, 
     Print["Calculating an array of partial derivatives of the",
           " expression: "];
     Print[expression];
     Print["The dependent variable input is ", dependVar];
  ];

  (* Returns False if the lengths of the input lists are not equal,    *)
  (* because the function needs an order for each independent variable *)
  If[Length[independVarList] =!= Length[expressionOrderList],
     Print["In partialDArray: The length of independVarList ",
           "is not equal to the length of expressionOrderList."];
     Print["Aborting calculations."];
     Return[False]
  ];

  (* Returns an empty set if exprnOrderList contains negative numbers    *)
  (* because you cannot have a negative order derivative. The expression *)
  (* order list will only contain a negative value if the original       *)
  (* expression did not contain the dependent variable. Thus, returning  *)
  (* {} gives the correct partial derivative array.                      *)
  If[Min[expressionOrderList] < 0,
     Print["In partialDArray: expressionOrderList contains ",
           "at least one negative number."];
     Return[{}];
  ];

  (* Creates a list of zeroes of the same length as independVarList, for *)
  (* use as the lower limits (array origin) of the array of partials     *)
  zeroList = Array[0 &, Length[independVarList]];

  (* Creates the array of partial derivatives:                           *)
  (* ## will give the internal independent variables                     *)
  (*   that Array[] iterates over                                        *)
  (* Sequence@@Transpose[List[...]] puts the derivatives in the form     *)
  (*   that D[] needs for its derivative arguments                       *)
  (* and expressionOrderList + 1 provides the number of entries in       *)
  (*   each dimension                                                    *)
  partialDResultArray = 
    Array[D[expression, 
            D[dependVar[independVarList, t], 
              Sequence@@Transpose[List[independVarList, {##}]]]] &,
          expressionOrderList + 1, zeroList
    ];

  (* Debug result of partialDArray*)
  If[debugPartialDResult, 
     Print["The array of partial derivatives for the expression, "];
     Print[expression];
     Print["is"];
     Print[partialDResultArray]
  ];

  (* Return partial derivative array *)
  Return[Expand[partialDResultArray]];

]; (* End of partialDArray Module *)

(* Print["partialDArray loaded successfully."]; *)

(* *************** END partialDArray.m ********** *)

(* *************** BEGIN zerothEulerMultiD.m ********** *)

(* Time-stamp: <Mon Jun  9 14:44:20 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* zerothEulerMultiD[expression, dependVar, independVarList_List]       *)
(* Purpose: To find the zeroth-order Euler operator of expression       *)
(* Authors: Kara Namanny and Ryan Sayers                                *)
(* Input: An expression: expression,                                    *)
(*        a dependent variable: dependVar,                              *)
(*        a list of independent variables: independVarList              *)
(* Output: The zeroth-order Euler operator of expression with respect   *)
(*         to the dependent variable using independent variables given  *)
(*         in independVarList                                           *)
(* Created: June 5, 2003 at CSM                                         *)
(* Last Modified: June 6, 2003 at 11:44 AM at CSM                       *)
(************************************************************************)

(* Debug Flags *)
debugZEulerMDInput = False;
debugZEulerMDOrder = False;
debugZEulerMDResult = False;

zerothEulerMultiD[expression_, dependVar_, independVarList_List] :=
Module[{zeroList, orderList, partialsArray, zerothEulerMDResult},

  (* Debug input *)
  If[debugZEulerMDInput,
     Print["The expression input to zerothEulerMultiD is"];
     Print[expression];
  ];

  (* Create an array of zeroes of length equal to that of independVarList *)
  (* This is used as the order of the expression *)
  zeroList = Array[0 &, Length[independVarList]];

  (* Create a list of the orders of the independent variables in *)
  (* the expression *)
  orderList = 
    Table[orderMultiD[expression, 
                      Head[dependVar], 
                      independVarList,
                      independVarList[[i]]
          ], 
          {i, 1, Length[independVarList]}
    ]; 

  (* Debug order list *)
  If[debugZEulerMDOrder,
    Print["The order list from zerothEulerMultiD is "];
    Print[orderList];
  ];

  (* Create an array of partial derivatives of the expression *)
  (* with respect to partials of the dependent variable with   *)
  (* respect to the independent variables *)
  partialsArray = 
    partialDArray[expression, dependVar, independVarList, orderList];

  (* Calculate the zeroth-order Euler operator *)
  zerothEulerMDResult = 
    Expand[
      eulerMultiD[partialsArray, independVarList, orderList, zeroList]
    ];

  (* Debug result *)
  If[debugZEulerMDResult,
     Print["The result of applying the zeroth-order Euler operator to "];
     Print[expression];
     Print["with respect to the dependent variable ", dependVar, " is"];
     Print[zerothEulerMDResult];
  ];

  (* Return result of zeroth order euler operator *)
  Return[zerothEulerMDResult];
];

(* Print["zerothEulerMultiD loaded successfully."]; *)

(* *************** BEGIN zerothEulerMultiD.m ********** *)

(* *************** BEGIN eulerMultiD.m ********** *)

(* Time-stamp: <Wed Jun 11 12:49:03 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* eulerMultiD[partialsArray_List, independVarList_List,                *) 
(*             exprnOrderList_List, operatorOrderList_List]             *)
(* Purpose: To compute the higher-order Euler operator of an expression *)
(*          whose partial derivatives are given                         *)
(* Authors: Ingo Kabirschke, Kara Namanny, Ryan Sayers                  *)
(* Input: An array of partial derivatives of                            *)
(*          an expression: partialsArray,                               *)
(*        the list of independent variables: independVarList,           *)
(*        the list of the order of the independent variables in the     *)
(*          original expression: exprnOrderList,                        *)
(*        the order of the operator in a list: operatorOrderList        *)
(*          Example: for independent variables {x, y, z}, and           *)
(*          the operator E^{xxyyyz}, the operatorOrderList              *) 
(*          would be {2, 3, 1}                                          *)
(* Output: The operatorOrderList-order Euler operator acting on the     *)
(*         expression with partial derivatives given by partialsArray   *)
(* Created: June 04, 2003 at CSM                                        *)
(* Last Modified: June 6, 2003 at 1:52 AM at CSM                        *)
(************************************************************************)

(* Debug Flags for eulerMultiD *)
debugEulerMDInput = False;
debugEulerMDTerms = False;
debugEulerMDOutput = False;

eulerMultiD[partialsArray_List, independVarList_List,
            exprnOrderList_List, operatorOrderList_List] :=
Module[
  {i, sumVarsList, sumVar, eulerMDResult},
 
  (* Debugging input to eulerMultiD *)
  If[debugEulerMDInput,
    Print["The partial derivatives given to eulerMultiD are: "];
    Print[partialsArray]
  ];

  (* Returns 0 if exprnOrderList contains negative numbers because *)
  (* you cannot have a negative order derivative. The expression   *)
  (* order list will only contain a negative value if the original *)
  (* expression was defined as 0 or a constant. Thus, returning 0  *)
  (* gives the correct Euler operator.                             *)
  If[Min[exprnOrderList] < 0,
     Print["eulerMultiD called with negative order numbers."];
     Print["Returning 0."];
     Return[0]
  ];


  (* Creates variables to sum over: one for each independent variable *)
  sumVarsList = Array[sumVar, Length[independVarList]];

  (* Using the formula given in Ryan's paper for higher-order    *)
  (*   multi-dimensional Euler operators:                        *)
  (* We have nested sums from the order of the operator          *)
  (*   to the order of the expression;                           *)
  (* The sign of the term is computed in terms of the            *)
  (*   sum variables and the order of the expression;            *)
  (* The coefficients are computed by multiplying the individual *)
  (*   binomial coefficients (Product[Binomial[...]])            *)
  (* The total derivatives are taken in terms of the             *)
  (*   sum variables and the order of the expression;            *)
  (* The correct partial derivatives are extracted from the      *)
  (*   partialsArray by using the (sumVarsList + 1) location -   *)
  (*   1 must be added because the partialsArray starts at 1,    *)
  (*   while the operatorOrderList could start at 0.             *)
  eulerMDResult =
    Expand[ 
      Sum[
        (-1)^(Tr[sumVarsList - operatorOrderList]) * 
        Product[
          Binomial[sumVarsList[[i]], operatorOrderList[[i]]],
          {i, 1, Length[independVarList]}
        ] * 
        D[partialsArray[[Sequence@@Flatten[sumVarsList + 1]]],
          Sequence@@Transpose[List[independVarList, 
                                   sumVarsList - operatorOrderList]]
        ],
        Evaluate[
          Sequence@@Transpose[List[sumVarsList, operatorOrderList, 
                                   exprnOrderList]
                    ]
        ]
      ]
    ]; 

  (* Debugging Terms produced by applying the Multi-Dimensional *)
  (*   Euler Operator                                           *)
  If[debugEulerMDTerms,
     Print["The terms produced by eulerMultiD are: "];
     Print[ 
       Table[
         (-1)^(Tr[sumVarsList - operatorOrderList]) * 
         Product[
           Binomial[sumVarsList[[i]], operatorOrderList[[i]]],
           {i, 1, Length[independVarList]}
         ] * 
         D[partialsArray[[Sequence@@Flatten[sumVarsList + 1]]],
           Sequence@@Transpose[List[independVarList, 
                                    sumVarsList - operatorOrderList]]
         ],
         Evaluate[
           Sequence@@Transpose[List[sumVarsList, operatorOrderList, 
                                    exprnOrderList]
                     ]
         ]
       ]
     ]
  ]; 
 
  (* Debugging eulerMultiD Output *)
  If[debugEulerMDOutput,
     Print["The result of applying the ", operatorOrderList,
           " order Euler operator is "];
     Print[eulerMDResult]
  ];

  (* Return result of eulerMultiD *)
  Return[eulerMDResult];

]; (* End of Module *)

(* Print["eulerMultiD loaded successfully."]; *)

(* *************** END eulerMultiD.m ********** *)

(* *************** BEGIN eulerArray.m ********** *)

(* Time-stamp: <Fri Jun  6 14:13:04 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* eulerArray[partialsArray, independVarList_List,                      *)
(*            expressionOrderList_List]                                 *)
(* Purpose: To produce an array of the results of Euler operators       *)
(*          acting on an expression with partial derivatives given by   *)
(*          partialsArray                                               *)
(* Authors: Frances Martin, Adam Ringler, Ryan Sayers                   *)
(* Input: An array of partial derivatives: partialsArray,               *)
(*        a list of independent variables: independVarList_List,        *)
(*        upper bounds for the array: expressionOrderList_List          *)
(* Output: An array of the results of Euler operators using partial     *)
(*         derivatives given in partialsArray                           *)
(* Created: June 03, 2003                                               *)
(* Last Modified: June 6, 2003 @ 1:17 AM @ CSM                          *)
(************************************************************************)

(* Debug Flags for eulerArray *)
debugEulerArrayInput = False;
debugEulerArrayResult = False;

eulerArray[partialsArray_List, independVarList_List, 
           expressionOrderList_List] :=
Module[
  {zeroList, eulerArrayResult},

  (* Debug input to eulerArray *)
  If[debugEulerArrayInput, 
     Print["The array of partial derivatives input to eulerArray is: "];
     Print[partialsArray]
  ];

  (* Returns False, aborting calculations, if the lengths of the input *)
  (* lists are not equal, because need an order for each independent   *)
  (* variable                                                          *)
  If[Length[independVarList] =!= Length[expressionOrderList],
     Print["In eulerArray: The length of independent variable list ",
           "is not equal to the length of the expression order list."];
     Print["Aborting calculations."];
     Return[False]
  ];

  (* Returns an empty set if expressionOrderList contains negative *)
  (* numbers. The order list will only contain negative values if  *)
  (* the original expression was defined as 0 or a constant. Thus, *)
  (* returning the empty set gives the correct array of Euler      *)
  (* operators.                                                    *)
  If[Min[expressionOrderList] < 0,
     Print["In eulerArray: expressionOrderList contains ",
           "at least one negative number."];
     Print["Returning an empty set for the array of Euler operators."];
     Return[{}];
  ];

  (* Create a list of zeros of the same length as independVarList, for *)
  (* use as the lower limits (index origin) of the array.              *)
  zeroList = Array[0 &, Length[independVarList]]; 

  (* Creates an array of Euler operators using eulerMultiD[...] to     *)
  (*  calculate the individual Euler Operators                         *)
  (* {##} are the internal array variables                             *)
  eulerArrayResult = 
    Array[
      eulerMultiD[
        partialsArray, independVarList, expressionOrderList, {##}] &,
      (* Limits of the array in the form: dimensions, index origin *)
      expressionOrderList + 1, zeroList
    ];

  (* Debug eulerArray result *)
  If[debugEulerArrayResult,
    Print["The array of euler operators is: "];
    Print[eulerArrayResult]
  ];

  (* Return the array of Euler operators *)
  Return[eulerArrayResult];

];

(* Print["eulerArray loaded successfully."]; *)

(* *************** END eulerArray.m ********** *)

(* *************** BEGIN interiorProductMD.m ********** *)

(* Time-stamp: <Fri Jun  6 12:54:54 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* interiorProductMD[eulerInputArray_List, dependVar,                   *)
(*                   independVarList_List, orderList_List, independVar] *)
(* Purpose: To produce the interior product given an Euler operator     *)
(*          array of an expression                                      *)
(* Authors: Ryan Sayers and Ingo Kabirschke                             *)
(* Input: An Euler operator array of an expresion: eulerInputArray,     *)
(*        a dependent variable: dependVar,                              *)
(*        a list of independent variables: independVarList,             *)
(*        a list of the order of the independent variables in           *)
(*          the original expression: orderList,                         *)
(*        the independent variable: independVar                         *)
(* Output: The interior product of the original expression with         *)
(*         respect to the independent variable and dependent variable   *)
(* Created: June 5, 2003 at CSM                                         *)
(* Last Modified: June 06, 2003 at 12:54 PM at CSM                      *)
(************************************************************************)

(* Debug Flags *)
debugIntProdMDInput = False;
debugIntProdMDSumVar = False;
debugIntProdMDIndexList = False;
debugIntProdMDTerms = False;
debugIntProdMDResult = False;

interiorProductMD[eulerInputArray_List, dependVar_,
                  independVarList_List, orderList_List, independVar_] :=
Module[
  {independVarPosition, indVarIndexList, zeroList, sumVarList, sumVar, 
   intProdMDResult},
  
  (* Debug interiorProductMD input *)
  If[debugIntProdMDInput,
     Print["The Euler array input to interiorProductMD is "];
     Print[eulerInputArray];
     Print["Calculating the Interior Product with respect to the",
           " dependent variable, ", dependVar, ", and the independent",
           " variable, ", independVar, "."];
  ];

  (* Return 0 if orderList contains a negative number, because you     *)
  (* cannot have a negative order derivative. There will only be a     *)
  (* negative value in orderList if the original expression was 0 or a *)
  (* constant. Thus, returning 0 gives the correct interior product.   *)
  If[Min[orderList] < 0,
     Print["In interiorProductMD: orderList contains a negative number."];
     Print["Returning zero."];
     Return[0]
  ];

  (* Find the position of independVar in independVarList *)
  independVarPosition = Position[independVarList, independVar][[1, 1]];

  (* Create a list of all zeroes except at independVarPosition *)
  indVarIndexList = independVarList /. 
                      independVar -> 1 /. 
                      Thread[independVarList -> 0];

  (* Debug index list *) 
  If[debugIntProdMDIndexList,
     Print["For independent variable list "];
     Print[independVarList];
     Print["And independent variable ", independVar, ", indVarIndexList is"];
     Print[indVarIndexList];
  ];

  (* Create a list of zeroes of the same length as independVarList.    *)
  (* The zeroes are used as the lower limits of the summing variables. *)
  zeroList = Array[0 &, Length[independVarList]];

  (* Create variables to sum over - one for each independent variable. *)
  sumVarList = Array[sumVar, Length[independVarList]];

  (* Debug list of summing variables *)
  If[debugIntProdMDSumVar,
     Print["The summing variable list is ", sumVarList]
  ];

  (* Using Eq. 43 from Ryan's paper, we sum over the independent   *)
  (* variables from 0 to the order, except for the independent     *)
  (* variable which goes from 0 to the order - 1 and is identified *)
  (* by its position in the list of independent variables.         *)
  (* We calculate the fraction, (i_k + 1)/(#I + 1), in terms       *)
  (* of the summing variables.                                     *) 
  intProdMDResult =
    Expand[ 
      Sum[
        ((sumVarList[[independVarPosition]] + 1) / (Tr[sumVarList] + 1)) * 
        D[
          dependVar[independVarList, t] * 
          eulerInputArray[[Sequence@@(sumVarList + 1 + indVarIndexList)]],
          Sequence@@Transpose[List[independVarList, sumVarList]]
        ],
        (* The following lines create the list of limits for the   *)
        (* summing variables                                       *)
        Evaluate[
          Sequence@@Transpose[List[
                                   sumVarList,
                                   zeroList, 
                                   orderList - indVarIndexList
                              ]
                    ]
        ]
      ]
    ];

  (* Debug the terms in the sum *)
  If[debugIntProdMDTerms,
     Print["The terms in interiorProductMD are "];
     Print[
       Table[
         ((sumVarList[[independVarPosition]] + 1) / (Tr[sumVarList] + 1)) * 
         D[
           dependVar * 
           eulerInputArray[[Sequence@@(sumVarList + 1 + indVarIndexList)]],
           Sequence@@Transpose[List[independVarList, sumVarList]]
         ],
         Evaluate[
           Sequence@@Transpose[List[
                                    sumVarList,
                                    zeroList, 
                                    orderList - indVarIndexList
                               ]
                     ]
         ]
       ]
     ]
  ];

  (* Debug result of Interior Product in Multiple Dimensions *)
  If[debugIntProdMDResult,
     Print["The result of applying interiorProductMD using the Euler array "];
     Print[eulerInputArray];
     Print[" with respect to independent variable ", independVar];
     Print[" from independent variable list "];
     Print[independVarList];
     Print[" using dependent variable ", dependVar, " is "];
     Print[intProdMDResult];
  ];

  (* Return interior product *)
  Return[intProdMDResult];

];(* End module *) 

(* Print["interiorProductMD loaded successfully."]; *)

(* *************** END interiorProductMD.m ********** *)

(* *************** BEGIN lambdaReplaceMD.m ********** *)

(* Time-stamp: <Fri Jun  6 12:07:18 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* lambdaReplaceMD[expression_, dependVarHead_, independVarList_]       *)
(* Purpose: To replace the dependent variable with lambda * dependVar   *)
(* Authors: Kara Namanny and Ingo Kabirschke                            *)
(* 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 list: independVarList             *)
(* Output: The expression with the dependent variable replaced          *)
(*         by lambda * dependent variable                               *)
(* Created: June 5, 2003 @ 11:41 AM @ CSM                               *)
(* Last Modified: June 6, 2003 @ 11:51 AM @ CSM                         *)
(************************************************************************)

(* Debug Flags for lambdaReplaceMD *)
debugLambdaReplaceMDInput = False;
debugLambdaReplaceMDResult = False;

lambdaReplaceMD[expression_, dependVarHead_, independVarList_List] :=
  Module[{i, j, lambdaReplaceMDRule, lambdaReplaceMDResult},

    (* Debugging Lambda Replacement input *)
    If[debugLambdaReplaceMDInput, 
      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              *)
    lambdaReplaceMDRule = {
      dependVarHead[i_][independVarList, t] -> 
        lambda * dependVarHead[i][independVarList, t], 
      Derivative[i_, j_][dependVarHead[k_]][independVarList, t] -> 
        lambda * Derivative[i, j][dependVarHead[k]][independVarList, t]};

    (* Apply replacement rule to expression *)
    lambdaReplaceMDResult = Expand[expression /. lambdaReplaceMDRule];

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

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

(* Print["lambdaReplaceMD loaded successfully."]; *)

(* *************** END lambdaReplaceMD.m ********** *)

(* *************** BEGIN homotopyOperatorMD.m ********** *)

(* Time-stamp: <Mon Jun  9 14:45:46 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* homotopyOperatorMD[expression_, dependVarHead_, numDependVar_ ,      *)
(*                    independVarList_List]                             *)
(* Purpose: To calculate the Total Homotopoy Operator of an expression  *)
(*          in n dimensions                                             *)
(* Authors: Ingo Kabirschke, Frances Martin, Kara Namanny, Adam Ringler *)
(*          Ryan Sayers                                                 *)
(* 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 (dependent variables) in the vector:   *)
(*           numDependVar,                                              *) 
(*        an independent variable list: independVarList                 *)
(* Output: A list of the results of applying the homotopy operator to   *)
(*         the expression with respect to each independent variable     *)
(* Created: June 5, 2003 at CSM                                         *)
(* Last Modified: June 6, 2003, at 12:10 PM at CSM                      *)
(************************************************************************)

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

homotopyOperatorMD[expression_, dependVarHead_, numDependVar_, 
                   independVarList_List] := 
Module[
  {i, k, orderList, partialsArrayTable, eulerArrayTable, 
   fullMDIntProductTable, lambdaExpressionMDTable,
   integrandMDTable, homotopyOpResultMDTable},

  (* Debug input to the Homotopy Operator *)
  If[debugHomotopyOpInput,
     Print["Applying the Total Homotopy Operator to "];
     Print[expression];
     Print["with independent variable list, ", independVarList];
  ];

  (* Check if the expression is integrable (Zeroth Euler Operator = 0)    *)
  (* eg. If the expression is closed                                      *)
  For[i = 1, i <= numDependVar, i++,
    If[zerothEulerMultiD[expression, dependVarHead[i], independVarList] != 0,
       (* Then (if not integrable) *)
       Print["The expression input, "];
       Print[expression];
       Print[" is not integrable."];
       Print["Aborting calculations."];
       Return[False];
    ];
  ];

  (* Create a list of the order of the independent variables in the  *)
  (* expression. Example: for u_xxxy + u_zz, the orderList would be  *)
  (*                      {3, 1, 2}                                  *)
  orderList = 
    Table[orderMultiD[expression, dependVarHead, independVarList, 
                      independVarList[[i]]],
          {i, 1, Length[independVarList]}
    ];
                         
  (*******************************************************************)
  (*                                                                 *)
  (* Steps below taken from Ryan Sayers' Paper, "The Total Homotopy  *)
  (* Method of Integration," draft 5/21/03                           *)
  (*                                                                 *)
  (*******************************************************************)

  (* Create a table of partial derivative arrays. This allows for    *)
  (* one-time calculation of the partial derivatives, which will be  *)
  (* used in the calculation of each of the Euler operators.         *)
  partialsArrayTable = 
    Table[partialDArray[expression, dependVarHead[i], independVarList,
                        orderList],
          {i, 1, numDependVar}
  ];

  (* Create a table of Euler Operators arrays, allowing for one-time   *)
  (* calculation of the Euler Operators, which will be used in the     *)
  (* calculation of each of the interior products.                     *)
  eulerArrayTable = 
    Table[eulerArray[partialsArrayTable[[i]], independVarList,
                     orderList],
          {i, 1, numDependVar}
    ];

  (* Calculate the interior product (called the integrand in Ryan's    *) 
  (* paper) of the given expression for each independent variable,     *)
  (* summing over all of the dependent variables.                      *)
  fullMDIntProductTable = 
    Table[Sum[     
           interiorProductMD[eulerArrayTable[[i]], dependVarHead[i], 
                             independVarList, orderList, 
                             independVarList[[k]]],
              {i, 1, numDependVar}
          ],
          {k, 1, Length[independVarList]}
    ];

  (* Replace the dependent variable with lambda * dependent variable for   *)
  (* each expression in the table of interior products.                    *)
  lambdaExpressionMDTable = 
    lambdaReplaceMD[fullMDIntProductTable, dependVarHead, independVarList];

  (* Divide each expression in the table by lambda *)
  integrandMDTable = Expand[lambdaExpressionMDTable / lambda];

  (* Integrate each expression in the table with respect to lambda to      *)
  (* find the flux (the anti-divergence). This is given in tabular form    *)
  (* sorted by independent variable.                                       *)
  homotopyOpResultMDTable = 
    Expand[Integrate[integrandMDTable, {lambda, 0, 1}]];

  (* Debug the result of homotopyOperatorMD *) 
  If[debugHomotopyOpResult,
     Print["The result of applying the Total Homotopy Operator is "];
     Print[homotopyOpResultMDTable];
  ];

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

]; (* End of Module *)

(* Print["homotopyOperatorMD loaded successfully."]; *)

Print["Multi-Dimensional Continuous Homotopy Operator loaded successfully."];

(* *************** END homotopyOperatorMD.m ********** *)

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