home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / CALCULUS.PAK / DIRACDEL.M next >
Encoding:
Text File  |  1992-07-29  |  29.3 KB  |  827 lines
  1. (* :Mathematica Version: 2.1 *)
  2.  
  3. (* :Context: Calculus` *)
  4.  
  5. (* :Title: DiracDelta *)
  6.  
  7. (* :Author: E.C.Martin *)
  8.  
  9. (* :Summary:
  10.     Defines the UnitStep function and its generalized derivative, DiracDelta.
  11.     Provides Limit, Integrate, and Derivative rules for UnitStep and
  12.     DiracDelta.
  13. *)
  14.  
  15. (* :Keywords: distribution, generalized function, Dirac, delta, Heaviside, step *)
  16.  
  17. (* :Package Version: 1.0 *)
  18.  
  19. (* :History: original version 1/92 ECM, WRI. *)
  20.  
  21.  
  22. (* :Discussion:
  23.  
  24.   The Laplace transform is frequently defined as an integral from 0- to
  25.   infinity.  The definite integral rule defined in this package for integrals
  26.   containing DiracDelta and UnitStep implements instead an integral from 0+
  27.   to infinity:
  28.  
  29.     Integrate[DiracDelta[t] Exp[-s t], {t, 0, Infinity}]
  30.  
  31.   The expected result for the most common definition of the Laplace transform
  32.   of DiracDelta is 1, but this integral returns 0.  In order to implement the
  33.   common definition of the Laplace transform, the lower limit of the integral
  34.   must be approached from the left:
  35.  
  36.     int = Integrate[DiracDelta[t] Exp[-s t], t]
  37.     
  38.     Limit[int, t->Infinity] - Limit[int, t->0, Direction->1] 
  39.  
  40.   The behavior of the definite integral rule defined in this package is
  41.   directly related to the default behavior of the built-in Limit function.
  42.  
  43. *)
  44.  
  45. (* :Sources: "Generalized Functions", R.F.Hoskins, Ellis Horwood Limited, 1979.
  46.          "Theory of Distributions, the Sequential Approach", P.Antosik,
  47.         J.Mikusinski, R.Sikorski, Elsevier Scientific, 1973.
  48. *)
  49.         
  50. (* :Warning: Should the meaning of "Direction -> Automatic" change for
  51.     the built-in Limit, then corresponding changes need to be made to
  52.     Limit rules in this package.
  53.     Definitions for D, Derivative, Limit, Integrate, and Abs are affected by
  54.     loading this package.
  55. *)
  56.  
  57.  
  58. BeginPackage["Calculus`DiracDelta`"]
  59.  
  60. Unprotect[ UnitStep, DiracDelta ]
  61.  
  62. Map[Clear, {UnitStep, DiracDelta}]
  63.  
  64. UnitStep::usage =
  65. "UnitStep[x] is a function that is 1 for x > 0 and 0 for x < 0.
  66. UnitStep[x1, x2, ...] is 1 for (x1 > 0) && (x2 > 0) && ... and 0 for
  67. (x1 < 0) || (x2 < 0) || ... ."
  68.  
  69. DiracDelta::usage =
  70. "DiracDelta[x] is a distribution that is 0 for x != 0 and satisfies
  71. Integrate[DiracDelta[x], {x, -Infinity, Infinity}] = 1.
  72. DiracDelta[x1, x2, ...] is a distribution that is 0 for x1 != 0 ||
  73. x2 != 0 || ... and satisfies Integrate[DiracDelta[x1, x2, ...],
  74. {x1, -Infinity, Infinity}, {x2, -Infinity, Infinity}, ...] = 1."
  75.  
  76. ZeroValue::usage =
  77. "ZeroValue is an option used to specify the value of UnitStep when one or more
  78. of the arguments are zero.  In particular, UnitStep[x1, x2, ..., ZeroValue -> 1]
  79. is 1 for (x1 >= 0) && (x2 >= 0) && ... and UnitStep[x1, x2, ..., ZeroValue -> 0]
  80. is 0 for (x1 <= 0) || (x2 <= 0) || ... ."
  81.  
  82. SimplifyUnitStep::usage =
  83. "SimplifyUnitStep[expr] returns a form of expr with fewer instances of UnitStep."
  84.  
  85. SimplifyDiracDelta::usage =
  86. "SimplifyDiracDelta[expr] simplifies a zero distribution in expr to 0."
  87.  
  88. Begin["`Private`"]
  89.  
  90. Attributes[UnitStep] = Attributes[DiracDelta] = Orderless
  91.  
  92. (********************************* DiracDelta *********************************)
  93.  
  94. DiracDelta[x__] := 0 /; !Apply[And, Map[(#==0)&, {x}]]
  95.  
  96. Unprotect[Derivative] 
  97. Derivative[n__Integer?NonNegative][DiracDelta][x__] :=
  98.     0 /; !Apply[And, Map[(#==0)&, {x}]]
  99. Protect[Derivative] 
  100.  
  101. (********************************* UnitStep ***********************************)
  102.  
  103. Options[UnitStep] = {ZeroValue -> 1}
  104. ZeroValue::zerov =
  105. "Warning: `` must be None or a number between 0 and 1.  ZeroValue default
  106. assumed."
  107.  
  108. UnitStep[x__, ___Rule] := 0 /; Apply[Or, Negative[{x}]]
  109.  
  110. (* Handling an illegal ZeroValue. *)
  111.  
  112. UnitStep[x__, ZeroValue -> k_] :=
  113.     (Message[ZeroValue::zerov, k];
  114.      UnitStep[x]) /; !(k === None || (NumberQ[N[k]] && 0 <= N[k] <= 1))
  115.  
  116. (* Some zero arguments. *)
  117.  
  118. UnitStep[x__?RuleFreeQ, r___Rule] :=
  119.   Module[{k = ZeroValue /. {r} /. Options[UnitStep], nk, select, comp},
  120.     (comp = Complement[{x}, select];
  121.      k^Length[select] If[comp === {}, 1, Apply[UnitStep, Join[comp, {r}]]]) /;
  122.     NumberQ[nk = N[k]] && 0 <= nk <= 1 &&
  123.     Length[select = Select[{x}, N[#]==0&]] > 0
  124.   ]
  125.  
  126. (* Some positive arguments. *)
  127.  
  128. UnitStep[x__?RuleFreeQ, r___Rule] :=
  129.   Module[{select, comp},
  130.     (comp = Complement[{x}, select];
  131.      If[comp === {}, 1, Apply[UnitStep, Join[comp, {r}]]]
  132.     ) /; Length[select = Select[{x}, N[#] > 0&]] > 0
  133.   ]
  134.  
  135. (* Power rule for ZeroValue equal to 0 or 1. *)
  136.  
  137. UnitStep /: UnitStep[f__?RuleFreeQ,
  138.                 r___Rule]^a_?((NumberQ[#] && Positive[#])&) :=
  139.   Module[{nk = N[ZeroValue /. {r} /. Options[UnitStep]]},
  140.     UnitStep[f, r] /; (NumberQ[nk] && (nk === 1. || nk === 0))
  141.   ]
  142.  
  143.  
  144. (****************************** limits, linearity ***************************)
  145.  
  146. Unprotect[Limit]
  147.  
  148. Limit[f_, x_ -> c_, dir___] :=
  149.     Module[{result},
  150.         Limit[result, x -> c, dir] /;
  151.         (result = Collect[f, Union[Cases[{f}, UnitStep[__], Infinity],
  152.             Cases[{f}, DiracDelta[__], Infinity],
  153.             Cases[{f}, Derivative[__][DiracDelta][__], Infinity]]];
  154.          !SameQ[f, result])
  155.     ] /; !(FreeQ[f, UnitStep] && FreeQ[f, DiracDelta])
  156.  
  157. Limit[a_. Derivative[n__Integer?NonNegative][DiracDelta][f__] + b_.,
  158.         x_ -> c_, dir___] :=
  159.     Module[{lima, limb, limdeltader},
  160.       If[limdeltader == 0 && !FreeQ[lima, DirectedInfinity],
  161.         limb,
  162.             lima limdeltader + limb]  /;
  163.       (FreeQ[{lima, limb, limdeltader} =
  164.        Limit[{a, b, Derivative[n][DiracDelta][f]}, x->c, dir], Limit] &&
  165.        !(limdeltader != 0 &&
  166.          !FreeQ[lima, DirectedInfinity] && !FreeQ[limb, DirectedInfinity]))
  167.     ] /; !(TrueQ[a==1] && TrueQ[b==0])
  168.  
  169. Limit[a_. DiracDelta[f__] + b_., x_ -> c_, dir___] :=
  170.     Module[{lima, limb, limdelta},
  171.       If[limdelta == 0 && !FreeQ[lima, DirectedInfinity],
  172.         limb,
  173.             lima limdelta + limb]  /;
  174.        (FreeQ[{lima, limb, limdelta} =
  175.         Limit[{a, b, DiracDelta[f]}, x->c, dir], Limit] &&
  176.         !(limdelta != 0 &&
  177.           !FreeQ[lima, DirectedInfinity] && !FreeQ[limb, DirectedInfinity]))
  178.     ] /; !(TrueQ[a==1] && TrueQ[b==0])
  179.  
  180. Limit[a_. UnitStep[f__?RuleFreeQ, zero___Rule] + b_., x_ -> c_, dir___] :=
  181.     Module[{lima, limb, limunit},
  182.         If[limunit == 0 && FreeQ[lima, Limit] &&
  183.                         !FreeQ[lima, DirectedInfinity],
  184.         limb,
  185.             lima limunit + limb]  /;
  186.           ({lima, limb} = Limit[{a, b}, x->c, dir];
  187.           FreeQ[limunit = Limit[UnitStep[f, zero], x->c, dir], Limit] &&
  188.        !(limunit != 0 && FreeQ[{lima, limb}, Limit] &&
  189.          !FreeQ[lima, DirectedInfinity] && !FreeQ[limb, DirectedInfinity]))
  190.     ] /; !(TrueQ[a==1] && TrueQ[b==0])
  191.  
  192.  
  193. Off[Sum::itform]
  194. Limit[Sum[a_. Derivative[n__Integer?NonNegative][DiracDelta][f__] + b_., i__],
  195.     x_ -> c_, dir___] := 
  196.     Module[{lima, limb, limdeltader},
  197.       If[limdeltader == 0 && !FreeQ[lima, DirectedInfinity],
  198.         Sum[limb, i],
  199.         Sum[lima limdeltader + limb, i]] /;
  200.        (FreeQ[{lima, limb, limdeltader} =
  201.         Limit[{a, b, Derivative[n][DiracDelta][f]}, x->c, dir], Limit] &&
  202.         !(limdeltader != 0 &&
  203.           !FreeQ[lima, DirectedInfinity] && !FreeQ[limb, DirectedInfinity]))
  204.         ] /; FreeQ[{i}, x]
  205.  
  206. Limit[Sum[a_. DiracDelta[f__] + b_., i__], x_ -> c_, dir___] := 
  207.     Module[{lima, limb, limdelta},
  208.       If[limdelta == 0 && !FreeQ[lima, DirectedInfinity],
  209.         Sum[limb, i],
  210.         Sum[lima limdelta + limb, i]] /;
  211.        (FreeQ[{lima, limb, limdelta} =
  212.         Limit[{a, b, DiracDelta[f]}, x->c, dir], Limit] &&
  213.         !(limdelta != 0 &&
  214.          !FreeQ[lima, DirectedInfinity] && !FreeQ[limb, DirectedInfinity]))
  215.     ] /; FreeQ[{i}, x]
  216.  
  217. Limit[Sum[a_. UnitStep[f__?RuleFreeQ, zero___Rule] + b_., i__],
  218.     x_ -> c_, dir___] :=
  219.     Module[{lima, limb, limunit},
  220.       If[limunit == 0 && FreeQ[lima, Limit] &&
  221.                              !FreeQ[lima, DirectedInfinity],
  222.          Sum[limb, i],
  223.          Sum[lima limunit + limb, i]] /;
  224.        ({lima, limb} = Limit[{a, b}, x->c, dir];
  225.         FreeQ[limunit = Limit[UnitStep[f, zero], x->c, dir], Limit] &&
  226.         !(limunit != 0 && FreeQ[{lima, limb}, Limit] &&
  227.           !FreeQ[lima, DirectedInfinity] && !FreeQ[limb, DirectedInfinity]))
  228.     ] /; FreeQ[{i}, x]
  229.  
  230. Limit[a_. sum_?((Head[#] === Sum)&) + b_., x_ -> c_, dir___] :=
  231.     Module[{lima, limb, limsum}, 
  232.       If[limsum == 0 && FreeQ[lima, Limit] &&
  233.                         !FreeQ[lima, DirectedInfinity],
  234.         limb,
  235.         lima limsum + limb] /;
  236.     ({lima, limb} = Limit[{a, b}, x->c, dir];
  237.      FreeQ[limsum = Limit[sum,  x->c, dir], Limit] &&
  238.      !(limsum != 0 && FreeQ[{lima, limb}, Limit] &&
  239.        !FreeQ[lima, DirectedInfinity] && !FreeQ[limb, DirectedInfinity]))
  240.     ] /; !(FreeQ[sum, UnitStep] && FreeQ[sum, DiracDelta])
  241. On[Sum::itform]
  242.  
  243.  
  244. (********************** limits, UnitStep' and DiracDelta **********************)
  245.  
  246. Limit[Derivative[n__Integer?NonNegative][DiracDelta][f__],
  247.     x_ -> c_, dir___Rule] :=
  248.    Module[{limit},
  249.      If[Apply[And, Map[NumberQ[N[#]]&, limit]],
  250.        0, Apply[Derivative[n][DiracDelta], limit]
  251.      ]  /; FreeQ[limit = Limit[{f}, x->c, dir], Limit]
  252.    ]
  253.  
  254. Protect[Limit]
  255.  
  256. DiracDelta /: Limit[DiracDelta[f__], x_ -> c_, dir___Rule] :=
  257.    Module[{limit},
  258.      If[Apply[And, Map[NumberQ[N[#]]&, limit]], 
  259.        0, Apply[DiracDelta, limit]
  260.      ]  /; FreeQ[limit = Limit[{f}, x->c, dir], Limit]
  261.    ]
  262.  
  263. (****************************** limits, UnitStep ******************************)
  264.  
  265. (* Default Limit Direction is Automatic. *)
  266.  
  267. UnitStep /: Limit[UnitStep[f__?RuleFreeQ, zero___Rule], x_ -> c_] :=
  268.      Module[{limit},
  269.        limit  /; FreeQ[limit = Limit[UnitStep[f, zero], x -> c,
  270.                 Direction -> Automatic], Limit]
  271.      ]
  272.  
  273. (* For each of UnitStep's arguments that approach an expression not equal to
  274.     zero, take the Limit inside of the UnitStep. *)
  275.  
  276. UnitStep /: Limit[UnitStep[f__?RuleFreeQ, zero___Rule],
  277.     x_ -> c_, dir___Rule] :=
  278.    Module[{limit = Limit[{f}, x -> c, dir], newf},
  279.     Limit[ UnitStep @@ Join[newf, {zero}], x -> c, dir]
  280.      /;  (newf = Map[If[!FreeQ[#[[2]], Limit] || TrueQ[#[[2]]==0],
  281.                 #[[1]], #[[2]]]&,
  282.               Transpose[{{f}, limit}]];
  283.      !SameQ[{f}, newf])
  284.    ] 
  285.  
  286. (* If even one of UnitStep's arguments approaches a negative number, then 
  287.      UnitStep is zero. *)
  288.  
  289. UnitStep /: Limit[UnitStep[f__?RuleFreeQ, ___Rule],
  290.     x_ -> c_, dir___Rule] :=
  291.    Module[{limit = Limit[{f}, x -> c, dir]},
  292.      0 /; Apply[Or, Negative[limit]]
  293.    ] 
  294.  
  295. (* In cases where one or more of UnitStep's arguments approach zero, check
  296.     the derivative of those arguments before deciding what the limit
  297.     should be.  There are four cases in where one or more of UnitStep's
  298.     arguments approach zero... *)
  299.  
  300. UnitStep /: Limit[UnitStep[f__?RuleFreeQ, ___Rule],
  301.     x_ -> c_, Direction -> dir_] :=
  302.   Module[{limit = Limit[{f}, x -> c, Direction -> dir]},
  303.       1 /; Apply[And, MapIndexed[(Positive[#1] ||
  304.         (#1==0 && IncreasingQ[{f}[[#2[[1]]]], x, c, dir]))&,
  305.         limit]]
  306.   ] /; (dir === Automatic || Negative[dir])
  307.  
  308. UnitStep /: Limit[UnitStep[f__?RuleFreeQ, ___Rule],
  309.     x_ -> c_, Direction -> dir_] :=
  310.    Module[{limit = Limit[{f}, x -> c, Direction -> dir]},
  311.        0 /; Apply[Or, MapIndexed[(Negative[#1] ||
  312.          (#1==0 && DecreasingQ[{f}[[#2[[1]]]], x, c, dir]))&,
  313.          limit]]
  314.    ] /; (dir === Automatic || Negative[dir])
  315.  
  316. UnitStep /: Limit[UnitStep[f__?RuleFreeQ, ___Rule], x_ -> c_,
  317.     Direction -> dir_?Positive] :=
  318.    Module[{limit = Limit[{f}, x -> c, Direction -> dir]},
  319.        1 /; Apply[And, MapIndexed[(Positive[#1] ||
  320.          (#1==0 && DecreasingQ[{f}[[#2[[1]]]], x, c, dir]))&,
  321.          limit]]
  322.   ] 
  323.  
  324. UnitStep /: Limit[UnitStep[f__?RuleFreeQ, ___Rule], x_ -> c_,
  325.     Direction -> dir_?Positive] :=
  326.    Module[{limit = Limit[{f}, x -> c, Direction -> dir]},
  327.        0 /; Apply[Or, MapIndexed[(Negative[#1] ||
  328.          (#1==0 && IncreasingQ[{f}[[#2[[1]]]], x, c, dir]))&,
  329.          limit]]
  330.    ]
  331.  
  332. (* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *)
  333.  
  334. IncreasingQ[f_, x_, c_, dir_] :=  Module[{g = D[f, x]},
  335.     If[!FreeQ[g, Derivative] || g == 0, Return[False]];
  336.     While[ TrueQ[Limit[g, x->c, Direction -> dir] == 0],
  337.         g = D[g, x];
  338.         If[!FreeQ[g, Derivative] || g == 0, Return[False]]
  339.     ];
  340.     TrueQ[Positive[Limit[g, x->c, Direction -> dir]]] ]
  341.  
  342. DecreasingQ[f_, x_, c_, dir_] :=  Module[{g = D[f, x]},
  343.     If[!FreeQ[g, Derivative] || g == 0, Return[False]];
  344.     While[ TrueQ[Limit[g, x->c, Direction -> dir] == 0],
  345.         g = D[g, x];
  346.         If[!FreeQ[g, Derivative] || g == 0, Return[False]]
  347.     ];
  348.     TrueQ[Negative[Limit[g, x->c, Direction -> dir]]] ]
  349.  
  350.  
  351. (****************************** Derivatives ***********************************)
  352.  
  353. Unprotect[Derivative]
  354.  
  355. Literal[Derivative[n__?NonNegative][UnitStep][x__?RuleFreeQ]] :=
  356.   Module[{l = Transpose[{{n}, {x}}], nonzero, zero, u, d},
  357.     nonzero = Select[l, #[[1]] != 0&]; zero = Complement[l, nonzero];
  358.     u = If[zero === {}, 1, UnitStep @@ Map[#[[2]]&, zero]];
  359.     d = If[nonzero === {}, 1,
  360.     ((Derivative @@ Map[#[[1]]&, nonzero]-1) @@ DiracDelta) @@
  361.     Map[#[[2]]&, nonzero]];
  362.     u * d
  363.   ] /; Length[{n}] == Length[{x}]
  364.  
  365. Literal[Derivative[n__][UnitStep][x__?RuleFreeQ, r_Rule]] :=
  366.   Module[{l = Transpose[{Drop[{n}, -1], {x}}], nonzero, zero, u, d},
  367.     nonzero = Select[l, #[[1]] != 0&]; zero = Complement[l, nonzero];
  368.     u = If[zero === {}, 1, UnitStep @@ Join[Map[#[[2]]&, zero], {r}]];
  369.     d = If[nonzero === {}, 1,
  370.     ((Derivative @@ Map[#[[1]]&, nonzero]-1) @@ DiracDelta) @@
  371.     Map[#[[2]]&, nonzero]];
  372.     u * d
  373.   ] /; Length[{n}]-1 == Length[{x}] && Apply[And, NonNegative[Drop[{n}, -1]]] &&
  374.         Last[{n}] == 0
  375.  
  376. Protect[Derivative]
  377. Unprotect[D]
  378.  
  379.  (* D rules for UnitStep with ZeroValue option *)
  380.  
  381. Literal[D[a_. UnitStep[f__?RuleFreeQ, r_Rule] + b_.,
  382.     x__?((Head[#] === Symbol || MatchQ[#, {_Symbol, _Integer}])&)]] :=
  383.    D[a UnitStep[f] + b, x] //. ZeroValueRules[r]
  384.  
  385.  (* D rules for Derivative[n__][UnitStep] with ZeroValue option *)
  386.  
  387. Literal[D[a_. Derivative[m__Integer?NonNegative][UnitStep][f__?RuleFreeQ,
  388.     r_Rule] + b_.,
  389.     x__?((Head[#] === Symbol || MatchQ[#, {_Symbol, _Integer}])&)]] :=
  390.    D[a Apply[Derivative, Drop[{m}, -1]][UnitStep][f] + b, x] //.
  391.         ZeroValueRules[r]
  392.  
  393. Protect[D]
  394.  
  395. ZeroValueRules[r_] := {Derivative[n__][UnitStep][g__?RuleFreeQ] :>
  396.                     Derivative[n, 0][UnitStep][g, r],
  397.             UnitStep[g__?RuleFreeQ] :> UnitStep[g, r]}
  398.  
  399. (**************************** integrals, linearity ****************************)
  400.  
  401. (* Note that Integrate[g[x] v[x] + h[x] v[x], x] ->
  402.          Integrate[(g[x] + h[x]) v[x], x] is built-in to Mathematica. *)    
  403.  
  404. Unprotect[Integrate]
  405.  
  406. Integrate[f_, x_] :=
  407.     Module[{result},
  408.         Integrate[result, x] /;
  409.         (result = Collect[f, Union[Cases[{f}, UnitStep[__], Infinity],
  410.             Cases[{f}, DiracDelta[__], Infinity],
  411.             Cases[{f}, Derivative[__][DiracDelta][__], Infinity]]];
  412.          !SameQ[f, result])
  413.     ] /; !(FreeQ[f, UnitStep] && FreeQ[f, DiracDelta])
  414.  
  415. Integrate[c1_. UnitStep[y1__?RuleFreeQ, r1___Rule] + c2_, x_] :=
  416.    Module[{i1},
  417.       SimplifyUnitStep[ i1 + Integrate[c2, x] ] /;
  418.     Head[i1 = Integrate[c1 UnitStep[y1, r1], x]] =!= Integrate ||
  419.     Apply[List, i1] =!= {c1 UnitStep[y1, r1], x} 
  420.    ]
  421.  
  422. Integrate[c1_. DiracDelta[y1__] + c2_, x_] :=
  423.    Module[{i1},
  424.       SimplifyUnitStep[ i1 + Integrate[c2, x] ] /;
  425.     Head[i1 = Integrate[c1 DiracDelta[y1], x]] =!= Integrate ||
  426.     Apply[List, i1] =!= {c1 DiracDelta[y1], x}
  427.    ]
  428.  
  429. Integrate[a_ UnitStep[y_, r___Rule], x_Symbol] :=
  430.     a Integrate[UnitStep[y, r], x] /; FreeQ[a, x]     
  431.  
  432. Integrate[a_ f_ UnitStep[y_, r___Rule], x_Symbol] :=
  433.        a Integrate[f UnitStep[y, r], x] /; FreeQ[a, x] && !FreeQ[f, x]
  434.  
  435. Off[Sum::itform]
  436. Integrate[a_. Sum[f_, i__] + b_., x_Symbol] :=
  437.     Module[{result = Module[{head, F},
  438.                 head = Head[F = Integrate[a f, x]];
  439.                 If[head =!= Integrate,
  440.                     Sum[F, i] + Integrate[b, x],
  441.                     $Failed]]},
  442.         result /; result =!= $Failed
  443.     ] /; !(FreeQ[f, DiracDelta] && FreeQ[f, UnitStep]) && FreeQ[{i}, x] &&
  444.         !FreeQ[f, x]
  445. On[Sum::itform]
  446.  
  447.  
  448. (**************************** multiple integrals ****************************)
  449.  
  450.     (* Use UnitStep to trim integration limits. *)
  451.  
  452. Integrate[f_ UnitStep[y__?RuleFreeQ, r___Rule]^(n_.) + g_., z__List] :=
  453.   Module[{x, a, b, na, nb, nc, f2, t, newz, newf = 
  454.       (f UnitStep[y, r]^n // SimplifyUnitStep) /. ToUniRules},
  455.      ( newz = 
  456.        Map[({x, a, b} = #; {na, nb} = N[{a, b}];
  457.         While[MatchQ[newf,
  458.              f1_ UnitStep[c1_. x + d1_., ___Rule]^(m_.) /;
  459.                 (NumberQ[N[c1]] && NumberQ[N[d1]] &&
  460.                  ({nc, t, f2} = {N[c1], -d1/c1, f1};
  461.                   If[nc > 0, N[t] > na, N[t] < nb]))],
  462.           If[nc > 0,
  463.         If[nb <= N[t], newf = 0,
  464.                    newf = f2; a = t; na = N[a]],
  465.         If[na >= N[t], newf = 0,
  466.                    newf = f2; b = t; nb = N[b]]]
  467.             ];
  468.             {x, a, b})&,
  469.            {z}];
  470.        Apply[Integrate, Join[{newf}, newz]] + Integrate[g, z] ) /;
  471.  
  472.        Apply[Or, Map[
  473.     ( MatchQ[#, {x1_Symbol, a1_, b1_} /;
  474.         ( {x, a, b} = {x1, a1, b1};
  475.           (NumberQ[N[a]] || Head[a] === DirectedInfinity) &&
  476.           (NumberQ[N[b]] || Head[b] === DirectedInfinity) ) ] &&
  477.       MemberQ[{y}, c_. x + d_. /;
  478.            ( NumberQ[N[c]] && NumberQ[N[d]] &&
  479.               (t = -d/c; If[N[c] > 0, N[t] > N[a], N[t] < N[b]]) ) ] )&,
  480.          {z}]]
  481.  
  482.   ]
  483.  
  484.     (* Eliminate UnitStep terms that are unity. *)
  485.  
  486. Integrate[f_ UnitStep[y__?RuleFreeQ, r___Rule]^(n_.) + g_., z__List] :=
  487.   Module[{x, a, b, unit, newunit = UnitStep[y, r]^n},
  488.     Scan[({x, a, b} = #;
  489.       If[(unit = (newunit //. {x -> a})) === (newunit //. {x -> b}),
  490.         newunit = unit])&,
  491.      {z}];
  492.     Integrate[f newunit + g, z] /;
  493.  
  494.     Apply[Or, Map[
  495.      ( MatchQ[#, {x1_Symbol, a1_, b1_} /;
  496.         ( {x, a, b} = {x1, a1, b1};
  497.           (NumberQ[N[a]] || Head[a] === DirectedInfinity) &&
  498.           (NumberQ[N[b]] || Head[b] === DirectedInfinity) ) ] &&
  499.            !FreeQ[{y}, x] &&
  500.            ((UnitStep[y, r]^n) //. {x -> a}) ===
  501.         ((UnitStep[y, r]^n) //. {x -> b}) )&,
  502.           {z}]]
  503.   ]
  504.  
  505.     (* Expand if it appears that the UnitStep will allow the
  506.         integration limits to be trimmed. *)
  507.  
  508. Integrate[f_ (UnitStep[y__?RuleFreeQ, r___Rule]^(n_.) + h_) + g_., z__List] :=
  509.   Module[{x, a, b, t},
  510.        Integrate[Expand[f (UnitStep[y, r]^n + h)] + g, z]  /;        
  511.  
  512.        Apply[Or, Map[
  513.     ( MatchQ[#, {x1_Symbol, a1_, b1_} /;
  514.         ( {x, a, b} = {x1, a1, b1};
  515.           (NumberQ[N[a]] || Head[a] === DirectedInfinity) &&
  516.           (NumberQ[N[b]] || Head[b] === DirectedInfinity) ) ] &&
  517.       MemberQ[{y}, c_. x + d_. /;
  518.            ( NumberQ[N[c]] && NumberQ[N[d]] &&
  519.               (t = -d/c; If[N[c] > 0, N[t] > N[a], N[t] < N[b]]) ) ] )&, 
  520.          {z}]]
  521.  
  522.   ]
  523.  
  524.  
  525.     (* Nest multiple integrals. *)
  526.  
  527. Integrate[f_, iseq__] :=
  528.   Module[{result =
  529.         Module[{integrand = f, scan, head},
  530.          scan = Scan[(head = Head[integrand = Integrate[integrand, #]];
  531.                   If[!(head =!= Integrate || integrand[[2]] =!= #),
  532.                   Return[$Failed]])&,
  533.              Reverse[{iseq}]];
  534.          If[scan === $Failed, $Failed, integrand]]},
  535.      result  /; result =!= $Failed
  536.   ]  /;     !(FreeQ[f, DiracDelta] && FreeQ[f, UnitStep]) && Length[{iseq}] > 1
  537.  
  538.  
  539. (***************************** definite integrals *****************************)
  540.  
  541. Integrate[f_, {x_, a_, b_}] :=        
  542.   Module[{F, lima, limb, pos},
  543.      limb - lima  /;
  544.    FreeQ[F = Integrate[f, x], Integrate] &&
  545.    FreeQ[{lima = Limit[F, x->a], limb = Limit[F, x->b]}, Limit] &&
  546.    (FreeQ[lima, DirectedInfinity] || FreeQ[limb, DirectedInfinity])
  547.   ] /; !(FreeQ[f, DiracDelta] && FreeQ[f, UnitStep]) && FreeQ[f, Integrate]
  548.  
  549.  
  550. (**************************** indefinite integrals ****************************)
  551.  
  552. (* Integral of a univariate step function with linear argument. *)
  553.  
  554. Integrate[f_. UnitStep[y_, r___Rule]^(_.), x_Symbol] :=
  555.   Module[{F, F0, a, b},  
  556.       SimplifyUnitStep[ (F - F0) UnitStep[a x + b, r] ]
  557.       /;  MatchQ[Collect[y, x], c_. x + d_. /; FreeQ[{c, d}, x]] &&
  558.      ({b, a} = CoefficientList[y, x];
  559.           DiracDeltaFreeQ[f, {{x -> -b/a}}] &&
  560.       FreeQ[F = Integrate[f, x], Integrate]) &&
  561.      (F0 = Limit[F, x -> -b/a];
  562.       If[!FreeQ[F0, Limit],
  563.          If[MemberQ[Messages[Power], Literal[Power::infy] :> $Off[_]],
  564.         F0 = F /. x -> -b/a,
  565.         Off[Power::infy];
  566.         F0 = F /. x -> -b/a;
  567.         On[Power::infy] ] ];
  568.       FreeQ[F0, DirectedInfinity] && FreeQ[F0, Indeterminate])
  569.   ] 
  570.  
  571.  
  572. (* Integral of a univariate delta function with arbitrary argument.  Delta's
  573.    located off the real line contribute 0. *)
  574.  
  575. Integrate[f_. DiracDelta[y_], x_Symbol] :=
  576.   Module[{solution, select, derivative, fselect},
  577.      Apply[Plus, Map[(UnitStep[x - (x /. #)])&, select] *
  578.             fselect / Abs[derivative] ]
  579.     /; Head[solution = Solve[Map[(# == 0)&, {y}], {x}]] == List && 
  580.       solution =!= {{}} && If[FreeQ[y, x], True, solution =!= {}] &&
  581.    (select = Select[solution, FreeQ[#, Complex]&]; (select === {} ||
  582.     Apply[And, Map[!TrueQ[# == 0]&, (derivative = D[y, x] /. select)]])) &&
  583.     DiracDeltaFreeQ[f, select] && SmoothQ[f, select] &&
  584.     (fselect = Limit[f, Flatten[select]];
  585.      If[!FreeQ[fselect, Limit],
  586.     If[MemberQ[Messages[Power], Literal[Power::infy] :> $Off[_]],
  587.        fselect = f /. select,
  588.        Off[Power::infy];
  589.        fselect = f /. select;
  590.        On[Power::infy] ] ];
  591.      FreeQ[fselect, DirectedInfinity] && FreeQ[fselect, Indeterminate])
  592.   ]
  593.  
  594.  
  595. (* Integral of a derivative of a univariate delta function with linear
  596.     argument. *)
  597.  
  598. Integrate[f_. Derivative[n_Integer?Positive][DiracDelta][y_], x_Symbol] :=
  599.   Module[{f0, a, b},
  600.       ( f0 Derivative[n-1][DiracDelta][a x + b] / a -
  601.       (1/a) Integrate[D[f,x] Derivative[n-1][DiracDelta][a x + b], x]
  602.       ) /; MatchQ[Collect[y, x], a1_. x + b1_. /;
  603.           ({a, b} = {a1, b1}; FreeQ[{a1, b1}, x])] &&
  604.            DiracDeltaFreeQ[f, {{x -> -b/a}}] && SmoothQ[f, {{x->-b/a}}] &&
  605.        (f0 = Limit[f, x-> -b/a];
  606.         If[!FreeQ[f0, Limit],
  607.            If[MemberQ[Messages[Power], Literal[Power::infy] :> $Off[_]],
  608.           f0 = f /. x -> -b/a,
  609.           Off[Power::infy];
  610.           f0 = f /. x -> -b/a;
  611.           On[Power::infy] ] ];
  612.             FreeQ[f0, DirectedInfinity] && FreeQ[f0, Indeterminate])
  613.   ]
  614.  
  615.  
  616. (* Integral of a multivariate step function with linear arguments. *)
  617.  
  618. Integrate[f_. UnitStep[yseq__?RuleFreeQ, r___Rule]^(_.), x_Symbol] :=
  619.   Module[{ylist, arg},
  620.       Apply[UnitStep, Join[Select[ylist, FreeQ[#, x]&], {r}]] *
  621.         Integrate[f Apply[UnitStep, Join[arg, {r}]], x]
  622.        /;  Length[ylist = {yseq}] > 1 &&
  623.        Length[arg = Select[ylist, (!FreeQ[#, x])&]] == 1 &&
  624.        MatchQ[Collect[arg, x], {a_. x + b_.} /; FreeQ[{a, b}, x]]
  625.   ] /; FreeQ[f, Integrate]
  626.     
  627.  
  628. (* Integral of a multivariate delta function with arbitrary arguments. *)
  629.  
  630. Integrate[f_. DiracDelta[yseq__], x_Symbol] :=
  631.   Module[{ylist, arg},
  632.       Apply[DiracDelta, Select[ylist, FreeQ[#, x]&]] *
  633.         Integrate[f Apply[DiracDelta, arg], x]
  634.        /;  Length[ylist = {yseq}] > 1 &&
  635.        Length[arg = Select[ylist, (!FreeQ[#, x])&]] == 1
  636.   ] /; FreeQ[f, Integrate] 
  637.     
  638.  
  639. (* Integral of a multivariate delta function derivative with arbitrary
  640.    arguments. *)
  641.  
  642. Integrate[f_. Derivative[n__Integer][DiracDelta][yseq__], x_Symbol] :=
  643.   Module[{ylist, pos},
  644.       Apply[Apply[Derivative, Drop[{n}, pos]] [DiracDelta], Drop[ylist, pos]] *
  645.       Integrate[f Apply[Apply[Derivative, {n}[[pos]]] [DiracDelta],
  646.         ylist[[pos]]], x]    
  647.        /; Length[ylist = {yseq}] > 1 &&
  648.       Length[pos = Flatten[Position[ylist, z_ /; !FreeQ[z, x]]]] == 1
  649.   ] /; FreeQ[f, Integrate] 
  650.  
  651. Protect[Integrate]
  652.  
  653. DiracDeltaFreeQ[f_, solution_List] :=
  654.   !Apply[Or, Map[TrueQ[# == 0]&, Flatten[
  655.    Join[Cases[{f}, DiracDelta[__], Infinity] /. DiracDelta -> List,
  656.     Cases[{f}, Derivative[__][DiracDelta][__], Infinity] /. 
  657.             Derivative[__][DiracDelta] -> List] /. solution
  658.                     ] ]]
  659.  
  660. SmoothQ[f_, solution_List] := 
  661.      !Apply[Or, Map[Module[{limleft = Limit[f, #, Direction -> 1],
  662.                         limright = Limit[f, #, Direction -> -1]},
  663.                    (FreeQ[limleft, Limit] || FreeQ[limright, Limit]) &&
  664.                limleft =!= limright]&,
  665.      solution]] 
  666.  
  667. SimplifyDiracDelta[expr_] :=
  668.   ( expr //. { f_ DiracDelta[___, a_. y_Symbol + b_.] :> 0 /;
  669.     (0 == (f /. y -> -b/a)) } ) /;
  670.   MemberQ[Messages[Power], Literal[Power::infy] :> $Off[_]]
  671.  
  672. SimplifyDiracDelta[expr_] :=
  673.   Module[{new},
  674.     Off[Power::infy];
  675.     new = expr //. { f_ DiracDelta[___, a_. y_Symbol + b_.] :> 0 /;
  676.             (0 == (f /. y -> -b/a)) };
  677.     On[Power::infy];
  678.     new
  679.   ] /; !MemberQ[Messages[Power], Literal[Power::infy] :> $Off[_]]
  680.  
  681. SimplifyDiracDelta[expr_] := expr /; FreeQ[expr, DiracDelta]
  682.     
  683. SimplifyUnitStep[expr_] :=
  684.    Module[{result = 
  685.      Module[{new = expr //. ToUniRules, original = Length[cases[expr]],
  686.         convertTomulti},
  687.     convertTomulti = original < Length[cases[new]];
  688.     new = Expand[new] //. UnitStepRules;
  689.     new = Collect[new, Union[cases[new]]];
  690.     If[convertTomulti, new = new //. ToMultiRules];
  691.     If[original < Length[cases[new]], expr, new]
  692.      ]},
  693.    result
  694.    ] /; !(FreeQ[expr, UnitStep] && FreeQ[expr, DiracDelta])
  695.  
  696. SimplifyUnitStep[expr_] := expr /; FreeQ[expr, UnitStep] &&
  697.                    FreeQ[expr, DiracDelta]
  698.  
  699. cases[expr_] := Join[Cases[{expr}, UnitStep[__], Infinity],
  700.         Cases[{expr}, DiracDelta[__], Infinity],
  701.         Cases[{expr}, Derivative[__][DiracDelta][__], Infinity]]
  702.  
  703.  
  704. UnitStepRules = {UnitStep[c1_. x_Symbol + d1_., r1___Rule] * 
  705.          UnitStep[c2_. x_ + d2_., r2___Rule] :>
  706.             Module[{t1 = N[d1/c1], t2 = N[d2/c2]},
  707.               If[N[c1] > 0 ,
  708.                 If[t1 > t2, UnitStep[c2 x + d2, r2],
  709.                         UnitStep[c1 x + d1, r1]],
  710.                     If[t1 > t2, UnitStep[c1 x + d1, r1],
  711.                         UnitStep[c2 x + d2, r2]]
  712.               ]] /; Apply[And, Map[NumberQ, N[{c1, c2, d1, d2}]]] &&
  713.                 N[d2/c2] != N[d1/c1] && Sign[c1] == Sign[c2],
  714.                  UnitStep[c1_. x_Symbol + d1_., r1___Rule] * 
  715.          UnitStep[c2_. x_ + d2_., r2___Rule] :> 0 /;
  716.                Apply[And, Map[NumberQ, N[{c1, c2, d1, d2}]]] &&
  717.                N[d2/c2] != N[d1/c1] && Sign[c1] != Sign[c2] &&
  718.                ((N[c1] > 0 && N[d1/c1] < N[d2/c2]) ||
  719.             (N[c2] > 0 && N[d2/c2] < N[d1/c1])),
  720.                  UnitStep[c1_. x_Symbol + d1_., r___Rule] *
  721.          DiracDelta[c2_. x_ + d2_.] :>
  722.             If[(N[c1] > 0 && N[d2/c2] < N[d1/c1]) ||
  723.                (N[c1] < 0 && N[d2/c2] > N[d1/c1]),
  724.                DiracDelta[c2 x + d2], 0] /;
  725.                Apply[And, Map[NumberQ, N[{c1, c2, d1, d2}]]] &&
  726.                N[d2/c2] != N[d1/c1],
  727.                  UnitStep[c1_. x_Symbol + d1_., r___Rule] *
  728.          Derivative[n_Integer?Positive][DiracDelta][c2_. x_ + d2_.] :>
  729.             If[(N[c1] > 0 && N[d2/c2] < N[d1/c1]) ||
  730.                (N[c1] < 0 && N[d2/c2] > N[d1/c1]),
  731.                Derivative[n][DiracDelta][c2 x + d2], 0] /;
  732.                Apply[And, Map[NumberQ, N[{c1, c2, d1, d2}]]] &&
  733.                N[d2/c2] != N[d1/c1],
  734.         c_. UnitStep[a_. x_ + b_., r1___Rule] +
  735.         d_. UnitStep[a_. x_ + b_., r2_Rule] :>
  736.           Module[{k1, k2, k},
  737.             ((c+d) If[SameQ[k, ZeroValue /. Options[UnitStep]],
  738.                   UnitStep[a x + b],
  739.                   UnitStep[a x + b, ZeroValue -> k]]
  740.             ) /; !SameQ[k1 = ZeroValue /. {r1} /. Options[UnitStep],
  741.                 k2 = ZeroValue /. r2] &&
  742.              !TrueQ[N[c + d /. x->-b/a] == 0] &&
  743.              NumberQ[N[k = (k1 c + k2 d)/(c+d) /. x->-b/a]]
  744.                   ] /; !SameQ[{r1}, {r2}],
  745.         c_. UnitStep[a_. x_ + b_., r1___Rule] +
  746.         d_. UnitStep[a_. x_ + b_., r2_Rule] :>
  747.           (c+d) UnitStep[a x + b] /; !SameQ[{r1}, {r2}] &&
  748.                 SameQ[ZeroValue /. {r1} /. Options[UnitStep],
  749.                       ZeroValue /. r2]
  750.         }
  751.  
  752. ToUniRules = {UnitStep[x__?RuleFreeQ, r___Rule] :>
  753.             Apply[Times, Map[UnitStep[#, r]&, {x}]] /; Length[{x}] > 1,
  754.           DiracDelta[x__] :>
  755.             Apply[Times, Map[DiracDelta, {x}]] /; Length[{x}] > 1,
  756.           Derivative[n__][DiracDelta][x__] :>
  757.         Apply[Times, Map[Derivative[#[[1]]][DiracDelta][#[[2]]]&,
  758.             Transpose[{{n}, {x}}] ]] /; Length[{x}] > 1}
  759.  
  760. ToMultiRules = {UnitStep[x__?RuleFreeQ, r___Rule] *
  761.         UnitStep[y__?RuleFreeQ, r___Rule] :>
  762.             Apply[UnitStep, Join[{x}, {y}, {r}]],
  763.         DiracDelta[x__] * DiracDelta[y__] :>
  764.             Apply[DiracDelta, Join[{x}, {y}]],
  765.         Derivative[m__][DiracDelta][x__] *
  766.         Derivative[n__][DiracDelta][y__] :>
  767.             Apply[Apply[Derivative, Join[{m}, {n}]][DiracDelta],
  768.                   Join[{x}, {y}]],
  769.         UnitStep[x_?RuleFreeQ, r___Rule]^n_Integer?Positive :>
  770.             Apply[UnitStep[#, r]&, Table[x, {n}]],
  771.         DiracDelta[x_]^n_Integer?Positive :>
  772.             Apply[DiracDelta, Table[x, {n}]],
  773.         (Derivative[m_][DiracDelta][x_])^n_Integer?Positive :>
  774.             Apply[Apply[Derivative, Table[m, {n}]][DiracDelta],
  775.                 Table[x, {n}]]}
  776.  
  777. RuleFreeQ[expr_] := FreeQ[expr, Rule]
  778.  
  779. (* zero distribution, Abs', and simplifying the sign of the arg of DiracDelta
  780.    or Derivative[_][DiracDelta] *)
  781.  
  782. DiracDelta /: (y_)^n_. DiracDelta[___, y_] := 0 /; IntegerQ[n] && n > 0
  783.  
  784. DiracDelta /: UnitStep[z__?RuleFreeQ, r___Rule] *
  785.   DiracDelta[x___, a_. y_Symbol + b_.] :=
  786.     Module[{k = UnitStep[z, r] /. y -> -b/a},
  787.     k DiracDelta[x, a y + b]
  788.     ] /; !FreeQ[{z}, y] && NumberQ[N[a]] && NumberQ[N[b]]
  789.  
  790. DiracDelta[a_ b_., y___] := DiracDelta[Abs[a] b, y] /; Negative[N[a]]
  791. DiracDelta[x_Plus, y___] :=
  792.     DiracDelta[ Apply[Plus, Map[(-#)&, Apply[List, x] ]] ] /;
  793.     Apply[And, Map[MatchQ[#, a_ b_. /; Negative[N[a]]]&, Apply[List, x] ]]
  794.  
  795. Unprotect[Derivative]
  796. Derivative /: (y_)^m_Integer?Positive *
  797.    Derivative[n__Integer?NonNegative][DiracDelta][x1___, y_, x2___] :=
  798.    0 /; Length[{n}] == Length[{x1, y, x2}] && m > {n}[[ Length[{x1}] + 1 ]]
  799.  
  800. Derivative /: y_ Derivative[1][Abs][y_] := Abs[y]
  801.  
  802. Derivative[n__Integer?NonNegative][DiracDelta][y1___, a_ b_., y2___] :=
  803.     (-1)^({n}[[ Length[{y1}] + 1 ]]) *
  804.     Derivative[n][DiracDelta][y1, Abs[a] b, y2] /; Negative[N[a]] &&
  805.         Length[{n}] == Length[{y1, a b, y2}]
  806.  
  807. Derivative[n__Integer?NonNegative][DiracDelta][y1___, x_, y2___] :=
  808.     (-1)^({n}[[ Length[{y1}] + 1 ]]) *
  809.     Apply[ Derivative[n][DiracDelta],
  810.            Join[{y1},      
  811.                     {Apply[Plus, Map[(-#)&, Apply[List, x] ]]},
  812.             {y2}] ] /;
  813.    Apply[And, Map[MatchQ[#, a_ b_. /; Negative[N[a]]]&, Apply[List, x] ]] &&
  814.    Length[{n}] == Length[{y1, a b, y2}]
  815. Protect[Derivative]
  816.  
  817. (*****************************************************************************)
  818.  
  819. Protect[UnitStep, DiracDelta]
  820.  
  821. End[]
  822.  
  823. EndPackage[]
  824.  
  825.  
  826.  
  827.