Initial population.

1 files changed, 395 insertions, 0 deletions

diff --git a/src/input/wester.input.pamphlet b/src/input/wester.input.pamphlet
new file mode 100644
index 00000000..5e08228b
--- /dev/null
+++ b/src/input/wester.input.pamphlet
@@ -0,0 +1,395 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input wester.input}
+\author{Michael Wester}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{License}
+<<license>>=
+--Copyright The Numerical Algorithms Group Limited 1996.
+@
+<<*>>=
+<<license>>
+
+
+-- ----------[ A x i o m ]----------
+-- ---------- Initialization ----------
+)set messages autoload off
+)set messages time on
+)set quit unprotected
+)set streams calculate 7
+-- ---------- Numbers ----------
+-- Let's begin by playing with numbers: infinite precision integers
+factorial(50)
+factor(%)
+-- Infinite precision rational numbers
+1/2 + 1/3 + 1/4 + 1/5 + 1/6 + 1/7 + 1/8 + 1/9 + 1/10
+-- Arbitrary precision floating point numbers
+digits(50);
+-- This number is nearly an integer
+exp(sqrt(163.)*%pi)
+digits(20);
+-- Special functions
+besselJ(2, 1 + %i)
+-- Complete decimal expansion of a rational number
+decimal(1/7)
+-- Continued fractions
+continuedFraction(%pi)
+-- Simplify an expression with nested square roots
+s1:=sqrt(2*sqrt(3) + 4)
+p:POLY FRAC INT:= (ratPoly(s1::Expression Integer)::SUP FRAC INT).'z
+solp:=radicalSolve p
+rhs select (z+-> real abs (complexNumeric rhs z - complexNumeric s1) < 1.E-19,solp).1
+simplify(s1)
+-- Try a more complicated example (from the Putnam exam)
+s1:=sqrt(14 + 3*sqrt(3 + 2*sqrt(5 - 12*sqrt(3 - 2*sqrt(2)))))
+p:POLY FRAC INT:= (ratPoly(s1::Expression Integer)::SUP FRAC INT).'z
+solp:=radicalSolve p
+rhs select (z+-> real abs (complexNumeric rhs z - complexNumeric s1) < 1.E-19,solp).1
+simplify(s1)
+-- Cardinal numbers
+2*Aleph(0) - 3
+-- ---------- Statistics ----------
+-- ---------- Algebra ----------
+-- Numbers are nice, but symbols allow for variability---try some high school
+-- algebra: rational simplification
+(x**2 - 4)/(x**2 + 4*x + 4)
+-- This example requires more sophistication
+(%e**x - 1)/(%e**(x/2) + 1)
+normalize(%)
+-- Expand and factor polynomials
+(x + 1)**20
+D(%, x)
+factor(%)
+x**100 - 1
+factor(%)
+-- Factor polynomials over finite fields and field extensions
+p:= x**4 - 3*x**2 + 1
+factor(p)
+phi:= rootOf(phi**2 - phi - 1);
+factor(p, [phi])
+factor(p :: Polynomial(PrimeField(5)))
+expand(%)
+-- Partial fraction decomposition
+(x**2 + 2*x + 3)/(x**3 + 4*x**2 + 5*x + 2)
+padicFraction(
+   partialFraction(numerator(%) :: UnivariatePolynomial(x, Fraction Integer),
+                   factor(denominator(%) :: Polynomial Integer) ::
+                      Factored UnivariatePolynomial(x, Fraction Integer)))
+-- ---------- Inequalities ----------
+-- ---------- Trigonometry ----------
+-- Trigonometric manipulations---these are typically difficult for students
+r:= cos(3*x)/cos(x)
+real(complexNormalize(%))
+real(normalize(simplify(complexNormalize(r))))
+-- Use rewrite rules
+sincosAngles:= rule _
+  (cos((n | integer?(n)) * x) == _
+      cos((n - 1)*x) * cos(x) - sin((n - 1)*x) * sin(x); _
+   sin((n | integer?(n)) * x) == _
+      sin((n - 1)*x) * cos(x) + cos((n - 1)*x) * sin(x) )
+sincosAngles r
+-- ---------- Determining Zero Equivalence ----------
+-- The following expressions are all equal to zero
+sqrt(997) - (997**3)**(1/6)
+sqrt(999983) - (999983**3)**(1/6)
+s1:=(2**(1/3) + 4**(1/3))**3 - 6*(2**(1/3) + 4**(1/3)) - 6
+simplify(%)
+p:POLY FRAC INT:= (ratPoly(s1::Expression Integer)::SUP FRAC INT).'z
+solp:=radicalSolve p
+rhs select (z+-> real abs (complexNumeric rhs z - complexNumeric s1) < 1.E-19,solp).1
+-- This expression is zero for x, y > 0 and n not equal to zero
+x**(1/n)*y**(1/n) - (x*y)**(1/n)
+normalize(%)
+-- See Joel Moses, ``Algebraic Simplification: A Guide for the Perplexed'',
+-- CACM, Volume 14, Number 8, August 1971
+expr:= log(tan(1/2*x + %pi/4)) - asinh(tan(x))
+complexNormalize(%)
+-- Use a roundabout method---show that expr is a constant equal to zero
+D(expr, x)
+normalize(rootSimp(expand(simplify(%))))
+normalize(eval(expr, x = 0))
+expr:=log((2*sqrt(r) + 1)/sqrt(4*r + 4*sqrt(r) + 1))
+D(expr, x)
+eval(expr, x = 0)
+(4*r + 4*sqrt(r) + 1)**(sqrt(r)/(2*sqrt(r) + 1)) _
+   * (2*sqrt(r) + 1)**(1/(2*sqrt(r) + 1)) - 2*sqrt(r) - 1
+normalize(%)
+-- ---------- The Complex Domain ----------
+-- Complex functions---separate into their real and imaginary parts
+rectform(z) == real(z) + %i*imag(z)
+rectform(log(3 + 4*%i))
+simplify(rectform(tan(x + %i*y)))
+-- Check for branch abuse.  See David R. Stoutemyer, ``Crimes and Misdemeanors
+-- in the Computer Algebra Trade'', Notices of the AMS, Volume 38, Number 7,
+-- September 1991.  This first expression can simplify to sqrt(x y)/sqrt(x),
+-- but no further in general (consider what happens when x, y = -1).
+sqrt(x*y*abs(z)**2) / (sqrt(x)*abs(z))
+rootSimp %
+-- If z = -1, sqrt(1/z) is not equal to 1/sqrt(z)
+sqrt(1/z) - 1/sqrt(z)
+-- If z = 3 pi i, log(exp(z)) is not equal to z
+log(%e**z)
+normalize(%)
+-- The principal value of this expression is (10 - 4 pi) i
+log(%e**(10*%i))
+normalize(%)
+-- If z = pi, arctan(tan(z)) is not equal to z
+atan(tan(z))
+-- If z = 2 pi i, sqrt(exp(z)) is not equal to exp(z/2)
+sqrt(%e**z) - %e**(z/2)
+-- ---------- Equations ----------
+-- Manipulate an equation using a natural syntax
+(x = 0)/2 + 1
+-- Solve various nonlinear equations---this cubic polynomial has all real roots
+radicalSolve(3*x**3 - 18*x**2 + 33*x - 19 = 0, x)
+map(e +-> lhs(e) = rectform(rhs(e)), %)
+-- Some simple seeming problems can have messy answers
+eqn:= x**4 + x**3 + x**2 + x + 1 = 0
+radicalSolve(eqn, x)
+-- Check one of the answers
+eval(eqn, %.1)
+%e**(2*x) + 2*%e**x + 1 = z
+solve(%, x)
+-- This equation is already factored and so *should* be easy to solve
+(x + 1) * (sin(x)**2 + 1)**2 * cos(3*x)**3 = 0
+solve(%, x)
+-- The following equations have an infinite number of solutions (let n be an
+-- arbitrary integer): z = 0 [+ n 2 pi i]
+solve(%e**z = 1, z)
+-- x = pi/4 [+ n pi]
+solve(sin(x) = cos(x), x)
+solve(tan(x) = 1, x)
+-- x = 0, 0 [+ n pi, + n 2 pi]
+solve(sin(x) = tan(x), x)
+-- This equation has no solutions
+solve(sqrt(x**2 + 1) = x - 2, x)
+-- Solve a system of linear equations
+eq1:=   x +   y +   z =  6
+eq2:= 2*x +   y + 2*z = 10
+eq3:=   x + 3*y +   z = 10
+-- Note that the solution is parametric
+solve([eq1, eq2, eq3], [x, y, z])
+-- Solve a system of nonlinear equations
+eq1:= x**2*y + 3*y*z - 4 = 0
+eq2:= -3*x**2*z + 2*y**2 + 1 = 0
+eq3:= 2*y*z**2 - z**2 - 1 = 0
+-- Solving this by hand would be a nightmare
+solve([eq1, eq2, eq3], [x, y, z])
+-- ---------- Matrix Algebra ----------
+m:= matrix([[a, b], [1, a*b]])
+-- Invert the matrix
+minv:= inverse(m)
+m * minv
+-- Define a Vandermonde matrix (useful for doing polynomial interpolations)
+matrix([[1,    1,    1,    1   ], _
+        [w,    x,    y,    z   ], _
+        [w**2, x**2, y**2, z**2], _
+        [w**3, x**3, y**3, z**3]])
+determinant(%)
+-- The following formula implies a general result
+factor(%)
+-- Compute the eigenvalues of a matrix from its characteristic polynomial
+m:= matrix([[ 5, -3, -7], _
+            [-2,  1,  2], _
+            [ 2, -3, -4]])
+characteristicPolynomial(m, lambda)
+solve(% = 0, lambda)
+-- ---------- Tensors ----------
+-- ---------- Sums and Products ----------
+-- Sums: finite and infinite
+summation(k**3, k = 1..n)
+sum(k**3, k = 1..n)
+limit(sum(1/k**2 + 1/k**3, k = 1..n), n = %plusInfinity)
+-- Products
+product(k, k = 1..n)
+-- ---------- Calculus ----------
+-- Limits---start with a famous example
+limit((1 + 1/n)**n, n = %plusInfinity)
+limit((1 - cos(x))/x**2, x = 0)
+-- Apply the chain rule---this is important for PDEs and many other
+-- applications
+y:= operator('y);
+x:= operator('x);
+D(y(x(t)), t, 2)
+)clear properties x y
+-- ---------- Indefinite Integrals ----------
+1/(x**3 + 2)
+-- This would be very difficult to do by hand
+integrate(%, x)
+D(%, x)
+-- This example involves several symbolic parameters
+integrate(1/(a + b*cos(x)), x)
+map(simplify, map(f +-> D(f, x), %))
+-- Calculus on a non-smooth (but well defined) function
+D(abs(x), x)
+integrate(abs(x), x)
+-- Calculus on a piecewise defined function
+a(x) == if x < 0 then -x else x
+D(a(x), x)
+integrate(a(x), x)
+)clear properties a
+-- The following two integrals should be equivalent.  The correct solution is
+-- [(1 + x)^(3/2) + (1 - x)^(3/2)] / 3
+integrate(x/(sqrt(1 + x) + sqrt(1 - x)), x)
+integrate((sqrt(1 + x) - sqrt(1 - x))/2, x)
+-- ---------- Definite Integrals ----------
+-- The following two functions have a pole at zero
+integrate(1/x, x = -1..1)
+integrate(1/x**2, x = -1..1)
+-- Different branches of the square root need to be chosen in the intervals
+-- [0, 1] and [1, 2].  The correct results are 4/3, [4 - sqrt(8)]/3,
+-- [8 - sqrt(8)]/3, respectively.
+integrate(sqrt(x + 1/x - 2), x = 0..1)
+integrate(sqrt(x + 1/x - 2), x = 0..1, "noPole")
+integrate(sqrt(x + 1/x - 2), x = 1..2)
+integrate(sqrt(x + 1/x - 2), x = 1..2, "noPole")
+integrate(sqrt(x + 1/x - 2), x = 0..2)
+integrate(sqrt(x + 1/x - 2), x = 0..2, "noPole")
+-- Contour integrals
+integrate(cos(x)/(x**2 + a**2), x = %minusInfinity..%plusInfinity)
+integrate(cos(x)/(x**2 + a**2), x = %minusInfinity..%plusInfinity, "noPole")
+-- Integrand with a branch point
+integrate(t**(a - 1)/(1 + t), t = 0..%plusInfinity)
+integrate(t**(a - 1)/(1 + t), t = 0..%plusInfinity, "noPole")
+-- Multiple integrals: volume of a tetrahedron
+integrate(integrate(integrate(1, z = 0..c*(1 - x/a - y/b)), _
+                    y = 0..b*(1 - x/a)), _
+          x = 0..a)
+-- ---------- Series ----------
+-- Taylor series---this first example comes from special relativity
+1/sqrt(1 - (v/c)**2)
+series(%, v = 0)
+1/%**2
+tsin:= series(sin(x), x = 0)
+tcos:= series(cos(x), x = 0)
+-- Note that additional terms will be computed as needed
+tsin/tcos
+series(tan(x), x = 0)
+-- Look at the Taylor series around x = 1
+)set streams calculate 1
+log(x)**a*exp(-b*x)
+series(%, x = 1)
+)set streams calculate 7
+-- Compare the Taylor series of two different formulations of a function
+taylor(log(sinh(z)) + log(cosh(z + w)), z = 0)
+% - taylor(log(sinh(z) * cosh(z + w)), z = 0)
+series(log(sinh(z)) + log(cosh(z + w)), z = 0)
+% - series(log(sinh(z) * cosh(z + w)), z = 0)
+-- Power series (compute the general formula)
+log(sin(x)/x)
+series(%, x = 0)
+exp(-x)*sin(x)
+series(%, x = 0)
+-- Derive an explicit Taylor series solution of y as a function of x from the
+-- following implicit relation
+y:= operator('y);
+x = sin(y(x)) + cos(y(x))
+seriesSolve(%, y, x = 1, 0)
+)clear properties y
+-- Pade (rational function) approximation
+pade(1, 1, taylor(exp(-x), x = 0))
+-- ---------- Transforms ----------
+-- Laplace and inverse Laplace transforms
+laplace(cos((w - 1)*t), t, s)
+inverseLaplace(%, s, t)
+-- ---------- Difference and Differential Equations ----------
+-- Second order linear recurrence equation
+r:= operator('r);
+r(n + 2) - 2 * r(n + 1) + r(n) = 2
+[%, r(0) = 1, r(1) = m]
+)clear properties r
+-- Second order ODE with initial conditions---solve first using Laplace
+-- transforms
+f:= operator('f);
+ode:= D(f(t), t, 2) + 4*f(t) = sin(2*t)
+map(e +-> laplace(e, t, s), %)
+-- Now, solve the ODE directly
+solve(ode, f, t = 0, [0, 0])
+-- First order linear ODE
+y:= operator('y);
+x**2 * D(y(x), x) + 3*x*y(x) = sin(x)/x
+solve(%, y, x)
+-- Nonlinear ODE
+D(y(x), x, 2) + y(x)*D(y(x), x)**3 = 0
+solve(%, y, x)
+-- A simple parametric ODE
+D(y(x, a), x) = a*y(x, a)
+solve(%, y, x)
+D(y(x), x) = a*y(x)
+solve(%, y, x)
+-- ODE with boundary conditions.  This problem has nontrivial solutions
+-- y(x) = A sin([pi/2 + n pi] x) for n an arbitrary integer.
+solve(D(y(x), x, 2) + k**2*y(x) = 0, y, x)
+-- bc(%, x = 0, y = 0, x = 1, D(y(x), x) = 0)
+-- System of two linear, constant coefficient ODEs
+x:= operator('x);
+system:= [D(x(t), t) = x(t) - y(t), D(y(t), t) = x(t) + y(t)]
+solve(system,[x,y],t)
+-- Check the answer
+-- Triangular system of two ODEs
+system:= [D(x(t), t) = x(t) * (1 + cos(t)/(2 + sin(t))), _
+          D(y(t), t) = x(t) - y(t)]
+-- Try solving this system one equation at a time
+solve(system.1, x, t)
+genericx:=C1*%.basis.1
+eval(lhs rightZero system.2,x,genericx,t)
+solve(%,y,t)
+genericy:=simplify (%.particular)+K1*(%.basis.1)
+eval(lhs rightZero system.1,x,genericx,t)
+eval(lhs rightZero system.2,[x,y],[genericx,genericy],t)
+)clear properties x y
+-- ---------- Operators ----------
+-- Linear differential operator
+DD:= operator("D") :: Operator(Expression Integer)
+evaluate(DD, e +-> D(e, x))$Operator(Expression Integer)
+L:= (DD - 1) * (DD + 2)
+g:= operator('g)
+L(f(x))
+subst(L(subst(g(y), y = x)), x = y)
+subst(L(subst(A * sin(z**2), z = x)), x = z)
+-- Truncated Taylor series operator
+T:= (f, xx, a) +-> subst((DD**0)(f(x)), x = a)/factorial(0) * (xx - a)**0 + _
+                   subst((DD**1)(f(x)), x = a)/factorial(1) * (xx - a)**1 + _
+                   subst((DD**2)(f(x)), x = a)/factorial(2) * (xx - a)**2
+T(f, x, a)
+T(g, y, b)
+Sin:= operator("sin") :: Operator(Expression Integer)
+evaluate(Sin, x +-> sin(x))$Operator(Expression Integer)
+T(Sin, z, c)
+-- ---------- Programming ----------
+-- Write a simple program to compute Legendre polynomials
+)clear properties p
+p(n, x) == 1/(2**n*factorial(n)) * D((x**2 - 1)**n, x, n)
+for i in 0..4 repeat {  output("");    output(concat(["p(", string(i), ", x) = "]));    output(p(i, x))}
+eval(p(4, x), x = 1)
+-- Now, perform the same computation using a recursive definition
+pp(0, x) == 1
+pp(1, x) == x
+pp(n, x) == ((2*n - 1)*x*pp(n - 1, x) - (n - 1)*pp(n - 2, x))/n
+for i in 0..4 repeat {   output("");   output(concat(["pp(", string(i), ", x) = "]));   output(pp(i, x))}
+)clear properties p pp
+-- ---------- Translation ----------
+-- Horner's rule---this is important for numerical algorithms
+a:= operator('a)
+sum(a(i)*x**i, i = 1..5)
+p:= factor(%)
+-- Convert the result into FORTRAN syntax
+)set fortran ints2floats off
+outputAsFortran('p = p)
+-- ---------- Boolean Logic ----------
+-- Simplify logical expressions
+true and false
+x or (not x)
+--x or y or (x and y)
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}