# file maple6

dot:= proc(fx,gx,lims)        int( expand(fx*gx),x=lims )        end:
# Hill Section 6.1 Legendre Polynomials
with(orthopoly):
# The first three Legendre polynomials (to within a constant factor) are
P(0,x), P(1,x), P(2,x);
# Check their orthogonality using the dot product given in Section 4.1:
dot(P(0,x),P(1,x),-1..1), dot(P(0,x),P(2,x),-1..1), dot(P(1,x),P(2,x),-1..1);
quit

# Hill Section 6.1 Example on p400-402.
# (See above for the definition of the procedure  dot )
one_one:= dot(1,1,0..1):     one_x:= dot(1,x,0..1):     x_x:= dot(x,x,0..1):
with(linalg):
A:= matrix([ [one_one,one_x], [one_x,x_x] ]);
# The line of best fit for the function  sqrt(x)  over the range 0..1 is
b_m:= linsolve( A, vector( [dot(1,sqrt(x),0..1),dot(x,sqrt(x),0..1)] ) ):
y = b_m[1] + b_m[2] * x;
# Alternatively, we find the orthogonal polynomials (using Gram Schmidt):
v0:= 1:                                              v1:= x - one_x * v0;
# Since  v0, v1  are orthogonal, the line of best fit for any function f is
best_fit:= proc(f)
    dot(f,v0,0..1)/dot(v0,v0,0..1)*v0 + dot(f,v1,0..1)/dot(v1,v1,0..1)*v1
end:
# So the line of best fit for the function  sqrt(x)  over the range 0..1 is
best_fit( sqrt(x) );
quit

with(linalg):
pseudoinverse:= proc(A)
    local U,V,Sigma_plus,singvals,i:
    Sigma_plus:= array( sparse,1..coldim(A),1..rowdim(A) ):
    singvals:= evalf(Svd(A,V,U));
    for i from 1 to vectdim(singvals) do
        if singvals[i]>0 then  Sigma_plus[i,i]:= evalf( 1/singvals[i] )  fi:
    od:
    map( evalf, evalm(U&*Sigma_plus&*transpose(V)) );
end:
#
# Hill Section 4.3, Example 5, p264/5
A:= matrix([ [1,3], [2,1], [0,1], [-2,1] ]):
# The matrix of the perpendicular projection function is P where
P:= evalm( A&*pseudoinverse(A) ):
map( evalf,evalm(33*P) );   # The previous answer had 33 in the denominators

# Hill Section 4.3, Example 6, p265/6
A:= matrix([ [40,1], [45,1], [50,1], [55,1], [60,1] ]):
y:= vector( [482,467,452,433,421] ):
# The least-squares solution to the equation  A x = y  is given by
map( evalf,evalm(pseudoinverse(A)&*y) );
quit

# Hill Section 6.3   Examples of Iterative Methods pp 414-418
# Example 1: Solve  cos(x) = x
x:= 0.7:             # initial value
for i from 1 to 4 do        x:= cos(x):        print(x);        od;

# Example 2: Find the square root of 2.
x:= 1.5:
for i from 1 to 3 do     x:= evalf( (x+2/x)/2 ):      print(x);      od;
    local n,DD,L,U,invLplusD,x,i,j:

with(linalg):
# Example 3: Jacobi's method for solving the linear equations  A x = b  where
A:= matrix([ [2,-1],
             [-1,2] ]):
b:= vector([3,6]):                                         m:= 6:
x:= vector([0,0]):
for i from 1 to m do
    newx1:= ( b[1] - A[1,2]*x[2] )/A[1,1]:
    newx2:= ( b[2] - A[2,1]*x[1] )/A[2,2]:
    x[1]:= evalf(newx1):            x[2]:= evalf(newx2):            print(x);
od:

# Example 4: Gauss-Seidel method    (Note also how much neater the code is)
x:= vector([0,0]):
for i from 1 to m do
    x[1]:= evalf( (b[1]-A[1,2]*x[2])/A[1,1] ):
    x[2]:= evalf( (b[2]-A[2,1]*x[1])/A[2,2] ):         print(x);
od:

# The matrix form of the Gauss-Seidel method (written as a procedure) is
GSsolve:= proc(A,b,m)     # m is the number of iterations
    local n,DD,L,U,invLplusD,x,i,j:
    n:= rowdim(A):
    L:= array( sparse,1..n,1..n ):      U:= array( sparse,1..n,1..n ):
    DD:= array( sparse,1..n,1..n ):
    for i from 1 to n do
        DD[i,i]:= A[i,i]:
        for j from 1 to i-1 do    L[i,j]:= A[i,j]:    U[j,i]:= A[j,i]:    od:
    od:
    invLplusD:= inverse(L+DD):
    x:= vector( [0$n] ):     #  0$n  denotes a sequence of n 0's
    for i from 1 to m do
        print(x);
        x:= map( evalf,evalm( invLplusD &* (b-U&*x) ) ):
    od;
end:
#
# So Example 4 (above) becomes
GSsolve(A,b,6);
quit

# Hill Section 6.4 Example 5 p424
with(linalg): 
A:= matrix([ [1,1], [0,2] ]):
singularvals(A);
# Note that, before Release 3, Maple did not simplify the answer to this
norm(A,infinity);
norm(A,2);
# Before Release 2, Maple could not decide that this was the larger one
evalf(");
quit

# Hill Section 6.5 Examples 7-10, pp 435-437
with(linalg):
alpha1:=vector([1,0,0]):  alpha2:=vector([1,1,0]):  alpha3:=vector([1,1,1]):
ALPHA:= augment( alpha1,alpha2,alpha3 );
# Since the field is  Z2,  use the function  modp( - , 2)
# Example 7
# {alpha1,alpha2,alpha3}  is a basis if and only if  ALPHA  is invertible:
modp( det(ALPHA),2 );
# Example 8
v:= map( modp, evalm(alpha1+alpha3), 2 );
# The coordinate vector of v relative to the basis ALPHA is  vAlpha  where
vAlpha:= map( modp, linsolve(ALPHA,v), 2);
# Hill's `Remarkable Fact 2' p435:    Note that  v <> 0,  but
modp( dotprod(v,v),2 );
# Example 9    (compare Section 3.8)
# Find the matrix of the linear transformation  T  where
T:= proc(x)    vector( [x[1]+x[3], x[2]+x[3] ])        end:
#    relative to the bases  ALPHA and BETA  where
beta1:= vector( [0,1] ):       beta2:= vector( [1,1] ):
BETA:= augment( beta1,beta2 );
T_BetaAlpha:= augment( map( modp, linsolve(BETA,T(alpha1)), 2),
                       map( modp, linsolve(BETA,T(alpha2)), 2),
                       map( modp, linsolve(BETA,T(alpha3)), 2) );
# Example 10
# The coordinates of T(v) relative to the basis  BETA  is  Tv_Beta  where
Tv_Beta:= map( modp, linsolve(BETA,T(v)), 2);
# By Theorem 6.97, this should be the same as the following
map( modp, evalm( T_BetaAlpha &* vAlpha ), 2);

det(ALPHA) mod 2;
quit

