#% FILE antonNM6.tex

interface(echo=2);

# A simpler example than AR6 Section 9.4 Example 3, p458-460
# We set the value of Digits to 5, and try to solve the equations
 eqn1:= 10^(-8)*x+y=1;                                   eqn2:= x+y=2;
 vars:= {x,y}:            Digits:= 5:            fsolve( {eqn1,eqn2},vars );
# If we pivot on the first coefficient of  eqn2,  we get
eqn:= evalf(eqn1-eqn2/10^8);                      fsolve( {eqn,eqn2},vars );
# which is the correct answer to five significant figures.
#
# However, if we try to pivot on the first coefficient of  eqn1,  we get
eqn:= evalf(eqn2-10^8*eqn1);                     fsolve( {eqn1,eqn},vars );
# which is not the correct answer.
# Note that Maple gives the correct answer (in both cases) with the value
Digits:= 10:

# with n prescribed, e.g.
n:=3:
A:= array(symmetric,1..n,1..n);


# THIS ITEM IS IN Maple and numerical linear algebra document
#
# It can be shown that if A is a symmetric banded matrix A, with
# determinants of the submatrices centered on the diagonal all non-zero, then
#                    A =  L D transpose(L) 
# where L is a banded unit lower triangular matrix, and D is diagonal.
# We illustrate this with symmetric banded matrices where the band is
#    no longer than 5 entries.
with(linalg):
A2:= matrix([ [a11,a12],
              [a12,a22] ]):
L:= matrix([ [1,0], [L21,1] ]):      
Dg:= diag(d1,d2):      U:= transpose(L):
B:= evalm( L&*Dg&*U - A2 );
L21:= solve(  B[1,2], L21 );
# Note that this does not make sense if  d1 = 0.
#
# We repeat this for 3 by 3 matrices.
a11:= d1:                            a22:= L21^2*d1 + d2:
A3:= matrix([ [a11,a12,a13],
              [a12,a22,a23],
              [a13,a23,a33] ]):
L:= matrix([[1,0,0], [L21,1,0], [L31,L32,1] ]):
Dg:= diag(Dg,d3):    U:= transpose(L):   
B:= map( simplify,evalm(L&*Dg&*U-A3) );
# The previous values of L21,d1,d2 have forced  B[1,1] = B[1,2] = B[2,2] = 0.
ans:= solve( {B[1,3],B[2,3]}, {L31,L32} );
# Again  d1 <> 0;  and also  d2 = det(A2) <> 0.
# Note also that  d3 = det(A3)/(d1*d2) <<<computationally useless
#
# Now look at 4 by 4 banded matrices,
# where the only new values found are the values for L42 and L43 and
 L42:= a24/d2:                          assign(ans):
 a33:= L31^2*d1 + L32^2*d2 + d3:        a44:= L42^2*d2 + L43^2*d3 + d4:
 A4:= matrix([ [a11,a12,a13,  0],
               [a12,a22,a23,a24],
               [a13,a23,a33,a34],
               [0,  a24,a34,a44] ]):
 L:= matrix([ [1,0,0,0], [L21,1,0,0], [L31,L32,1,0], [0,L42,L43,1] ]):
 Dg:= diag(D,d4):    U:= transpose(L):   
B:= map( simplify,evalm(L&*Dg&*U-A4) );
 solve( B[3,4], L43 );

interface(echo=1);
