summaryrefslogtreecommitdiff
path: root/75/CH9
diff options
context:
space:
mode:
Diffstat (limited to '75/CH9')
-rwxr-xr-x75/CH9/EX9.1/ex_1.sce11
-rwxr-xr-x75/CH9/EX9.11/ex_11.sce27
-rwxr-xr-x75/CH9/EX9.12/ex_12.sce6
-rwxr-xr-x75/CH9/EX9.13/ex_13.sce6
-rwxr-xr-x75/CH9/EX9.14/ex_14.sce8
-rwxr-xr-x75/CH9/EX9.15/ex_15.sce13
-rwxr-xr-x75/CH9/EX9.16/ex_16.sce23
-rwxr-xr-x75/CH9/EX9.18/ex_18.sce17
-rwxr-xr-x75/CH9/EX9.19/ex_19.sce29
-rwxr-xr-x75/CH9/EX9.2/ex_2.sce25
-rwxr-xr-x75/CH9/EX9.3/ex_3.sce49
-rwxr-xr-x75/CH9/EX9.4/ex_4.sce26
-rwxr-xr-x75/CH9/EX9.5/ex_5.sce25
-rwxr-xr-x75/CH9/EX9.7/ex_7.sce27
-rwxr-xr-x75/CH9/EX9.8/ex_8.sce17
-rwxr-xr-x75/CH9/EX9.9/ex_9.sce13
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