diff options
Diffstat (limited to 'src/input/wester.input.pamphlet')
| -rw-r--r-- | src/input/wester.input.pamphlet | 395 |
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} |