This is InvLaPlace.m in view mode; [Download] [Up]
(* :Title: Inverse Laplace Transform *)
(* :Authors: Wally McClure, Brian Evans, James McClellan *)
(*
:Summary: Multidimensional, bilateral inverse Laplace transform
for signal processing expressions.
*)
(* :Context: SignalProcessing`Analog`InvLaPlace` *)
(* :PackageVersion: 2.6 *)
(*
:Copyright: Copyright 1990-1991 by Brian L. Evans
Georgia Tech Research Corporation
Permission to use, copy, modify, and distribute this software
and its documentation for any purpose and without fee is
hereby granted, provided that the above copyright notice
appear in all copies and that both that copyright notice and
this permission notice appear in supporting documentation,
and that the name of the Georgia Tech Research Corporation,
Georgia Tech, or Georgia Institute of Technology not be used
in advertising or publicity pertaining to distribution of the
software without specific, written prior permission. Georgia
Tech makes no representations about the suitability of this
software for any purpose. It is provided "as is" without
express or implied warranty.
*)
(*
:History: date started -- April 25, 1990 (adapted from "InvZTransform.m")
added Dialogue ability -- September, 1990
handles real-valued polynomials -- March, 1991
*)
(* :Keywords: Laplace transform, region of convergence *)
(* :Source: Operational Mathematics by Churchill *)
(*
:Warning: Works only on one dimensional and separable
multidimensional transforms.
*)
(* :Mathematica Version: 1.2 or 2.0 *)
(* :Limitation: *)
(*
:Discussion: 1-D rules base has 70 rules:
I. rational transforms 26 rules *
II. non-rational transforms 28 rules
III. transform properties 10 rules *
IV. transform SP structures 0 rules
V. transform strategies 6 rules
* a partial fractions strategy rule exists in section
Unlike some symbolic transform packages, these rules are maintained
in a list called InvLaPlaceRules. This was necessary because
Mathematica reordered these rules if they were coded as a set of
transformation functions. Another benefit of the list form is that
it allowed more control over how rules are applied.
The InvLaPlace object (function) first calls the interface rules
(InvLaPlaceInterfaceRules), which handle default arguments, error
messages, and multidimensional transforms. These rules simply call
MyInvLaPlace once per dimension of the transform. This inverse
transform is biased toward causal systems. However, by specifying
the proper region of convergence (ROC), anti-causal functions can
be obtained (see InvLaPlace::usage).
The driver for the one-dimensional rule base is the six argument
version of myinvlaplace. myinvlaplace[X, s, t, rm, rp, st, op]
either returns the inverse transform as a function of time (t) or
an approximation to the actual inverse transform. If the rule base
can do neither, which should not happen, the interface function
cleanup will report any errors. Arguments of myinvlaplace:
X laplace-transform function rm Rminus of ROC
s laplace-transform variable rp Rplus of ROC
t "time"-domain variable st local state semaphores
At each step in the transform rule base, the current expression
has a local state associated with it. This state consists of a
list of six boolean values. Each boolean value is associated with
a strategy. If an element is True, then that strategy has not been
tried yet; if False, then that strategy has already been tried, and
it will not be tried again. Thus, local state is used to prevent
infinite loops which would be caused by the repetitive application
of certain strategy rules. See the section below called S T A T E
D E F I N I T I O N and see section VI of the rules.
*)
(* :Functions: InvLaPlace *)
If [ TrueQ[ $VersionNumber >= 2.0 ],
Off[ General::spell ];
Off[ General::spell1 ] ];
(* B E G I N P A C K A G E *)
BeginPackage [ "SignalProcessing`Analog`InvLaPlace`",
"SignalProcessing`Analog`LSupport`",
"SignalProcessing`Support`TransSupport`",
"SignalProcessing`Support`SigProc`",
"SignalProcessing`Support`ROC`",
"SignalProcessing`Support`SupCode`",
"SignalProcessing`Support`FilterSupport`" ]
(* B E G I N P A C K A G E *)
(* U S A G E I N F O R M A T I O N *)
InvLaPlace::usage =
"InvLaPlace[f, s] and InvLaPlace[f, s, t] gives the multidimensional \
bilateral inverse Laplace transform of f. \
A region of convergence can be specified by using \
InvLaPlace[{f, rm, rp}, s, t], where rm is R- and rp is R+ \
in the region (strip) of convergence: R- < Re(s) < R+. \
Note that InvLaPlaceTransform is an alias for InvLaPlace."
(* E N D U S A G E I N F O R M A T I O N *)
Begin[ "`Private`" ]
(* U S E R I N T E R F A C E *)
(* InvL operator *)
Unprotect[InvL]
InvL/: TheFunction[ InvL[ s_, t_ ][ f_ ] ] := InvLaPlace[ f, s, t ]
Protect[InvL]
(* InvLaPlace *)
Options[InvLaPlace] :=
{ Apart -> Rational, Definition -> False,
Dialogue -> True, Simplify -> True,
Terms -> 10, TransformLookup -> {} }
(* bridge from outside world to interface to rule base *)
InvLaPlace[ x_ ] :=
invlaplacedriver[ Options[InvLaPlace], x ]
InvLaPlace[ x_, s_ ] :=
invlaplacedriver[ Options[InvLaPlace], x, s ]
InvLaPlace[ x_, s_, t_ ] :=
invlaplacedriver[ Options[InvLaPlace], x, s, t ]
InvLaPlace[ x_, s_, t_, op__ ] :=
invlaplacedriver[ ToList[op] ~Join~ Options[InvLaPlace], x, s, t ]
(* driver for interface to rule base *)
invlaplacedriver[ op_, rest__ ] :=
Block [ {},
tlist = {}; (* global variable *)
cleanup[ invlaplace[op, rest],
Replace[Dialogue, op],
TrueQ[Replace[Simplify, op]] ] ]
(* interface to rule base *)
invlaplace[ op_, args__ ] :=
Block [ {},
Replace[ linverse[op, args], InvLaPlaceInterfaceRules ] ]
InvLaPlaceInterfaceRules = {
linverse[ op_, e_ ] :>
invlaplace[ op, e, LVariables[e] ] /; LTransformQ[e],
linverse[ op_, e_ ] :>
Message[ InvLaPlace::novariables ],
linverse[ op_, e_, s_ ] :>
invlaplace[ op, e, s, DummyVariables[Length[s], Global`t] ] /;
validvarQ[s],
linverse[ op_, e_, s_, rest___ ] :>
Message[ Transform::badvar, "s", InvLaPlace, s ] /;
! validvarQ[s],
linverse[ op_, e_, s_, t_, rest___ ] :>
Message[ Transform::badvar, "time", InvLaPlace, t ] /;
! validvarQ[t],
linverse[ op_, e_List, s_, t_ ] :>
invlaplace[ op, MakeLObject[e, s], s, t ],
linverse[ op_, e_, s_Symbol, t_ ] :>
MyInvLaPlace[ e, s, t, op ],
linverse[ op_, e_, s_List, t_ ] :>
multiDInvLaPlace[ e, s, t, op ]
}
cleanup[ trans_, dialogue_, simplify_ ] :=
Block [ {adjtrans, lt, sdialogue},
sdialogue = SameQ[dialogue, All];
lt = If [ simplify,
SPSimplify[trans,
Dialogue -> sdialogue,
Variables -> tlist],
trans ];
If [ InformUserQ[dialogue], Scan[explain, lt, Infinity] ];
lt ]
explain[ invlaplace[ f_, s_, rest__ ] ]:=
Message[ Transform::incomplete, "inverse Laplace transform", f, s ]
validvarQ[ x_Symbol ] := True
validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
validvarQ[ x_ ] := False
(* M E S S A G E S *)
InvLaPlace::badROC =
"Improper region of convergence in ``."
InvLaPlace::badterms =
"Option Terms is not an integer (Terms -> ``)."
InvLaPlace::novariables =
"No Laplace variables specified in inverse Laplace transform."
(* S T A T E D E F I N I T I O N *)
Lfactorfield = 1 (* factor denominator *)
Lexpandfield = 2 (* apply Expand[] to expression *)
Lexpandallfield = 3 (* apply ExpandAll[] to expression *)
Lthefunfield = 4 (* apply TheFunction to expression *)
Lapartfield = 5 (* apply Apart[] to expression *)
Lmyapartfield = 6 (* apply MyApart[] to rational poly. *)
Lnormalizefield = 7 (* normalize denominator *)
Lseriesfield = 8 (* approximate function with power series *)
Lsplitfield = 9 (* split rat. fun. into rat. poly & other *)
Ldefinitionfield = 10 (* apply defintion of inv. Laplace trans. *)
Lstatevariables = 10
initInvLstate[] := Table[True, {Lstatevariables}]
nullInvLstate[] := Table[False, {Lstatevariables}]
(* S U P P O R T I N G R O U T I N E S *)
(* Causal or anti-causal *)
leftsided[ p_, t_ ] :=
( p /. { CStep[ b_. t + c_. ] -> CStep[ -b t + c ], t -> - t } )
rightsided[ p_, t_ ] := p
leftOrRightSided[ a_, p_, t_, rp_ ] := rightsided[ p, t ] /; InfinityQ[rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ] /; SameQ[Re[a], rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ] /; SameQ[-Re[a], rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ] /; N[Re[a] > rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := rightsided[ p, t ]
(* MyInvLaPlace *)
MyInvLaPlace[ e_, s_, t_, op_ ] :=
Block [ {rminus, rplus, valid},
rminus = SPSimplify[ GetRMinus[e] ];
rplus = SPSimplify[ GetRPlus[e] ];
valid = Apply[And, Thread[ToList[rminus] < ToList[rplus]]];
If [ TrueQ[! valid],
Message[ InvLaPlace::badROC, e ],
MyInvLaPlace[ TheFunction[e], s, t, rminus,
rplus, initInvLstate[], op ] ] ] /;
LTransformQ[e]
MyInvLaPlace[ e_, s_, t_, op_ ] :=
MyInvLaPlace[ e, s, t, -Infinity, Infinity, initInvLstate[], op ]
MyInvLaPlace[ e_, s_, t_, rm_, rp_, st_, op_ ] :=
Block [ {laste, newe, newinvrules, norme, trace},
(* add the time variable to the global list of time vars *)
AppendTo[tlist, t];
(* generate the rules *)
newinvrules = TransformFixUp[ s, t, op, myinvlaplace, False,
InvLaPlace, Null, Null ];
InvLaPlaceRules = Join[ newinvrules, OriginalInvLaPlaceRules ];
(* determine the level of dialogue and pre-process f(t) *)
trace = SameQ[ Replace[Dialogue, op], All ];
norme = normalizeExpression[e, s];
If [ trace && ! SameQ[norme, e],
Print[e]; Print["is first rewritten as"];
Print[norme]; Print["Then, its inverse transform is"] ];
(* take the 1-D inverse Laplace transform *)
newe = myinvlaplace[ norme, s, t, rm, rp, st, fixUp[op] ];
While [ ! SameQ[laste, newe],
If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
laste = newe;
newe = MapAll[transform, laste] ];
(* fix up the result *)
newe = postinvtransform[laste];
While [ ! SameQ[laste, newe],
If [ trace, Print[ newe]; Print[ "which becomes" ] ];
laste = newe;
newe = postinvtransform[laste] ];
If [ trace, Print[ newe] ];
(* Simplification rule based on t being real variable *)
(* Mathematica 2.0 always assumes var to be complex *)
If [ TrueQ[ $VersionNumber >= 2.0 ],
newe = newe /. ( ((t^l_.) x_)^k_ :> t^(l k) x^k /;
( Positive[l] || Negative[l] ) &&
( Positive[k] || Negative[k] ) ) ];
newe /. ( f_ Delta[t + a_.] :> (f /. t -> -a) Delta[t + a] ) ]
(* Kludge around the way Mathematica does Apart *)
(* If the user has selected Apart -> All, then use approximate roots *)
(* Otherwise, inverse transform it as a filter structure *)
(* Enable denominator normalization just in case Apart doesn't do it *)
breakUpRealValuedPoly[x_, s_, t_, rm_, rp_, st_, op_] :=
Block [ {dialogue, input, numer, state,},
state = SetStateField[st, Lmyapartfield, False];
state = SetStateField[state, Lnormalizefield, True];
dialogue = InformUserQ[op];
If [ dialogue,
Print[ "( since the Apart option is not All," ];
Print[ " the inverse Laplace transform of" ];
Print[ " ", x ];
Print[ " will be an IIR filter whose input is" ];
Print[ " the inverse transform of the numerator. )" ] ];
numer = Numerator[x];
input = If [ realValuedPolynomialQ[numer, s] &&
realValuedCoefficientsQ[numer, s],
CFIR[t, CoefficientList[numer, s]],
myinvlaplace[numer, s, t, rm, rp, state, op] ];
CIIR[t, CoefficientList[Denominator[x], s]][input] ] /;
! SameQ[ Replace[Apart, op], All ]
breakUpRealValuedPoly[x_, s_, t_, rm_, rp_, st_, op_] :=
partialFractionsKludge[x, s, t, rm, rp, st, op, N, "approximate"]
(* definitionDialogue *)
definitionDialogue[ oldexpr_, trans_, operator_, options_ ] :=
Block [ {},
If [ InformUserQ[options],
Print[ "( after using ", operator,
" to apply the definition to" ];
Print[ " ", oldexpr ];
Print[ " to find its transform"];
Print[ " ", trans, " )" ] ];
trans ]
(* fixUp *)
fixUp[ op_ ] := { Apart -> Replace[Apart, op],
Definition -> Replace[Definition, op],
Dialogue -> Replace[Dialogue, op],
Simplify -> Replace[Simplify, op],
Terms -> Replace[Terms, op] }
(* multiDInvLaPlace *)
multiDInvLaPlace[ e_, s_, t_, op_ ] :=
MultiDInvTransform [ e, s, t, op, LTransformQ,
MyInvLaPlace, MakeLObject, Global`t ]
(* normalizePoly and normalizeExpression *)
truepoly[x_, s_] := (! FreeQ[x, s]) && PolynomialQ[x, s]
unnormalizedq[x_, s_] := ! SameQ[1, Coefficient[x, s, Exponent[x, s]]]
normalizePoly[poly_, s_, k_] :=
Block [ {coeff, maxpower},
maxpower = Exponent[poly, s];
coeff = Coefficient[poly, s, maxpower];
coeff^k ( Distribute[poly / coeff]^k ) ]
normalizeExpression[x_, s_] :=
Block [ {},
x /. { (poly_Plus)^k_. :> normalizePoly[Expand[poly], s, k] /;
truepoly[poly, s],
poly_ :> Expand[poly] /; truepoly[poly, s] } ]
(* normalizeRatPoly *)
normalizeRatPoly[x_, s_] := normalizeRatPoly[Numerator[x], Denominator[x], s]
normalizeRatPoly[numer_, denom_^k_., s_] :=
Block [ {coeff, maxpower},
maxpower = Exponent[denom, s];
coeff = Coefficient[denom, s, maxpower];
coeff^k numer / ( Distribute[denom / coeff]^k ) ]
(* partialFractionsDialogue *)
partialFractionsDialogue[p_, s_, options_] :=
partialFractionsDialogue[ p, s, options,
normalizeExpression[Apart[p, s], s] ]
partialFractionsDialogue[p_, s_, options_, partfrac_, rootinfo_:"exact"] :=
Block [ {dialogue},
dialogue = InformUserQ[options];
(* Tell the user that we're applying partial fractions *)
If [ dialogue && ! SameQ[p, partfrac],
Print["( after performing partial fraction expansion on"];
Print[" ", p];
Print[" using its ", rootinfo, " roots:" ];
Print[" ", partfrac, " . )" ] ];
partfrac ]
(* partialFractionsKludge -- use MyApart because Apart has failed us *)
partialFractionsKludge[ x_, s_, t_, rm_, rp_, st_, op_,
filter_:Identity, roots_:"exact" ] :=
Block [ {partfrac, result, state},
state = SetStateField[st, Lapartfield, False];
state = SetStateField[state, Lmyapartfield, False];
state = SetStateField[state, Lnormalizefield, False];
partfrac = normalizeExpression[MyApart[x, s, filter], s];
result = partialFractionsDialogue[x, s, op, partfrac, roots];
myinvlaplace[ result, s, t, rm, rp, state, op ] ]
(* postinvtransform *)
postinvtransform[f_] := ( f /. postinvrules )
postinvrules = {
derivative[f_, t_Symbol, k_, False] :>
D[ postinvtransform[f], {t, k} ],
derivative[f_, t_Symbol, k_, v_] :>
Block [ {curderivative, ffun, i, limit, result},
ffun = postinvtransform[f];
result = Unit[k-1][t] Limit[ffun, t -> v];
curderivative = ffun;
For [ i = 2, i <= k, i++,
curderivative = D[curderivative, t];
limit = Limit[curderivative, t -> v];
result += Unit[k - i][t] limit ];
result + D[curderivative, t] ],
integrate[a_ + b_, t_, dummy_Symbol] :>
integrate[postinvtransform[a], t, dummy] +
integrate[postinvtransform[b], t, dummy],
integrate[f_. CStep[dummy_ + b_.], t_, dummy_Symbol ] :>
Integrate[f, {dummy, -b, t}] CStep[t + b],
integrate[f_. CStep[a_. dummy_ + b_.], t_, dummy_Symbol ] :>
( CStep[a t + b] / a ) *
Integrate[f /. t -> t / a,
{dummy, Max[-b, -Infinity a], a t}],
substitutefort[f_, t_, newt_] :>
( postinvtransform[f] /. t -> newt ),
(* These two rules are needed to support a series *)
(* expansion which returns a constant term plus *)
(* terms of the form approx * s^r . *)
myinvlaplace[ a_, s_, t_, rm_, rp_, st_, op_] :> a Delta[t] /;
FreeQ[a, s], (* Converges all ROC *)
myinvlaplace[ a_. s_^r_., s_, t_, rm_, rp_, st_, op_] :>
a Unit[r][t] /;
FreeQ[{a,r}, s], (* Converges all ROC *)
(* Attempt a series expansion about s = 0 *)
(* This is the strategy when all else has failed. *)
(* Please keep this as the last rule. *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
SeriesStrategy[ x, s, t, rm, rp, st, op ] /;
GetStateField[st, Lseriesfield],
x_ :> x
}
(* realAndPosPart *)
realposcheckq[x_] := Implies[ NumberQ[x], RealQ[x] && x > 0 ]
realnegcheckq[x_] := Implies[ NumberQ[x], RealQ[x] && x < 0 ]
realAndPosPart[0] := 0
realAndPosPart[x_Times] := Map[realAndPosPart, x]
realAndPosPart[Complex[a_Integer, b_Integer]] := realAndPosPart[GCD[a,b]]
realAndPosPart[Complex[b_, b_]] := realAndPosPart[b]
realAndPosPart[Complex[0, b_]] := realAndPosPart[b]
realAndPosPart[x_Symbol] := x
realAndPosPart[x_] := x /; realposcheckq[ N[x] ]
realAndPosPart[x_] := -x /; realnegcheckq[ N[x] ]
realAndPosPart[x_] := 1
(* realValuedCoefficientsQ *)
realValuedCoefficientsQ[z_, t_] :=
Apply[Or, Map[RealQ, CoefficientList[z, t]]]
(* realValuedPolynomialQ *)
realValuedPolynomialQ[z_] := realValuedPolynomialQ[z, Global`x]
realValuedPolynomialQ[c_. t_^r_., t_] := IntegerQ[r] && r > 0 && RealValuedQ[c]
realValuedPolynomialQ[c_, t_] := RealValuedQ[c]
realValuedPolynomialQ[x_ + y_, t_] :=
realValuedPolynomialQ[x, t] && realValuedPolynomialQ[y, t]
(* seriesDialogue *)
seriesDialogue[x_, expansion_, terms_, options_] :=
Block [ {dialogue},
dialogue = InformUserQ[options];
If [ dialogue,
Print["( after breaking up the expression"];
Print[" ", x];
Print[" into its series representation"];
Print[" ", expansion];
Print[" the inverse Laplace transform will now be",
" applied to the first ", terms, " terms. )"] ] ]
(* SeriesExpansion *)
SeriesExpansion[ True, x_, s_, power_ ] := Series[x, {s, 0, power}]
SeriesExpansion[ False, x_, s_, power_ ] :=
Block [ {negx, expansion},
negx = x /. s -> s^-1;
expansion = Series[negx, {s, 0, power}];
expansion /. s -> s^-1 ]
(* SeriesStrategy *)
SeriesStrategy[ x_, s_, t_, rm_, rp_, st_, op_ ] :=
Block [ {expansion, exponents, posexpandflag, state, terms},
terms = Replace[Terms, op];
state = SetStateField[st, Lseriesfield, False];
(* If the Terms option is disabled, then take no action. *)
If [ ! terms,
Return[ myinvlaplace[x, s, t, rm, rp, state, op] ] ];
If [ ! IntegerQ[terms],
Message[ InvLaPlace::badterms, terms ];
Return[ myinvlaplace[x, s, t, rm, rp, state, op] ] ];
(* See if we should expand in pos. or neg. powers *)
exponents = GetAllExponents[x, s];
posexpandflag = If [ Apply[And, Thread[exponents >= 0]],
True, False, False ];
(* First expansion *)
expansion = SeriesExpansion[posexpandflag, x, s, terms - 1];
(* Second expansion if first failed *)
If [ ! SameQ[Head[expansion], SeriesData],
expansion = SeriesExpansion[Not[posexpandflag], x, s, terms - 1] ];
(* If this fails, call MyInvLaPlace again *)
(* but with series expansion disabled *)
If [ ! SameQ[Head[expansion], SeriesData],
myinvlaplace[ x, s, t, rm, rp, state, op ],
seriesDialogue[x, expansion, terms, op];
state = nullInvLstate[];
Map [ myinvlaplace[#1, s, t, rm, rp, state, op]&,
Normal[expansion] ] ] ]
(* similarityQ *)
similarityScale[f_, s_] := realAndPosPart[ ScalingFactor[f, s] ]
similarityQ[f_, s_] :=
Block [ {scale},
scale = similarityScale[f, s];
(! SameQ[scale, 1]) && (! SameQ[scale, 0]) ]
(* splitratfun *)
keepterm[ x_ + y_, s_ ] := keepterm[x, s] + keepterm[y, s]
keepterm[ a_. s_^n_., s_ ] := a s^n /; FreeQ[{a, n}, s]
keepterm[ a_, s_ ] := a /; FreeQ[a, s]
keepterm[ a_, s_ ] := 0
splitratfun[x_, s_] := splitratfun[ Numerator[x], Denominator[x], s ]
splitratfun[ numer_, denom_, s_ ] :=
Block [ {expnumer, polyexpr, restexpr},
expnumer = Expand[numer];
polyexpr = keepterm[expnumer, s];
restexpr = expnumer - polyexpr;
{ polyexpr / denom, restexpr / denom } ]
(* transform *)
transform[ expr_ ] :=
If [ SameQ[Head[expr], myinvlaplace],
Replace[expr, InvLaPlaceRules],
expr ]
Format[ myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] ] :=
SequenceForm[ ColumnForm[{"L" Superscript[-1],
" " ~StringJoin~ ToString[s]}],
{ x } ]
Format[ derivative[x_, t_Symbol, k_, flag_] ] :=
SequenceForm[ ColumnForm[{"D" Superscript[k],
" " ~StringJoin~ ToString[t]}],
{ x } ]
Format[ substitutefort[f_, t_, newt_] ] :=
SequenceForm[ {f}, Subscript[t -> newt] ]
Format[ integrate[f_. CStep[a_. dummy_ + b_.], t_, dummy_Symbol ] ] :=
SequenceForm[ "Integrate[", f, ", ", {dummy, 0, t}, "] ",
CStep[a t + b] ]
Format[ integrate[f_, t_, dummy_Symbol ] ] :=
SequenceForm[ "Integrate[", f, ", ", {dummy, 0, t}, "]" ]
(* B E G I N R U L E S *)
InvLaPlaceRules = {}
OriginalInvLaPlaceRules = {
(* I. Rational Transform Pairs *)
(* A. constant functions *)
myinvlaplace[ a_, s_, t_, rm_, rp_, st_, op_] :> a Delta[t] /;
FreeQ[a, s], (* Converges all ROC *)
(* B. constant times exponential *)
myinvlaplace[ a_. Exp[b_. s_], s_, t_, rm_, rp_, st_, op_ ] :>
a Delta[t + b] /;
FreeQ[{a, b}, s],
(* C. inverse of sinusoidal forms *)
myinvlaplace[ c_. ( w_. Cos[w_] + s_ Sin[w_] ) / (s_^2 + w_^2), s_, t_,
rm_, rp_, st_, op_ ] :>
leftOrRightSided[ 0, c Sin[w + w t] CStep[t], t, rp ] /;
FreeQ[{c, w}, s],
myinvlaplace[ c_. ( s_ Cos[b_] - w_. Sin[b_] ) / ( s_^2 + w_^2 ), s_, t_,
rm_, rp_, st_, op_ ] :>
leftOrRightSided[ 0, c Cos[b + w t] CStep[t], t, rp ] /;
FreeQ[{w,b,c}, s],
(* D. sections that are functions of single power of s *)
(* 1. single pole sections *)
myinvlaplace[ b_. (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ a,
b t^-(n+1) Exp[- a t] CStep[t] / Gamma[-n],
t, rp ] /;
FreeQ[{a,b,n}, s] && Negative[n],
myinvlaplace[ b_. / (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ a,
b t^(n-1) Exp[- a t] CStep[t] / Gamma[n],
t, rp ] /;
FreeQ[{a,b,n}, s],
(* 2. single zero sections *)
myinvlaplace[ b_. (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
b Exp[- a t] Unit[n][t] /;
FreeQ[{a,b,n}, s],
(* D. second-order sections *)
(* only encode the cosine/sine pairs because Mma *)
(* will automatically convert Sin[I x] to I Sinh[x] *)
(* and Cos[I x] to Cosh[x]. *)
myinvlaplace[ b_. / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[Sqrt[a], b Sin[Sqrt[a] t] CStep[t] / Sqrt[a], t, rp] /;
FreeQ[{a,b}, s],
myinvlaplace[ b_. s_ / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ Sqrt[a], b Cos[Sqrt[a] t] CStep[t], t, rp ] /;
FreeQ[{a,b}, s],
myinvlaplace[ (s_ + a_)/(s_^2 + c_. s_ + d_), s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {b, invtrans},
b = Sqrt[d - a^2]; (* d = a^2 + b^2 *)
invtrans = Exp[-a t] Cos[b t] CStep[t];
leftOrRightSided[-a, invtrans, t, rp] ] /;
FreeQ[{a,c,d}, s] && SameQ[c, 2 a],
myinvlaplace[ k_. / ( s_^2 + c_. s_ + d_ ), s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {a},
a = c/2; (* complete the square *)
Exp[-a t] myinvlaplace[ k / ( s^2 + d - a^2 ), s, t,
rm, rp, st, op ] ] /;
FreeQ[{c,d,k}, s],
(* E. third order sections *)
myinvlaplace[ b_. /( s_^3 + a_), s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ Max[-Re[(-a)^(1/3)/2], -Re[a^(1/3)]],
b / (3 a^(2/3)) (Exp[- a^(1/3) t] - Exp[a^(1/3) t/2]
( Cos[ a^(1/3) t Sqrt[3] / 2 ] -
Sqrt[3] Sin[ a^(1/3) t Sqrt[3] / 2 ] ) )
CStep[t],
t,
rp ] /;
FreeQ[{a, b}, s], (* [Churchill] *)
(* F. fourth order sections *)
myinvlaplace[ b_. / ( s_^2 + a_ )^2, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {k},
k = Sqrt[a] t;
leftOrRightSided[ Sqrt[-a],
b (Sin[k] - k Cos[k]) CStep[t] / (2 a^(3/2)),
t, rp ] ] /;
FreeQ[{a, b}, s],
myinvlaplace[ b_. s_ / ( s_^2 + a_ )^2, s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ Sqrt[-a],
b t Sin[Sqrt[a] t] CStep[t] / (2 Sqrt[a]),
t, rp ] /;
FreeQ[{a, b}, s],
myinvlaplace[ b_. / ( s_^4 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {k},
k = (a/4)^(1/4);
leftOrRightSided[ Re[ a^(1/4) / Sqrt[2] ],
b (Sin[k t] Cosh[k t] - Cos[k t] Sinh[k t])
CStep[t] / ( 4 k^3 ),
t,
rp ] ] /;
FreeQ[{a, b}, s],
myinvlaplace[ b_. s_ / ( s_^4 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {k},
k = (a/4)^(1/4);
leftOrRightSided[ Re[ a^(1/4) / Sqrt[2] ],
b Sin[k t] Sinh[k t] CStep[t] / ( 2 k^2 ),
t,
rp ] ] /;
FreeQ[{a, b}, s],
(* G. higher order sections *)
myinvlaplace[ b_. s_^n_., s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ 0, b t^-(n + 1) / Gamma[-n] CStep[t], t, rp ] /;
Negative[n] && FreeQ[b, s],
myinvlaplace[ b_. s_^n_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {index = Unique["index"]},
leftOrRightSided[ 0,
b CStep[t] 2^(-n-1/2) t^(-n-1) /
( Product[index, {index, 1, -2 n - 1}] Sqrt[Pi] ),
t, rp ] ] /;
IntegerQ[ - n - 1/2 ] && FreeQ[b, s],
myinvlaplace[ b_. ( s_^2 + a_ )^k_., s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ a,
b Sqrt[Pi] t^(-k - 1/2) / ((2 Sqrt[a])^(-k - 1/2))
BesselJ[-k - 1/2, Sqrt[a] t ] CStep[t] / Gamma[-k],
t, rp ] /;
FreeQ[{a, b}, s] && Negative[k],
myinvlaplace[ b_. / ( s_^2 + a_ )^k_., s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ a,
b Sqrt[Pi] t^(k -1/2) / ((2 Sqrt[a])^(k - 1/2))
BesselJ[k - 1/2, Sqrt[a] t ] CStep[t] / Gamma[k],
t, rp ] /;
FreeQ[{a,b,k}, s],
(* H. partial fractions of rational polynomials *)
(* 1. split expression into rational polys and rest *)
(* this must occur before partial fractions *)
(* because Apart ONLY handles ratio of polys *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {splitlist, state},
splitlist = splitratfun[x, s];
state = SetStateField[st, Lsplitfield, False];
myinvlaplace[ First[splitlist], s, t, rm, rp, state, op] +
myinvlaplace[ Second[splitlist], s, t, rm, rp, state, op] ] /;
GetStateField[st, Lsplitfield] && GetStateField[st, Lapartfield] &&
truepoly[Denominator[x], s] && (! PolynomialQ[Numerator[x], s]),
(* 2. Reduce all rational polynomials into *)
(* a sum of single pole, second order sections, *)
(* and sections like 1 / ( s^4 + a^4 ) *)
(* Enable denominator normalization just in case *)
(* partial fractions does not do it automatically *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {state},
state = SetStateField[st, Lapartfield, False];
state = SetStateField[state, Lnormalizefield, True];
myinvlaplace[ partialFractionsDialogue[x, s, op], s, t,
rm, rp, state, op ] ] /;
GetStateField[st, Lapartfield] && RationalPolynomialQ[x, s],
(* I. partial fractions that Apart could not resolve *)
(* perform MyApart to redo partial fractions for *)
(* denominators with real-valued coefficients *)
(* that are not handled by Apart. *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
breakUpRealValuedPoly[x, s, t, rm, rp, st, op] /;
GetStateField[st, Lmyapartfield] &&
realValuedPolynomialQ[Denominator[x], s] &&
realValuedCoefficientsQ[Denominator[x], s],
(* II. Non-rational Transform Pairs *)
(* A. square root forms *)
myinvlaplace[ 1 / Sqrt[s_ + b_.], s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ 0, Exp[- b t] CStep[t] / (Sqrt[t] Sqrt[Pi]),
t, Infinity ] /;
FreeQ[b, s],
myinvlaplace[ s_ / ( s_ + a_ )^(3/2), s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ -Re[a], Exp[- a t] (1 - 2 a t) CStep[t] / Sqrt[Pi],
t, rp ] /;
FreeQ[a, s],
myinvlaplace[ 1 / ((s_ + a_.)^(3/2) Sqrt[s_ + b_.]), s_, t_,
rm_, rp_, st_, op_ ] :>
leftOrRightSided[ a,
t Exp[- (a + b) t / 2] ( BesselI[0, (a - b) t / 2] -
BesselI[1, (a - b) t / 2] ) CStep[t],
t, rp ] /;
FreeQ[{a,b}, s] && ! SameQ[a, b],
myinvlaplace[ 1 / ( Sqrt[s_] + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
(1 / (Sqrt[Pi] Sqrt[t]) - a Exp[a^2 t](1 - Erf[a Sqrt[t]])) CStep[t] /;
FreeQ[a, s],
(* general form of #38 and #39 [Churchill, 326], *)
(* Sqrt[s] / ( s + a ), augmented by s -> s + b *)
myinvlaplace[ Sqrt[s_ + b_.] / ( s_ + c_ ), s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[c,
( 1 / (Sqrt[Pi] Sqrt[t]) - Sqrt[b - c] Exp[(b - c) t]
Erf[ Sqrt[(b - c) t] ] ) CStep[t],
t, rp ] /;
FreeQ[{b,c}, s],
(* general form of #40 and #41 [Churchill, 326], *)
(* 1 / ( Sqrt[s] (s + a) ), augmented by s -> s + b *)
myinvlaplace[ 1 / (Sqrt[s_ + b_.] (s_ + c_)), s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[c,
Exp[- c t] Erf[ Sqrt[(b - c) t] ] CStep[t] /
Sqrt[b - c],
t, rp ] /;
FreeQ[{b,c}, s],
(* forms that inverse transform to Bessel functions *)
(* 1/Sqrt[s^2 + a] covered as a higher order section *)
myinvlaplace[ ( Sqrt[s_^2 + a_] - s_ )^k_, s_, t_, rm_, rp_, st_, op_ ] :>
leftOrRightSided[ a, k Sqrt[a]^k BesselJ[k, Sqrt[a] t] CStep[t] / t,
t, rp ] /;
FreeQ[{a,k}, s] && Positive[k],
myinvlaplace[ ( Sqrt[s_^2 + a_] - s_ )^v_. / Sqrt[s_^2 + a_], s_, t_,
rm_, rp_, st_, op_ ] :>
Block [ {timefun},
Assuming[ Positive[v + 1], op ];
timefun = Sqrt[a]^v BesselJ[v, Sqrt[a] t] CStep[t];
leftOrRightSided[a, timefun, t, rp] ] /;
FreeQ[{a,v}, s] && Implies[NumberQ[N[v]], N[v > -1]],
myinvlaplace[ 1 / ( Sqrt[s_ + a_.] Sqrt[s_ + b_.] ), s_, t_,
rm_, rp_, st_, op_ ] :>
Exp[- (a + b) t / 2] BesselI[0, (a - b) t / 2] CStep[t] /;
FreeQ[{a, b}, s],
myinvlaplace[ ( s_ - Sqrt[s_^2 + a_] )^v_. / Sqrt[s_^2 + a_], s_, t_,
rm_, rp_, st_, op_ ] :>
Block [ {timefun},
Assuming[ Positive[v + 1], op];
timefun = Sqrt[-a]^v BesselI[v, Sqrt[-a] t] CStep[t];
leftOrRightSided[ a, timefun, t, rp ] ] /;
FreeQ[{a,v}, s] && Implies[NumberQ[N[v]], N[v > -1]],
(* B. hyperbolic forms *)
myinvlaplace[ Tanh[k_. s_] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {n, context},
context = $Context;
$Context = "Global`";
n = Unique["n"]; (* n is an integer *)
$Context = context;
(-1)^(n-1) ( CStep[t - 2 k (n - 1)] - CStep[t - 2 k n] ) ] /;
FreeQ[k, s],
myinvlaplace[ Coth[b_. s_] / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
Abs[Sin[Sqrt[a] t] ] / Sqrt[a] CStep[t] /;
b == ( Pi / ( 2 Sqrt[a] ) ),
(* C. exponential forms *)
myinvlaplace[ Exp[a_. Sqrt[s_]] / Sqrt[s_], s_, t_,
rm_, rp_, st_, op_ ] :>
Exp[-a^2 / (4 t)] CStep[t] / (Sqrt[Pi] Sqrt[t]) /;
FreeQ[a, s],
myinvlaplace[ Exp[k_. / s_] s_^u_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {arg, defarg, exp, order},
defarg = 2 Sqrt[k t];
arg = 2 Sqrt[k] Sqrt[t];
order = -u - 1;
exp = order/2;
Which [ Negative[k],
t^exp BesselJ[order, arg] CStep[t] / k^exp,
Positive[k],
t^exp BesselI[order, arg] CStep[t] / k^exp,
True,
t^exp BesselI[order, defarg] CStep[t] / k^exp ] ] /;
Negative[u],
myinvlaplace[ Exp[k_. / s_], s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {arg, scale},
arg = 2 Sqrt[-k] Sqrt[t];
scale = Sqrt[-k] / ( 2 Sqrt[t] );
Delta[t] + scale ( BesselJ[-1, arg] - BesselJ[1, arg] ) ] /;
FreeQ[k, s],
myinvlaplace[ Exp[ k_ Sqrt[s_] ], s_, t_, rm_, rp_, st_, op_ ] :>
- k Exp[- k^2 / (4 t)] CStep[t] / (2 Sqrt[Pi] Sqrt[t^3]) /;
Negative[k],
(* D. logarithmic forms *)
myinvlaplace[ Log[c_. s_ + b_.] (c_. s_ + b_.)^k_., s_, t_,
rm_, rp_, st_, op_ ] :>
Exp[- b t / Abs[c]] (t / Abs[c])^(-k - 1)
( Gamma[-k] PolyGamma[-k] / Abs[Gamma[-k]]^2 -
Log[t / Abs[c]] / Gamma[-k] ) CStep[t] / Abs[c] /;
FreeQ[{b,c}, s] && Negative[k],
myinvlaplace[ Log[c_. s_ + b_.] / (c_. s_ + b_.)^k_., s_, t_,
rm_, rp_, st_, op_ ] :>
Exp[- b t / Abs[c]] (t / Abs[c])^(k - 1)
( Gamma[k] PolyGamma[k] / Abs[Gamma[k]]^2 -
Log[t / Abs[c]] / Gamma[k] ) CStep[t] / Abs[c] /;
FreeQ[{b,c,k}, s],
myinvlaplace[ Log[b_. s_ + a_] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
( Log[a] - ExpIntegralEi[- a t / Abs[b]] ) CStep[t] /;
FreeQ[{a,b}, s],
myinvlaplace[ Log[b_. s_ + a_.], s_, t_, rm_, rp_, st_, op_ ] :>
- Abs[b] Exp[- a t / Abs[b]] CStep[t] / t /;
FreeQ[{a,b}, s],
myinvlaplace[ Log[a_. + b_. s_^2 ] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
2 ( Log[Sqrt[a]] - CosIntegral[Sqrt[a] t / Sqrt[b]] ) CStep[t] /;
FreeQ[{a,b}, s],
myinvlaplace[ Log[b_. s_] / (c_. s_^2 + 1), s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {time = t / Abs[b]},
( Cos[time] SinIntegral[time] -
Sin[time] CosIntegral[time] ) CStep[t] / Abs[b] ] /;
FreeQ[{b,c}, s] && c == b^2,
myinvlaplace[ s_ Log[b_. s_] / (c_. s_^2 + 1), s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {time = t / Abs[b]},
( - Cos[time] CosIntegral[time] -
Sin[time] SinIntegral[time] ) CStep[t] / (b Abs[b]) ] /;
FreeQ[{b,c}, s] && c == b^2,
myinvlaplace[ Log[ a_ / b_ ], s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ Log[a] - Log[b], s, t, rm, rp, st, op ],
myinvlaplace[ Log[ a_ b_ ], s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ Log[a] + Log[b], s, t, rm, rp, st, op ],
(* E. arc tangent *)
myinvlaplace[ ArcTan[ k_ / s_ ], s_, t_, rm_, rp_, st_, op_ ] :>
Sin[ k t ] CStep[t] / t /; FreeQ[k,s],
(* F. BesselK *)
myinvlaplace[ BesselK[ 0, k_. s_ ], s_, t_, rm_, rp_, st_, op_ ] :>
1 / Sqrt[t^2 - k^2] CStep[t - k] /;
FreeQ[k, s],
myinvlaplace[ BesselK[-1, v_. Sqrt[s_]] / Sqrt[s_], s_, t_,
rm_, rp_, st_, op_ ] :>
Exp[-v^2/(4 t)] CStep[t] / v /;
FreeQ[v, s],
(* III. Transform Properties *)
(* A. homogeneity -- pick off constants *)
myinvlaplace[ a_ f_, s_, t_, rm_, rp_, st_, op_ ] :>
a myinvlaplace[ f, s, t, rm, rp, st, op ] /;
FreeQ[a, s],
(* -. perform MyApart to redo partial fractions for *)
(* denominators with repeated roots that were not *)
(* properly handled by Apart. kludge *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
partialFractionsKludge[x, s, t, rm, rp, st, op] /;
GetStateField[st, Lmyapartfield] && RationalPolynomialQ[x, s],
(* D. additivity *)
myinvlaplace[ x_ + y_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ x, s, t, rm, rp, st, op ] +
myinvlaplace[ y, s, t, rm, rp, st, op ],
myinvlaplace[ (x_ + y_) / c_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ x / c, s, t, rm, rp, st, op ] +
myinvlaplace[ y / c, s, t, rm, rp, st, op ],
myinvlaplace[ a_. Summation[i_, ib_, ie_, inc_][x_], s_, t_,
rm_, rp_, st_, op_ ] :>
Summation[i, ib, ie, inc][ myinvlaplace[ a x, s, t, rm, rp, st, op ] ],
(* E. multiplication by an exponential *)
myinvlaplace[ a_. Exp[b_. s_ + c_.], s_, t_, rm_, rp_, st_, op_ ] :>
substitutefort[ myinvlaplace[ a Exp[c], s, t, rm, rp, st, op ],
t, t + b ] /;
FreeQ[b, s],
(* F. similarity *)
(* only use real, positive scaling factors *)
myinvlaplace[ f_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {newf, newrm, newrp, scale},
scale = similarityScale[f, s];
Assuming[ Positive[scale], op];
newf = f /. s -> s / scale;
newrm = If [ SameQ[rm, -Infinity], rm, scale rm ];
newrp = If [ SameQ[rp, Infinity], rp, scale rp ];
substitutefort[ myinvlaplace[newf, s, t, newrm, newrp, st, op],
t, t / scale ] / scale ] /;
similarityQ[f, s],
(* G. shift *)
myinvlaplace[ f_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {newrm, newrp, normexpr, shift},
{shift, normexpr} = GetShiftFactor[f, s];
normexpr = normexpr /. s -> (s - shift);
newrm = rm + Re[shift];
newrp = rp + Re[shift];
Exp[- shift t] *
myinvlaplace[normexpr, s, t, newrm, newrp, st, op] ] /;
! SameQ[ First[GetShiftFactor[f, s]], 0 ],
(* H. pole in denominator *)
myinvlaplace[ f_ / ( s_ + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ partialFractionsDialogue[f / (s + a), s, op], s, t,
rm, rp, SetStateField[st, Lapartfield, False], op ] /;
GetStateField[st, Lapartfield] && FreeQ[a, s],
myinvlaplace[ f_ / ( s_ + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
Exp[- a t] myinvlaplace[(f /. s -> s - a) / s, s, t, rm, rp, st, op] /;
FreeQ[a, s],
(* I. multiplication by s is progressive *)
(* differentiation in the time domain; *)
(* progressiveness handled by postinvrules *)
myinvlaplace[ s_^k_. a_., s_, t_, -Infinity, Infinity, st_, op_ ] :>
Block [ {dialogue},
dialogue = SameQ[ Replace[Dialogue, op], All ];
Assuming[ k > 0, dialogue ];
If [ dialogue && ! IntegerQ[k],
Print[ "where ", k, " is an integer" ] ];
derivative[ myinvlaplace[a, s, t, rm, rp, st, op],
t, k, False ] ] /;
FreeQ[k, s] && Implies[ NumberQ[k], IntegerQ[k] && ( k > 0 )],
myinvlaplace[ s_^k_. a_., s_, t_, 0, Infinity, st_, op_ ] :>
Block [ {dialogue},
dialogue = SameQ[ Replace[Dialogue, op], All ];
Assuming[ k > 0, dialogue ];
If [ dialogue && ! IntegerQ[k],
Print[ "where ", k, " is an integer" ] ];
derivative[myinvlaplace[a, s, t, rm, rp, st, op], t, k, 0] ] /;
FreeQ[k, s] && Implies[ NumberQ[k], IntegerQ[k] && ( k > 0 )],
(* J. division by s is integration in the time domain *)
myinvlaplace[ f_ s_^n_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {tau = Unique["tau"]},
integrate[ myinvlaplace[f s^(n + 1), s, tau, rm, rp, st, op],
t, tau ] ] /;
IntegerQ[n] && (n < 0),
(* V. Strategies *)
(* A. partial fractions *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ partialFractionsDialogue[x, s, op], s, t, rm, rp,
SetStateField[st, Lapartfield, False], op ] /;
GetStateField[st, Lapartfield],
(* B. Normalize denominator polynomials *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ normalizeRatPoly[x, s], s, t, rm, rp,
SetStateField[st, Lnormalizefield, False], op ] /;
GetStateField[st, Lnormalizefield] && unnormalizedq[Denominator[x], s],
(* C. factor denominator *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {denom, numer},
denom = Denominator[x];
numer = Numerator[x];
denom = Factor[denom];
myinvlaplace[ numer / denom, s, t, rm, rp,
SetStateField[st, Lfactorfield, False], op ] ] /;
GetStateField[st, Lfactorfield],
(* D. distribute expression *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ Distribute[x], s, t, rm, rp, st, op ] /;
SameQ[Head[x], Times] && ! SameQ[x, Distribute[x]],
(* E. expand numerators of expression *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ Expand[x], s, t, rm, rp,
SetStateField[st, Lexpandfield, False], op ] /;
GetStateField[st, Lexpandfield],
(* F. expand all terms in expression *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ ExpandAll[x], s, t, rm, rp,
SetStateField[st, Lexpandallfield, False], op ] /;
GetStateField[st, Lexpandallfield],
(* G. replace new operators and functions *)
(* with their mathematical definition *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
myinvlaplace[ TheFunction[x], s, t, rm, rp,
SetStateField[st, Lthefunfield, False], op ] /;
GetStateField[st, Lthefunfield],
(* H. apply definition if the user has enabled that *)
(* the Definition option. *)
myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
Block [ {state, trans},
state = SetStateField[st, Ldefinitionfield, False];
trans = Integrate[x Exp[s t], {s, rm, rp} ];
If [ SameQ[Head[trans], Integrate],
myinvlaplace[x, s, t, rm, rp, state, op],
definitionDialogue[x, trans/(2 Pi), Integrate, op] ] ] /;
GetStateField[st, Ldefinitionfield] && Replace[Definition, op]
}
(* E N D R U L E S *)
(* E N D P A C K A G E *)
End[]
EndPackage[]
If [ TrueQ[ $VersionNumber >= 2.0 ],
On[ General::spell ];
On[ General::spell1 ] ];
(* A L I A S E S *)
Unprotect[ { InverseLaPlace, InverseLaPlaceTransform } ]
InverseLaPlace = InvLaPlace
InverseLaPlaceTransform = InvLaPlace
Protect[ { InverseLaPlace, InverseLaPlaceTransform } ]
(* H E L P I N F O R M A T I O N *)
Combine[ SPfunctions, { InvLaPlace } ]
Protect[ InvLaPlace ]
(* E N D I N G M E S S A G E *)
Print[ "The inverse Laplace transform rule base InvLaPlace has been loaded." ]
Null
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.