This is LSolve.m in view mode; [Download] [Up]
(* :Title: Solution of Differential Equations via Laplace ransforms. *)
(* :Authors: Wally McClure, Brian Evans, James McClellan *)
(*
:Summary: This package will solve linear constant coefficient
differential equations by means of the Laplace transform.
*)
(* :Context: SignalProcessing`Analog`LSolve` *)
(* :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 -- September 28, 1990 (adapted from ZSolve.m)
Handle arbitrary initial conditions -- May, 1991
*)
(* :Keywords: linear constant-coefficient differential equations *)
(* :Source: {Operational Mathematics} by Ruel Churchill *)
(* :Warning: *)
(* :Mathematica Version: 1.2 or 2.0 *)
(* :Limitation: *)
(* :Discussion: *)
(* :Functions: PartialL LSolve *)
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`LSolve`",
"SignalProcessing`Analog`LaPlace`",
"SignalProcessing`Analog`InvLaPlace`",
"SignalProcessing`Analog`LSupport`",
"SignalProcessing`Support`TransSupport`",
"SignalProcessing`Support`ROC`",
"SignalProcessing`Support`SigProc`",
"SignalProcessing`Support`SupCode`" ]
(* U S A G E I N F O R M A T I O N *)
LSolve::usage =
"LSolve[ diffequ == drivingfun, y[t] ] solves the differential \
equation diffequ = drivingfun, where diffequ is a linear constant \
coefficient differential equation and drivingfun is the driving \
function (a function of t). \
Thus, diffequ has the form a0 y[t] + a1 y'[t] + .... \
One can specify initial values; e.g., \
LSolve[ y''[t] + 3/2 y'[t] + 1/2 y[t] == Exp[a t], \
y[t], y[0] -> 4, y'[0] -> 10 ]. \
A differential equation of N terms needs N-1 initial conditions. \
All unspecified conditions are considered to be zero. \
LSolve can justify its answers."
(* E N D U S A G E I N F O R M A T I O N *)
Begin[ "`Private`" ]
(* Messages *)
LSolve::nolaplace =
"The forcing function `` does not have a valid Laplace ransform."
LSolve::toomany =
"Initial/boundary conditions occur at more than two locations in ``."
PartialL::nolaplace =
"The differential equation `` does not have a partial Laplace transform."
(* completetransform *)
completetransform[ a_. y_ + b_., y_, x_ ] := ( x - b ) / a /; FreeQ[b, y]
(* ic2index -- extract location of initial/bondary condition *)
SetAttributes[ic2index, Listable]
ic2index[ (y_[index_] -> value_), y_ ] := index
ic2index[ (Derivative[n_][y_][index_] -> value_), y_ ] := index
(* introduction *)
introduction[dialogue_, diffequ_, drivingfun_, bvflag_, initlist_, y_[t_]] :=
Block [ {},
If [ dialogue || SameQ[dialogue, All],
Print["Solving for ",y[t]," in the differential equation"];
Print[" ", diffequ, " = ", drivingfun];
If [ bvflag,
Print[" subject to the boundary conditions"],
Print[" subject to the initial conditions"] ];
If [ EmptyQ[initlist],
Print[" being zero."],
Print[" ", initlist ] ] ] ]
(* preprocess *)
preprocess[ left_ == right_, y_[t_] ] :=
Block [ {diffequ, drivingfun},
diffequ = Select[ left - right, validterm[y[t]] ];
drivingfun = -(left - right - diffequ);
diffequ == drivingfun ]
(* printdialogue *)
printdialogue[dialogue_, y_[t_], Ypartial_, X_, s_, bvflag_] :=
Block [ {formula, ll, n, rules, term1, term2, term3, terml, ul, YYstr},
(* formula for Laplace transform of nth derivative of y(t) *)
YYstr = ToString[StringForm["L { `` }", y[t]]];
term1 = s^n YYstr;
If [ bvflag,
rules = { ul -> "ul", ll -> "ll" };
term2 = s^(n-1) ( y[ul] Exp[-s ul] - y[ll] Exp[-s ll] );
term3 = s^(n-2) ( y'[ul] Exp[-s ul] - y'[ll] Exp[-s ll] );
terml = ( Derivative[n-1][y][ul] Exp[-s ul] -
Derivative[n-1][y][ll] Exp[-s ll] ),
rules = { ll -> "t0" };
term2 = - s^(n-1) y[ll] Exp[-s ll];
term3 = - s^(n-2) y'[ll] Exp[-s ll];
terml = - Derivative[n-1][y][ll] Exp[-s ll] ];
rules = rules ~Join~ { s -> "s", n -> "n" };
term1 = term1 /. rules;
term2 = term2 /. rules;
term3 = term3 /. rules;
terml = terml /. rules; (* yeh, this was a lot of work *)
(* dialogue for user *)
If [ dialogue,
Print["After taking the Laplace transform of both sides"];
Print["and solving for the transform of ", y[t], ":" ] ];
If [ SameQ[dialogue, All],
Print[ " " ];
Print[ "The Laplace transform of the left side is:" ];
Print[ " ",
( Ypartial /. YY -> YYstr ) /. s -> "s" ];
Print[ "( In the general case, the unilateral Laplace" ];
Print[ " transform of the nth derivative of ",
y[t], " is:" ];
Print[ " L {", Derivative["n"][y][t], " } = ", term1,
" + ", term2, " + ", term3, " + ... + ", terml ];
If [ bvflag,
Print[" where t=ll is the lower limit and t=ul is"];
Print[" the upper limit of the boundary conditions. )"],
Print[" where t=t0 is the initial condition. )"] ];
Print[ " " ];
Print[ "The Laplace transform of the right side is:" ];
Print[ " ", X /. s -> "s" ];
Print[ "Solving for the unknown transform yields" ] ] ]
(* validic -- valid initial/bondary condition? *)
validic [y_[t_]] [ y_[index_] -> value_ ] := True
validic [y_[t_]] [ Derivative[n_][y_][index_] -> value_ ] := True
(* validterm -- valid term in left-hand side of a differential equation? *)
validterm [y_[t_]] [ a_. y_[t_ + t0_.] ] := True
validterm [y_[t_]] [ a_. Derivative[n_][y_][t_] ] := True
(* PartialL *)
PartialL[ x_ + v_, rest__ ] :=
PartialL[ x, rest ] + PartialL[ v, rest ]
PartialL[ a_ x_, y_[t_], s_, YY_, rest___ ] :=
a PartialL[ x, y[t], s, YY, rest ] /; FreeQ[a, t]
PartialL[ Derivative[1][y_][t_ + a_.], y_[t_], s_, YY_, ll_, Infinity ] :=
Exp[a s] ( s YY - y[ll] Exp[-s ll] ) /;
FreeQ[a, t]
PartialL[ Derivative[1][y_][t_ + a_.], y_[t_], s_, YY_, ll_, ul_ ] :=
Exp[a s] ( s YY + y[ul] Exp[-s ul] - y[ll] Exp[-s ll] ) /;
FreeQ[a, t] && ! SameQ[ul, Infinity]
PartialL[ Derivative[n_][y_][t_ + a_.], y_[t_], s_, YY_, ll_, Infinity ] :=
Block [ {i, otherterms, term, term1, term2},
term1 = s^n YY;
term2 = - s^(n-1) y[ll] Exp[-s ll];
term = - Derivative[i][y][ll] Exp[-s ll];
otherterms = Apply[Plus, Table[s^(n-i-1) term, {i, 1, n-1}]];
Exp[a s] ( term1 + term2 + otherterms ) ] /;
FreeQ[a, t] && n > 1
PartialL[ Derivative[n_][y_][t_ + a_.], y_[t_], s_, YY_, ll_, ul_ ] :=
Block [ {i, otherterms, term, term1, term2},
term1 = s^n YY;
term2 = s^(n-1) ( y[ul] Exp[-s ul] - y[ll] Exp[-s ll] );
term = ( Derivative[i][y][ul] Exp[-s ul] -
Derivative[i][y][ll] Exp[-s ll] );
otherterms = Apply[Plus, Table[s^(n-i-1) term, {i, 1, n-1}]];
Exp[a s] ( term1 + term2 + otherterms ) ] /;
FreeQ[a, t] && n > 1 && ! SameQ[ul, Infinity]
PartialL[ y_[t_ + a_.], y_[t_], s_, YY_, rest___ ] :=
Exp[a s] YY /; FreeQ[a, t]
PartialL[ a_, y_[t_], s_, YY_, rest___ ] :=
Message[ PartialL::nolaplace, a ] /; FreeQ[a, t]
(* LSolve *)
LSolve/: Options[LSolve] := { Dialogue -> True }
LSolve[ y_[t_ + t0_.] == drivingfun_, y_[t_] ] :=
( drivingfun /. { t -> t - t0 } ) /; FreeQ[t0, t]
LSolve[ left_ == right_, y_[t_], init___ ] :=
Block [ {auginitlist, bvflag, dialogue, diffequ, drivingfun,
equationlist, indexlist, initfun, lowerlimit, options,
rminus, rplus, s, sobj, upperlimit, X, Y, Ycollected,
Ypartial, YY, YYstr},
(* Parse options, determine the level of Dialogue requested *)
options = ToList[init] ~Join~ Options[LSolve];
dialogue = Replace[Dialogue, options];
initlist = Select[ options, validic[y[t]] ];
auginitlist = initlist ~Join~
{ y[t0_] :> 0, Derivative[n_][y][t0_] :> 0 };
(* Collect terms *)
equationlist = preprocess[left == right, y[t]];
diffequ = First[equationlist];
drivingfun = Second[equationlist];
(* Find the location in t of initial/boundary conditions *)
If [ EmptyQ[initlist],
bvflag = False;
lowerlimit = 0;
upperlimit = Infinity,
indexlist = Union[ ic2index[initlist, y] ];
If [ Length[indexlist] > 2,
Message[ LSolve::toomany, t ];
indexlist = Take[indexlist, 2] ];
indexlist = Sort[indexlist];
lowerlimit = First[indexlist];
bvflag = ( Length[indexlist] == 2 );
upperlimit = If [ bvflag, Second[indexlist], Infinity ] ];
(* Tell the user about the differential equation *)
introduction[ dialogue, diffequ, drivingfun,
bvflag, initlist, y[t] ];
(* Find complete Laplace transform of driving function *)
sobj = LaPlace[ drivingfun, t, s, Dialogue -> False ];
If [ ! SameQ[Head[sobj], LTransData],
Message[ LSolve::nolaplace, drivingfun ];
Return[] ];
rminus = GetRMinus[sobj]; (* R- of X(s) *)
rplus = GetRPlus[sobj]; (* R+ of X(s) *)
X = TheFunction[sobj]; (* X(s) from Laplace object *)
(* Find partial Laplace transform of differential equation *)
Ypartial = PartialL[ diffequ, y[t], s, YY,
lowerlimit, upperlimit ];
Ycollected = Collect[Ypartial, YY];
(* Dialogue taking Laplace transform of both sides to user *)
printdialogue[dialogue, y[t], Ycollected, X, s, bvflag];
(* Inverse transform complete Laplace transform of y, Y(s) *)
Y = completetransform[Ycollected, YY, X];
If [ SameQ[dialogue, All],
Print[ " ", Y /. s -> "s" ];
Print[ "which becomes" ] ];
(* Replace initial conditions in Y(s) with values *)
Y = Y /. auginitlist;
If [ dialogue || SameQ[dialogue, All],
Print[ " ", Y /. s -> "s" ];
Print[ "Inverse transforming this gives ", y[t], ":" ] ];
InvLaPlace[{Y, rminus, rplus}, s, t, Dialogue -> False] ]
(* E N D P A C K A G E *)
End[]
EndPackage[]
If [ TrueQ[ $VersionNumber >= 2.0 ],
On[ General::spell ];
On[ General::spell1 ] ];
(* H E L P I N F O R M A T I O N *)
Combine[ SPfunctions, { LSolve } ]
Protect[ LSolve ]
(* E N D I N G M E S S A G E *)
Print[ "LSolve, a differential equation solver, is loaded." ]
Null
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.