/* Original version of this file copyright 1999 by Michael Wester,
 * and retrieved from http://www.math.unm.edu/~wester/demos/Series/problems.macsyma
 * circa 2006-10-23.
 *
 * Released under the terms of the GNU General Public License, version 2,
 * per message dated 2007-06-03 from Michael Wester to Robert Dodier
 * (contained in the file wester-gpl-permission-message.txt).
 *
 * See: "A Critique of the Mathematical Abilities of CA Systems"
 * by Michael Wester, pp 25--60 in
 * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
 * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
 */
/* ----------[ M a c s y m a ]---------- */
/* ---------- Initialization ---------- */
showtime: all$
prederror: false$
/* ---------- Series ---------- */
/* Taylor series---this first example comes from special relativity
   => 1 + 1/2 (v/c)^2 + 3/8 (v/c)^4 + 5/16 (v/c)^6 + O((v/c)^8) */
1/sqrt(1 - (v/c)^2);
taylor(%, v, 0, 7);
1/%^2;
/* Note: sin(x) = x - x^3/6 + x^5/120 - x^7/5040 + O(x^9)
         cos(x) = 1 - x^2/2 + x^4/24 - x^6/720 + O(x^8)
         tan(x) = x + x^3/3 + 2/15 x^5 + 17/315 x^7 + O(x^9) */
tsin: taylor(sin(x), x, 0, 7);
tcos: taylor(cos(x), x, 0, 7);
/* Note that additional terms will be computed as needed */
tsin/tcos;
taylor(tan(x), x, 0, 7);
remvalue(tsin, tcos)$
/* => -x^2/6 - x^4/180 - x^6/2835 - O(x^8) */
taylor(log(sin(x)/x), x, 0, 7);
taylor(sin(x)/x, x, 0, 7);
taylor(log(%), x, 0, 7);
/* => [a f'(a d) + g(b d) + integrate(h(c y), y = 0..d)]
      + [a^2 f''(a d) + b g'(b d) + h(c d)] (x - d) */
diff(f(a*x), x) + g(b*x) + integrate(h(c*y), y, 0, x);
taylor(%, x, d, 1);
ev(%, diff);
/* Taylor series of nonscalar objects (noncommutative multiplication)
   => (B A - A B) t^2/2 + O(t^3)   [Stanly Steinberg] */
declare([A, B], nonscalar)$
%e^((A + B)*t) - %e^(A*t) * %e^(B*t);
ratsimp(%);
taylor(%e^((A + B)*t) - %e^(A*t) * %e^(B*t), t, 0, 4);
remove([A, B], nonscalar)$
/* Laurent series:
   => sum( Bernoulli[k]/k! x^(k - 2), k = 1..infinity )
      = 1/x^2 - 1/(2 x) + 1/12 - x^2/720 + x^4/30240 + O(x^6)
      [Levinson and Redheffer, p. 173] */
taylor(1/(x*(exp(x) - 1)), x, 0, 7);
/* Puiseux series (terms with fractional degree):
   => 1/sqrt(x - 3/2 pi) + (x - 3/2 pi)^(3/2) / 12 + O([x - 3/2 pi]^(7/2)) */
taylor(sqrt(sec(x)), x, 3/2*%pi, 4);
/* Generalized Taylor series => sum( [x log x]^k/k!, k = 0..infinity ) */
taylor(x^x, x, 0, 3);
/* Compare the generalized Taylor series of two different formulations of a
   function => log(z) + log(cosh(w)) + tanh(w) z + O(z^2) */
s1: taylor(log(sinh(z)) + log(cosh(z + w)), z, 0, 1);
s2: taylor(log(sinh(z) * cosh(z + w)), z, 0, 1);
radcan(s1 - s2);
remvalue(s1, s2)$
/* Look at the generalized Taylor series around x = 1
   => (x - 1)^a/e^b [1 - (a + 2 b) (x - 1) / 2 + O((x - 1)^2)] */
log(x)^a*exp(-b*x);
taylor(%, x, 1, 1);
/* Asymptotic expansions => sqrt(2) x + O(1/x) */
taylor(sqrt(2*x^2 + 1), x, inf, 0);
/* Wallis' product => 1/sqrt(pi n) + ...   [Knopp, p. 385] */
errcatch(taylor(1/2^(2*n) * binomial(2*n, n), n, inf, 0));
/* => 0!/x - 1!/x^2 + 2!/x^3 - 3!/x^4 + O(1/x^5)   [Knopp, p. 544] */
exp(x) * 'integrate(exp(-t)/t, t, x, inf);
errcatch(taylor(%, x, inf, 5));
/* Multivariate Taylor series expansion => 1 - (x^2 + 2 x y + y^2)/2 + O(x^4)
   */
taylor(cos(x + y), [x, y], 0, 3);
taylor(cos(x + y), [x, 0, 3], [y, 0, 3]);
expand(%);
/* Power series (compute the general formula) */
verbose: true$
powerseries(log(sin(x)/x), x, 0);
taylor(%, x, 0, 7);
powerseries(exp(-x)*sin(x), x, 0);
taylor(%, x, 0, 7);
/* Derive an explicit Taylor series solution of y as a function of x from the
   following implicit relation:
   y = x - 1 + (x - 1)^2/2 + 2/3 (x - 1)^3 + (x - 1)^4 + 17/10 (x - 1)^5 + ...
   */
x = sin(y) + cos(y);
taylor_revert(%, y, 7);
/* Pade (rational function) approximation => (2 - x)/(2 + x) */
pade(taylor(exp(-x), x, 0, 2), 1, 1);
/* Fourier series of f(x) of period 2 p over the interval [-p, p]
   => - (2 p / pi) sum( (-1)^n sin(n pi x / p) / n, n = 1..infinity ) */
assume(p > 0)$
fourier_series(x, x, p);
/* => p / 2
   - (2 p / pi^2) sum( [1 - (-1)^n] cos(n pi x / p) / n^2, n = 1..infinity ) */
fourier_series(abs(x), x, p);
forget(p > 0)$
/* ---------- Quit ---------- */
quit();
