(* 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
