tests/wester_problems/test_complex_domain.mac

   1 /* Original version of this file copyright 1999 by Michael Wester,
   2  * and retrieved from http://www.math.unm.edu/~wester/demos/ComplexDomain/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 /* ---------- The Complex Domain ---------- */
  19 /* Complex functions---separate into their real and imaginary parts.
  20    Here, variables default to REAL.
  21    [Re(x + i y), Im(x + i y)] => [Re(x) - Im(y), Im(x) + Re(y)]
  22    for x and y complex */
  23 [realpart(x + %i*y), imagpart(x + %i*y)];
  24 declare([x, y], complex)$
  25 [realpart(x + %i*y), imagpart(x + %i*y)];
  26 remove([x, y], complex)$
  27 /* => 1   [W. Kahan] */
  28 abs(3 - sqrt(7) + %i*sqrt(6*sqrt(7) - 15));
  29 ratsimp(%);
  30 /* => 1/sqrt(a^2 + (1/a + b)^2) for real a, b */
  31 abs(1/(a + %i/a + %i*b));
  32 /* => log 5 + i arctan(4/3) */
  33 rectform(log(3 + 4*%i));
  34 /* => [sin(x) cos(x) + i sinh(y) cosh(y)] / [cos(x)^2 + sinh(y)^2] */
  35 rectform(tan(x + %i*y));
  36 /* Check for branch abuse.  See David R. Stoutemyer, ``Crimes and Misdemeanors
  37    in the Computer Algebra Trade'', _Notices of the American Mathematical
  38    Society_, Volume 38, Number 7, September 1991, 778--785.  This first
  39    expression can simplify to sqrt(x y)/sqrt(x), but no further in general
  40    (consider what happens when x, y = -1).  sqrt(x y) = sqrt(x) sqrt(y) if
  41    either x >= 0 or y >= 0 or both x and y lie in the right-half plane
  42    (Re x, Re y > 0) [considering principal values]. */
  43 expr: sqrt(x*y*abs(z)^2) / (sqrt(x)*abs(z));
  44 ratsimp(%);
  45 declare([x, y, z], complex)$
  46 ratsimp(expr);
  47 remove(y, complex)$
  48 /* Special case: sqrt(x y |z|^2)/(sqrt(x) |z|) => sqrt(y) [PV] for y >= 0 */
  49 assume(y >= 0)$
  50 sqrt(x*y*abs(z)^2) / (sqrt(x)*abs(z));
  51 forget(y >= 0)$
  52 remove([x, z], complex)$
  53 /* sqrt(1/z) = 1/sqrt(z) except when z is real and negative, in which case
  54    sqrt(1/z) = - 1/sqrt(z) [considering principal values] */
  55 expr: sqrt(1/z) - 1/sqrt(z);
  56 ratsimp(%);
  57 declare(z, complex)$
  58 ratsimp(expr);
  59 remove(z, complex)$
  60 /* Special case: sqrt(1/z) - 1/sqrt(z) => 0 [PV] for z > 0 */
  61 assume(z > 0)$
  62 ratsimp(expr);
  63 forget(z > 0)$
  64 /* Special case: sqrt(1/z) + 1/sqrt(z) => 0 [PV] for z < 0 */
  65 assume(z < 0)$
  66 sqrt(1/z) + 1/sqrt(z);
  67 forget(z < 0)$
  68 /* sqrt(e^z) = e^(z/2) if and only if Im z is contained in the interval
  69    ((4 n - 1) pi, (4 n + 1) pi] for n an integer: ..., (-5 pi, -3 pi],
  70    (-pi, pi], (3 pi, 5 pi], ...; otherwise, sqrt(e^z) = - e^(z/2) [considering
  71    principal values] */
  72 declare(z, complex)$
  73 sqrt(%e^z) - %e^(z/2);
  74 ratsimp(%);
  75 remove(z, complex)$
  76 /* Special case: sqrt(e^z) - e^(z/2) => 0 [PV] for z real */
  77 sqrt(%e^z) - %e^(z/2);
  78 /* The principal value of this expression is - e^(3 i) = - cos 3 - i sin 3 */
  79 sqrt(%e^(6*%i));
  80 bfloat(%);
  81 /* log(e^z) = z if and only if Im z is contained in the interval (-pi, pi]
  82    [considering principal values] */
  83 declare(z, complex)$
  84 log(%e^z);
  85 remove(z, complex)$
  86 /* Special case: log(e^z) => z [PV] for z real */
  87 log(%e^z);
  88 /* The principal value of this expression is (10 - 4 pi) i */
  89 log(%e^(10*%i));
  90 /* (x y)^n = x^n y^n if either x > 0 or y > 0 or both x and y lie in the
  91    right-half plane (Re x, Re y > 0) or n is an integer [considering principal
  92    values] */
  93 expr: (x*y)^(1/n) - x^(1/n)*y^(1/n);
  94 ratsimp(%);
  95 declare([x, y], complex)$
  96 ratsimp(expr);
  97 remove(y, complex)$
  98 /* Special case: (x y)^(1/n) - x^(1/n) y^(1/n) => 0 [PV] for y > 0 */
  99 assume(y > 0)$
 100 ratsimp(expr);
 101 forget(y > 0)$
 102 /* Special case: (x y)^n - x^n y^n => 0 [PV] for integer n */
 103 declare(y, complex, n, integer)$
 104 (x*y)^n - x^n*y^n;
 105 remove([x, y], complex, n, integer)$
 106 /* arctan(tan(z)) = z for z real if and only if z is contained in the interval
 107    (-pi/2, pi/2] [considering principal values] */
 108 expr: atan(tan(z));
 109 ratsimp(%);
 110 declare(z, complex)$
 111 ratsimp(expr);
 112 remove(z, complex)$
 113 /* Special case: arctan(tan(z)) => z [PV] for -pi/2 < z < pi/2 */
 114 assume(-%pi/2 < z, z < %pi/2)$
 115 ratsimp(expr);
 116 forget(-%pi/2 < z, z < %pi/2)$
 117 ev(expr, triginverses: all);
 118 remvalue(expr)$
 119 /* The principal value of this expression is 10 - 3 pi */
 120 atan(tan(10));
 121 /* The principal value of this expression is 11 - 4 pi + 30 i = -1.56637 + 30 i
 122    */
 123 atan(tan(11 + 30*%i));
 124 atan(tan(11.0 + 30.0*%i));
 125 sfloat(%);
 126 /* This is a challenge problem proposed by W. Kahan: simplify the following
 127    expression for complex z.  Expanding out the expression produces
 128    (z^2 + 1)/(2 z) +- (z + 1)*(z - 1)/(2 z) => z or 1/z in each of its branches
 129    */
 130 declare(z, complex)$
 131 w: (z + 1/z)/2;
 132 expr: w + sqrt(w + 1)*sqrt(w - 1);
 133 ratsimp(expr);
 134 radcan(expr);
 135 remvalue(w, expr)$
 136 remove(z, complex)$