This is Laplace.m in view mode; [Download] [Up]
(* Laplace.m *)
(* Copyright 1990 Halchin and Fleming *)
(* by Charles G. Fleming *)
(* Created November 21, 1990 *)
(* last modified February 15, 1991 *)
(* This program solves the second order linear differential equation that *)
(* models a mass attached to a spring and subject to an external forcing *)
(* function f Sin[w t]. The solutions are obtained using the Laplace *)
(* transform method. *)
BeginPackage["Laplace`"]
Laplace::usage= "
Laplace[m, b, k, f, w, x[0], x'[0] ] will generate the solution to the
differential equation that models a mass attached to a spring and subject
to the external forcing function f Sin[w t]."
Begin["`Private`"]
Laplace[m_, b_, k_, f_, w_, i1_, i2_] :=
Block[{c1, c2, c3, w0, s, rule, LT, num, discrim, r1, r2, LT1,
AA, BB, CC, DD, num1, poly, d0, d1, d2, d3, rule2, x, e1, e2},
(* Some initializations. *)
c1 = b/m;
c2 = k/m;
c3 = f/m;
w0 = Sqrt[k/m];
x[0] = i1;
x'[0] = i2;
(* Compute the Laplace transform of the differential equation. *)
rule = Solve[s^2 LT - s x[0] - x'[0] + c1 (s LT - x[0]) + c2 LT -
c3 w /(s^2 + w^2) == 0, LT];
LT = LT /. rule[[1]] ;
LT = Simplify[LT];
num = Numerator[LT];
(* Determine the FORM of the partial fraction decomposition of *)
(* the Laplace transform and the FORM of the solution to the deq. *)
discrim = b^2 - 4 m k;
(* This is the case where a forcing function is present. *)
If[ f != 0,
Which[
(* Overdamped *)
discrim > 0,
r1 = (-b + Sqrt[discrim])/(2 m);
r2 = (-b - Sqrt[discrim]) / (2 m);
LT1 = AA/(s-r1) + BB/(s-r2) + (CC s + DD) /(s^2 + w^2);
x[Global`t] = AA E^(r1 Global`t) + BB E^(r2 Global`t) +
CC Cos[w Global`t] + DD/w Sin[w Global`t];,
(* Underdamped *)
discrim < 0,
If[ c1 == 0 && w0 == w,
LT1 = (AA s + BB)/(s^2 + w^2) + (CC s + DD)/(s^2 + w^2)^2;
x[Global`t] = AA Cos[w Global`t] + BB Sin[w Global`t] +
CC/(2 w) Global`t Sin[w Global`t] +
DD/(2 w^3) (Sin[w Global`t] -
w Global`t Cos[w Global`t]);,
e1 = -c1/2;
e2 = Sqrt[c2 - d1^2];
LT1 = (AA s + BB) / (c2 + c1 s + s^2) + (CC s + DD) /
(s^2 + w^2);
x[Global`t] = E^(e1 Global`t) (AA Cos[e2 Global`t] +
(BB + e1 AA)/e2 Sin[e2 Global`t]) +
CC Cos[w Global`t] + DD/w Sin[w Global`t];
],
(* Critically damped *)
True,
r1 = (-b + Sqrt[discrim])/(2 m);
LT1 = AA/(s-r1) + BB/(s-r1)^2 + (CC s + DD) / (s^2 + w^2);
x[Global`t] = AA E^(r1 Global`t) + BB Global`t E^(r1 Global`t) +
CC Cos[w Global`t] + DD/w Sin[w Global`t];
],
(* This is the case where no forcing function is present. *)
Which[
(* Overdamped *)
discrim > 0,
r1 = (-b + Sqrt[discrim])/(2 m);
r2 = (-b - Sqrt[discrim]) / (2 m);
LT1 = AA/(s-r1) + BB/(s-r2);
x[Global`t] = AA E^(r1 Global`t) + BB E^(r2 Global`t);,
(* Underdamped *)
discrim < 0,
e1 = -c1/2;
e2 = Sqrt[c2 - d1^2];
LT1 = (AA s + BB) / (c2 + c1 s + s^2);
x[Global`t] = E^(e1 Global`t) (AA Cos[e2 Global`t] +
(BB + e1 AA)/e2 Sin[e2 Global`t]);,
(* Critically damped *)
True,
r1 = (-b + Sqrt[discrim])/(2 m);
LT1 = AA/(s-r1) + BB/(s-r1)^2;
x[Global`t] = AA E^(r1 Global`t) + BB Global`t E^(r1 Global`t);
];
];
(* Calculate the partial fraction decomposition of the Laplace *)
(* transform.*)
LT1 = Together[LT1];
num1 = Numerator[LT1];
poly = Collect[ num - num1, s];
d0 = Coefficient[poly, s, 0];
d1 = Coefficient[poly, s, 1];
d2 = Coefficient[poly, s, 2];
d3 = Coefficient[poly, s, 3];
rule2 = Solve[{d0==0, d1==0, d2==0, d3==0}, {AA, BB, CC, DD}];
rule2 = Simplify[rule2];
(* Return the solution with the unknown constants determined. *)
AA = AA /. rule2[[1]];
BB = BB /. rule2[[1]];
CC = CC /. rule2[[1]];
DD = DD /. rule2[[1]];
Return[Simplify[x[Global`t]]];
]
End[]
EndPackage[]
Laplace
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.