initial commit / add all books

16 files changed, 322 insertions, 0 deletions

diff --git a/75/CH9/EX9.1/ex_1.sce b/75/CH9/EX9.1/ex_1.sce
new file mode 100755
index 000000000..a38f34fe3
--- /dev/null
+++ b/75/CH9/EX9.1/ex_1.sce
@@ -0,0 +1,11 @@
+        //    EXAMPLE 590    
+
+A = [4 1 0;1 0 -1;1 1 -4]
+[n,m] = size(A);
+
+if m<>n then
+	error('eigenvectors - matrix A is not square');
+	abort;
+end;
+
+lam = spec(A)'                   //Eigenvalues of matrix A
diff --git a/75/CH9/EX9.11/ex_11.sce b/75/CH9/EX9.11/ex_11.sce
new file mode 100755
index 000000000..84d26d659
--- /dev/null
+++ b/75/CH9/EX9.11/ex_11.sce
@@ -0,0 +1,27 @@
+        //    PG (613)
+
+A = [4 1 1;1 4 1;1 1 4]
+w1 = [0.985599 0.119573 0.119573]'
+P1 = eye() - 2*w1*w1'
+A2 = P1*A
+w2 = [0 0.996393 0.0848572]'
+P2 = eye() - 2*w2*w2'
+R = P2*A2
+Q = P1*P2
+Q*R
+
+//    A = Q * R
+
+abs(det(A))
+abs(det(Q)*det(R))
+
+//    |det(A)| = |det(Q)*det(R)| = |det(R)| = 54 (approx) 
+
+lam = spec(A)'
+lam1 = lam(1,1)
+lam2 = lam(1,2)
+lam3 = lam(1,3)
+lam1 * lam2 * lam3
+
+//    Product of eigen values also comes out to be 54
+
diff --git a/75/CH9/EX9.12/ex_12.sce b/75/CH9/EX9.12/ex_12.sce
new file mode 100755
index 000000000..94c3559bb
--- /dev/null
+++ b/75/CH9/EX9.12/ex_12.sce
@@ -0,0 +1,6 @@
+        //    PG (617)
+
+A = [1 3 4;3 1 2;4 2 1]
+w2 = [0 2/sqrt(5) 1/sqrt(5)]'
+P1 = eye() - 2*w2*w2'
+T = P1' * A * P1    //    Tridiagonal matrix
diff --git a/75/CH9/EX9.13/ex_13.sce b/75/CH9/EX9.13/ex_13.sce
new file mode 100755
index 000000000..3e7638bd7
--- /dev/null
+++ b/75/CH9/EX9.13/ex_13.sce
@@ -0,0 +1,6 @@
+        //    PG (619)
+
+x = %pi/4
+R = [cos(x) 0 sin(x);0 1 0;-sin(x) 0 cos(x)]
+
+//    Planner Rotation Orthogonal Matrix
\ No newline at end of file
diff --git a/75/CH9/EX9.14/ex_14.sce b/75/CH9/EX9.14/ex_14.sce
new file mode 100755
index 000000000..79f91ae2f
--- /dev/null
+++ b/75/CH9/EX9.14/ex_14.sce
@@ -0,0 +1,8 @@
+        //    PG (620)
+
+T = [2 1 0 0 0 0;1 2 1 0 0 0;0 1 2 1 0 0;0 0 1 2 1 0;0 0 0 1 2 1;0 0 0 0 1 2]
+lam = spec(T)'
+lam1 = lam(1,1)
+B = [2-lam1 1 0 0 0 0;1 2-lam1 1 0 0 0;0 1 2-lam1 1 0 0;0 0 1 2-lam1 1 0;0 0 0 1 2-lam1 1;0 0 0 0 1 2]
+f0 = abs(det(B))
+f1 = 2-lam1
diff --git a/75/CH9/EX9.15/ex_15.sce b/75/CH9/EX9.15/ex_15.sce
new file mode 100755
index 000000000..8ce1790a8
--- /dev/null
+++ b/75/CH9/EX9.15/ex_15.sce
@@ -0,0 +1,13 @@
+        //    PG (621)
+
+//    For the previous example, consider the sequence f0, f1....f6
+
+//    For lam = 3,
+
+//        (f0,.....f6) = (1,-1,0,1,-1,0,1)
+
+//    The corresponding sequence of signs is 
+
+//        (+,-,+,+,-,+,+)
+
+//    and s(3) = 2
\ No newline at end of file
diff --git a/75/CH9/EX9.16/ex_16.sce b/75/CH9/EX9.16/ex_16.sce
new file mode 100755
index 000000000..7d2024183
--- /dev/null
+++ b/75/CH9/EX9.16/ex_16.sce
@@ -0,0 +1,23 @@
+        //    PG (624)
+
+A1 = [2 1 0;1 3 1;0 1 4]
+lam = spec(A1)'
+[Q1,R1] = qr(A1);
+A2 = R1 * Q1
+[Q2,R2] = qr(A2);
+A3 = R2 * Q2
+[Q3,R3] = qr(A3);
+A4 = R3 * Q3
+[Q4,R4] = qr(A4);
+A5 = R4 * Q4
+[Q5,R5] = qr(A5);
+A6 = R5 * Q5
+[Q6,R6] = qr(A6);
+A7 = R6 * Q6
+[Q7,R7] = qr(A7);
+A8 = R7 * Q7
+[Q8,R8] = qr(A8);
+A9 = R8 * Q8
+[Q9,R9] = qr(A9);
+A10 = R9 * Q9
+[Q10,R10] = qr(A10);
\ No newline at end of file
diff --git a/75/CH9/EX9.18/ex_18.sce b/75/CH9/EX9.18/ex_18.sce
new file mode 100755
index 000000000..db558f531
--- /dev/null
+++ b/75/CH9/EX9.18/ex_18.sce
@@ -0,0 +1,17 @@
+        //    PG (631)
+
+A = [2 1 0;1 3 1;0 1 4]
+lam = spec(A)
+[L,U] = lu(A)
+y1 = [1 1 1]'
+w1 = [3385.2 -2477.3 908.20]'
+z1 = [w1/norm(w1,'inf')]'
+w2 = [20345 -14894 5451.9]'
+z2 = [w2/norm(w2,'inf')]'
+z3 = z2
+
+//    The true answer is
+
+x3 = [1 1-sqrt(3) 2-sqrt(3)]'
+
+//    z2 equals x3 to within the limits of rounding error accumulations.
\ No newline at end of file
diff --git a/75/CH9/EX9.19/ex_19.sce b/75/CH9/EX9.19/ex_19.sce
new file mode 100755
index 000000000..6b88a7ee4
--- /dev/null
+++ b/75/CH9/EX9.19/ex_19.sce
@@ -0,0 +1,29 @@
+        //    PG (633)
+
+A = [2 1 0;1 3 1;0 1 4]
+lam = spec(A)
+[L,U] = lu(A)
+y1 = [1 1 1]'
+w1 = [3385.2 -2477.3 908.20]'
+z1 = [w1/norm(w1,'inf')]'
+w2 = [20345 -14894 5451.9]'
+z2 = [w2/norm(w2,'inf')]'
+z3 = z2
+
+//    The true answer is
+
+x3 = [1 1-sqrt(3) 2-sqrt(3)]'
+
+//    z2 equals x3 to within the limits of rounding error accumulations.
+
+//    Consider lam = 1.2679
+
+//    0.7321*x1 + x2 = 0
+//    x1 + 1.7321*x2 + x3 = 0
+//    Taking x1= 1.0, we have the approximate eigenvector 
+
+//                    x = [1.0000 -0.73210 0.26807]
+
+
+//    Compared with the true answer obtained above, this is a slightly poorer
+//    result obtained by inverse iteration.
\ No newline at end of file
diff --git a/75/CH9/EX9.2/ex_2.sce b/75/CH9/EX9.2/ex_2.sce
new file mode 100755
index 000000000..2d28eadd5
--- /dev/null
+++ b/75/CH9/EX9.2/ex_2.sce
@@ -0,0 +1,25 @@
+        // PG 591
+
+n = 4
+A = [4 1 0 0;1 4 1 0;0 1 4 1;0 0 1 4]
+lam = spec(A)
+
+//    Since A is symmetric, all eigen values are real.
+//    The radii are all 1 or 2.
+//    The centers of all the circles are 4.
+//    All eigen values must all lie in the interval [2,6]
+//    Since the eigen values of inv(A) are the reciprocals of those of A,
+//    1/6 <= mu <= 1/2
+
+//    Let inv(A) = B
+
+B=inv(A);
+norm(B,2)
+n
+i = 1:n;
+j = 1:n;
+
+    //    for j~i
+    //        r = sum(abs(B(i,j)))
+
+//    norm(B,2) = r(B) <= o.5
\ No newline at end of file
diff --git a/75/CH9/EX9.3/ex_3.sce b/75/CH9/EX9.3/ex_3.sce
new file mode 100755
index 000000000..865d25af7
--- /dev/null
+++ b/75/CH9/EX9.3/ex_3.sce
@@ -0,0 +1,49 @@
+        //    PG 593
+
+disp("Consider Hilbert matrix of order three")
+
+n=3;         //    Order of the matrix
+A=zeros(n,n);//    a symmetric positive definite real or complex matrix.
+for i=1:n    //    Initializing 'for' loop
+    for j=1:n
+        A(i,j)=1/(i+j-1);
+    end
+end        //End of 'for' loop
+A
+
+[n,m] = size(A)
+
+if m<>n then
+	error('eigenvectors - matrix A is not square');
+	abort;
+end;
+
+lam = spec(A)'                   //Eigenvalues of matrix A
+
+lam1 = lam(1,1)
+lam2 = lam(1,2)
+lam3 = lam(1,3)
+
+    //    Rounding off to 4 decimal places
+
+A = A*10^4;
+    A = int(A);
+    A = A*10^(-4);
+    disp(A)        //    Final Solution
+
+lamr = spec(A)'
+
+lamr1 = lamr(1,1)
+lamr2 = lamr(1,2)
+lamr3 = lamr(1,3)
+
+    //    Errors
+
+lam - lamr
+
+    //    Relative Errors
+
+R1 = (lam1-lamr1)/lam1
+R2 = (lam2-lamr2)/lam2
+R3 = (lam3-lamr3)/lam3
+
diff --git a/75/CH9/EX9.4/ex_4.sce b/75/CH9/EX9.4/ex_4.sce
new file mode 100755
index 000000000..f7ac8f3f1
--- /dev/null
+++ b/75/CH9/EX9.4/ex_4.sce
@@ -0,0 +1,26 @@
+        //    PG 594
+
+A = [101 -90;110 -98]
+[n,m] = size(A)
+
+if m<>n then
+	error('eigenvectors - matrix A is not square');
+	abort;
+end;
+
+lam = spec(A)'                   //Eigenvalues of matrix A
+
+
+    //    A+E = [101-e -90-e;110 -98]
+    //    Let e = 0.001
+e = 0.001;
+    //    Let A+E = D
+D = [101-e -90-e;110 -98]
+
+[n,m] = size(D)
+
+if m<>n then
+	error('eigenvectors - matrix D is not square');
+	abort;
+end;
+lam = spec(D)'                   //Eigenvalues of matrix A
diff --git a/75/CH9/EX9.5/ex_5.sce b/75/CH9/EX9.5/ex_5.sce
new file mode 100755
index 000000000..6982ab9db
--- /dev/null
+++ b/75/CH9/EX9.5/ex_5.sce
@@ -0,0 +1,25 @@
+        //    PG 599
+
+    //    e = 0.001
+    //    From earlier example : 
+    //    eigen values of matrix A are 1 and 2. So,..
+
+    //    inv(P)*A*P = [1 0;0 2]
+
+A = [101 -90;110 -98]
+B = [-1 -1;0 0]
+    //    From the above equation, we get:
+
+P = [9/sqrt(181) -10/sqrt(221);10/sqrt(181) -11/sqrt(221)]
+inv(P)
+K = norm(P)*norm(inv(P))        //    K is condition number
+u1 = P(:,1)
+u2 = P(:,2)
+Q = inv(P)
+R = Q'
+w1 = R(:,1)
+w2 = R(:,2)
+s1 = 1/norm(w1,2)
+norm(B)
+
+//    abs(lam1(e) - lam1) <= sqrt(2)*e/0.005 + O(e^2) = 283*e + O(e^2)
diff --git a/75/CH9/EX9.7/ex_7.sce b/75/CH9/EX9.7/ex_7.sce
new file mode 100755
index 000000000..6bca22d5e
--- /dev/null
+++ b/75/CH9/EX9.7/ex_7.sce
@@ -0,0 +1,27 @@
+        //    (PG 607)
+
+A = [1 2 3;2 3 4;3 4 5]
+lam = spec(A)'
+lam1 = lam(1,3)
+lam2 = lam(1,1)
+lam3 = lam(1,2)
+
+    //    Theoretical ratio of convergence
+
+lam2/lam1
+
+b = 0.5*(lam2+lam3)
+B = A-b*eye(3,3)
+
+    //    Eigen values of A-bI are:
+
+lamb = spec(B)'
+lamb1 = lamb(1,3)
+lamb2 = lamb(1,2)
+lamb3 = lamb(1,1)
+
+    //    Ratio of convergence for the power method applied to A-bI will be:
+
+lamb2/lamb1
+
+    //    This is less than half the magnitude of the original ratio.
\ No newline at end of file
diff --git a/75/CH9/EX9.8/ex_8.sce b/75/CH9/EX9.8/ex_8.sce
new file mode 100755
index 000000000..25013bc14
--- /dev/null
+++ b/75/CH9/EX9.8/ex_8.sce
@@ -0,0 +1,17 @@
+        //    PG (608)
+
+A = [1 2 3;2 3 4;3 4 5]
+lam = spec(A)'        //    Eigen values of A
+lam1 = lam(1,3)
+lam2 = lam(1,1)
+lam3 = lam(1,2)
+
+//    Theoretical ratio of convergence
+
+lam2/lam1
+
+//    After extrapolating, we get
+        lame1 = 9.6234814
+
+//    Error:
+lam1-lame1
\ No newline at end of file
diff --git a/75/CH9/EX9.9/ex_9.sce b/75/CH9/EX9.9/ex_9.sce
new file mode 100755
index 000000000..92de3a6a6
--- /dev/null
+++ b/75/CH9/EX9.9/ex_9.sce
@@ -0,0 +1,13 @@
+        //    PG (610)
+
+w = [1/3 2/3 2/3]'
+w1 = w(1,1)
+w2 = w(2,1)
+w3 = w(3,1)
+
+U = [1-2*abs(w1)^2 -2*w1*w2' -2*w1*w3';-2*w1'*w2 1-2*abs(w2)^2 -2*w2*w3';-2*w1'*w3 -2*w2'*w3 1-2*abs(w3)^2]
+U
+inv(U)
+//    U = inv(U)------Hence, U is Hermitian
+U*U
+//    U*U = I---------Hence, U is orthogonal