tests/wester_problems/test_transforms.mac

   1 /* Original version of this file copyright 1999 by Michael Wester,
   2  * and retrieved from http://www.math.unm.edu/~wester/demos/Transforms/problems.macsyma
   3  * circa 2006-10-23.
   4  *
   5  * Released under the terms of the GNU General Public License, version 2,
   6  * per message dated 2007-06-03 from Michael Wester to Robert Dodier
   7  * (contained in the file wester-gpl-permission-message.txt).
   8  *
   9  * See: "A Critique of the Mathematical Abilities of CA Systems"
  10  * by Michael Wester, pp 25--60 in
  11  * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
  12  * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
  13  */
  14 /* ----------[ M a c s y m a ]---------- */
  15 /* ---------- Initialization ---------- */
  16 showtime: all$
  17 prederror: false$
  18 /* ---------- Transforms ---------- */
  19 /* Laplace and inverse Laplace transforms
  20    => s/[s^2 + (w - 1)^2]   (Re s > |Im(w - 1)|)
  21       [Gradshteyn and Ryzhik 17.13(33)] */
  22 laplace(cos((w - 1)*t), t, s);
  23 ilt(%, s, t);
  24 /* => w/(s^2 - 4 w^2)   (Re s > |Re w|)   [Gradshteyn and Ryzhik 17.13(84)] */
  25 laplace(sinh(w*t)*cosh(w*t), t, s);
  26 /* e^(-6 sqrt(s))/s   (Re s > 0)   [Gradshteyn and Ryzhik 17.13(102)] */
  27 laplace(erf(3/sqrt(t)), t, s);
  28 specint(exp(-s*t)*erf(3/sqrt(t)), t);
  29 /* Solve y'' + y = 4 [H(t - 1) - H(t - 2)], y(0) = 1, y'(0) = 0 where H is the
  30    Heaviside (unit step) function (the RHS describes a pulse of magnitude 4 and
  31    duration 1).  See David A. Sanchez, Richard C. Allen, Jr. and Walter T.
  32    Kyner, _Differential Equations: An Introduction_, Addison-Wesley Publishing
  33    Company, 1983, p. 211.  First, take the Laplace transform of the ODE
  34    => s^2 Y(s) - s + Y(s) = 4/s [e^(-s) - e^(-2 s)]
  35    where Y(s) is the Laplace transform of y(t) */
  36 atvalue(y(t), t = 0, 1)$
  37 atvalue(diff(y(t), t), t = 0, 0)$
  38 laplace(diff(y(t), t, 2) + y(t) =
  39         4*(unit_step(t - 1) - unit_step(t - 2)), t, s);
  40 expand(%);
  41 /* Now, solve for Y(s) and then take the inverse Laplace transform
  42    => Y(s) = s/(s^2 + 1) + 4 [1/s - s/(s^2 + 1)] [e^(-s) - e^(-2 s)]
  43    => y(t) = cos t + 4 {[1 - cos(t - 1)] H(t - 1) - [1 - cos(t - 2)] H(t - 2)}
  44    */
  45 solve(%, 'laplace(y(t), t, s));
  46 ilt(%[1], s, t);
  47 /* What is the Laplace transform of an infinite square wave?
  48    => 1/s + 2 sum( (-1)^n e^(- s n a)/s, n = 1..infinity )
  49       [Sanchez, Allen and Kyner, p. 213] */
  50 laplace(1 + 2*'sum((-1)^n*unit_step(t - n*a), n, 1, inf), t, s);
  51 /* Fourier transforms => sqrt(2 pi) delta(z)   [Gradshteyn and Ryzhik 17.23(1)]
  52    */
  53 errcatch(fourier_int_coeffs(1, x));
  54 intanalysis: false$
  55 errcatch(fourier_int_coeffs(1, x));
  56 intanalysis: true$
  57 /* => e^(-z^2/36) / [3 sqrt(2)]   [Gradshteyn and Ryzhik 17.23(13)] */
  58 fourier_int_coeffs(exp(-9*x^2), x);
  59 /* => sqrt(2 / pi) (9 - z^2)/(9 + z^2)^2   [Gradshteyn and Ryzhik 17.23(11)] */
  60 fourier_int_coeffs(abs(x)*exp(-3*abs(x)), x);
  61 /* Mellin transforms
  62    => pi cot(pi s)   (0 < Re s < 1)   [Gradshteyn and Ryzhik 17.43(5)] */
  63 MellinTransform(f, x, s):= integrate(f * x^(s - 1), x, 0, inf)$
  64 assume(0 < s, s < 1)$
  65 MellinTransform(1/(1 - x), x, s);
  66 /* => 2^(s - 4) gamma(s/2)/gamma(4 - s/2)   (0 < Re s < 1)
  67       [Gradshteyn and Ryzhik 17.43(16)] */
  68 MellinTransform(bessel_j[3](x)/x^3, x, s);
  69 forget(0 < s, s < 1)$
  70 /* Z transforms.  See _CRC Standard Mathematical Tables_, Twenty-first Edition,
  71    The Chemical Rubber Company, 1973, p. 518.
  72    Z[H(t - m T)] => z/[z^m (z - 1)]   (H is the Heaviside (unit step) function)
  73    */
  74 unit_step(t - 3);
  75 unit_step(t - m);