tests/wester_problems/test_numerical_analysis.mac

   1 /* Original version of this file copyright 1999 by Michael Wester,
   2  * and retrieved from http://www.math.unm.edu/~wester/demos/NumericalAnalysis/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 /* ---------- Numerical Analysis ---------- */
  19 /* This number should immediately simplify to 0.0 */
  20 0.0/sqrt(2);
  21 /* This number normally produces an underflow => 3.29683e-434295 */
  22 exp(-1000000.0);
  23 exp(-1000000.0b0);
  24 /* Arbitrary precision floating point numbers */
  25 oldprec: fpprec$
  26 fpprec: 50$
  27 /* This number is nearly an integer:
  28    26253741 2640768743.9999999999 9925007259 7198185688 ... */
  29 bfloat(exp(sqrt(163)*%pi));
  30 fpprec: oldprec$
  31 remvalue(oldprec)$
  32 /* => [-2, -1] */
  33 [floor(-5/3), ceil(-5/3)];
  34 /* Generate a cubic natural spline s from x = [1, 2, 4, 5] and y = [1, 4, 2, 3]
  35    and then compute s(3) => 27/8 */
  36 spline_coeff([1, 2, 4, 5], [1, 4, 2, 3], s, 4)$
  37 errcatch(spline_interpolate(3, [1, 2, 4, 5], [1, 4, 2, 3], s, 4));
  38 spline_interpolate(3., [1, 2, 4, 5], [1, 4, 2, 3], s, 4);
  39 /* Translation */
  40 p: sum(a[i]*x^i, i, 1, n);
  41 /* Convert into FORTRAN syntax */
  42 fortran(p)$
  43 gentran(p: eval(p))$
  44 /* Convert into C syntax */
  45 gentranlang: 'c$
  46 gentran(p: eval(p))$
  47 gentranlang: 'fortran$
  48 /* Horner's rule---this is important for numerical algorithms
  49    => (a[1] + (a[2] + (a[3] + (a[4] + a[5] x) x) x) x) x */
  50 p: sum(a[i]*x^i, i, 1, 5);
  51 p: horner(p, x);
  52 /* Convert the result into FORTRAN syntax
  53    => p = (a(1) + (a(2) + (a(3) + (a(4) + a(5)*x)*x)*x)*x)*x */
  54 fortran(p)$
  55 gentran(p: eval(p))$
  56 /* Convert the result into C syntax
  57    => p = (a[1] + (a[2] + (a[3] + (a[4] + a[5]*x)*x)*x)*x)*x ; */
  58 gentranlang: 'c$
  59 gentran(p: eval(p))$
  60 gentranlang: 'fortran$
  61 remvalue(p)$
  62 /* Count the number of (floating point) operations needed to compute an
  63    expression => {[+, n - 1], [*, (n^2 - n)/2], [f, (n^2 + n)/2]} */
  64 sum(product(f(i, k), i, 1, k), k, 1, n);
  65 count_ops(%);
  66 /* Interval analysis (interval polynomial example):
  67    ([-4, 2] x + [1, 3])^2 => [-8, 16] x^2 + [-24, 12] x + [1, 9] */
  68 /* Discretize a PDE: for example, forward differencing time (explicit Euler)
  69    and central differencing x on the heat equation =>
  70    (f[i, j+1] - f[i, j])/dt = (f[i+1, j] - 2 f[i, j] + f[i-1, j])/dx^2 */
  71 diff(f(x, t), t) = diff(f(x, t), x, 2);
  72 difference_pde(%, f(x, t), x, f, dx, 'central, t, dt, 'exp_euler);