diff --git a/ENUME/projectA/QRNoShifts.asv b/ENUME/projectA/QRNoShifts.asv new file mode 100644 index 00000000..4ee40386 --- /dev/null +++ b/ENUME/projectA/QRNoShifts.asv @@ -0,0 +1,62 @@ +function [eigenValues, whichIterationAreWeOn, Matrix] = QRNoShifts(Matrix) + [whichIterationAreWeOn, threshold] = initializeValues(); + while threshold > 1e-6 + [Q, R] = gramSchmidtAlgorithm(Matrix); + Matrix = R * Q; + whichIterationAreWeOn = whichIterationAreWeOn + 1; + + % iterate until all non-diagonal elements are below the threshold + nondiag = Matrix - diag(diag(Matrix)); % first diag converts Matrix + % into vector consisting of values on the diagonal of matrix, + % second diag converts this vector into matrix with zeros on + % everything except diagonal + threshold = max(abs(nondiag)); + end + + % convert eigenvalue matrix to vector + eigenValues = diag(Matrix)'; + disp("Iterations it took: "); + disp(whichIterationAreWeOn); + disp("Final matrix: "); + disp(Matrix); +end + +function [whichIterationAreWeOn, threshold] = initializeValues() + whichIterationAreWeOn = 0; + threshold = inf; +end + +% performs QR or QRdash decomposition of a matrix +function [Q, R] = gramSchmidtAlgorithm(Matrix) + [columns, Q, R, d] = initializeGramSchmid(Matrix); + [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d); + [Q, R] = normalizeColumns(columns, Q, R); +end + +function [columns, Q, R, d] = initializeGramSchmid(Matrix) + % We start with empty matrices + [rows, columns] = size(Matrix); + Q = zeros(rows, columns); + R = zeros(columns, columns); + d = zeros(1, columns); +end + +function [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d) + for i = 1 : columns + Q(:, i) = Matrix(:, i); + R(i, i) = 1; + d(i) = Q(:, i)' * Q(:, i); + for i2 = i + 1 : columns + R(i, i2) = (Q(:, i)' * Matrix(:, i2)) / d(i); + Matrix(:, i2) = Matrix(:, i2) - R(i, i2) * Q(:, i); + end + end +end + +function [Q, R] = normalizeColumns(columns, Q, R) + for i = 1 : columns + dd = norm(Q(:, i)); + Q(:, i) = Q(:, i) / dd; + R(i, i:columns) = R(i, i : columns) * dd; + end +end \ No newline at end of file diff --git a/ENUME/projectA/QRNoShifts.m b/ENUME/projectA/QRNoShifts.m new file mode 100644 index 00000000..8669d237 --- /dev/null +++ b/ENUME/projectA/QRNoShifts.m @@ -0,0 +1,64 @@ +function [eigenValues, whichIterationAreWeOn, Matrix] = QRNoShifts(Matrix) + [whichIterationAreWeOn, threshold] = initializeValues(); + while threshold > 1e-6 + [Q, R] = gramSchmidtAlgorithm(Matrix); + Matrix = R * Q; + whichIterationAreWeOn = whichIterationAreWeOn + 1; + + % iterate until all non-diagonal elements are below the threshold + matrixWithoutDiagonal = Matrix - diag(diag(Matrix)); % first diag converts Matrix + % into vector consisting of values on the diagonal of matrix, + % second diag converts this vector into matrix with zeros on + % everything except diagonal + % If we substract it from Matrix we get original Matrix with zeros + % on a diagonal + threshold = max(abs(matrixWithoutDiagonal)); + end + + % convert eigenvalue matrix to vector + eigenValues = diag(Matrix)'; + disp("Iterations it took: "); + disp(whichIterationAreWeOn); + disp("Final matrix: "); + disp(Matrix); +end + +function [whichIterationAreWeOn, threshold] = initializeValues() + whichIterationAreWeOn = 0; + threshold = inf; +end + +% performs QR or QRdash decomposition of a matrix +function [Q, R] = gramSchmidtAlgorithm(Matrix) + [columns, Q, R, d] = initializeGramSchmid(Matrix); + [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d); + [Q, R] = normalizeColumns(columns, Q, R); +end + +function [columns, Q, R, d] = initializeGramSchmid(Matrix) + % We start with empty matrices + [rows, columns] = size(Matrix); + Q = zeros(rows, columns); + R = zeros(columns, columns); + d = zeros(1, columns); +end + +function [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d) + for i = 1 : columns + Q(:, i) = Matrix(:, i); + R(i, i) = 1; + d(i) = Q(:, i)' * Q(:, i); + for i2 = i + 1 : columns + R(i, i2) = (Q(:, i)' * Matrix(:, i2)) / d(i); + Matrix(:, i2) = Matrix(:, i2) - R(i, i2) * Q(:, i); + end + end +end + +function [Q, R] = normalizeColumns(columns, Q, R) + for i = 1 : columns + dd = norm(Q(:, i)); + Q(:, i) = Q(:, i) / dd; + R(i, i:columns) = R(i, i : columns) * dd; + end +end \ No newline at end of file diff --git a/ENUME/projectA/QRShifts.asv b/ENUME/projectA/QRShifts.asv new file mode 100644 index 00000000..4691f5e7 --- /dev/null +++ b/ENUME/projectA/QRShifts.asv @@ -0,0 +1,121 @@ +% finds the eigenValues of a matrix using QR with shifts +function [eigenValues, whatIterationAreWeOn, Matrix] = QRShifts(Matrix) + [eigenValues, whatIterationAreWeOn, matrixSize, minThreshold] = initiateValues(Matrix); + [Matrix, whatIterationAreWeOn, eigenValues] = QRShiftLoop(matrixSize, Matrix, eigenValues, minThreshold, whatIterationAreWeOn); +end + +function [eigenValues, whatIterationAreWeOn, matrixSize, minThreshold] = initiateValues(Matrix) + eigenValues = double.empty(1, 0); + whatIterationAreWeOn = 0; + matrixSize = size(Matrix, 1); + minThreshold = 1e-6; +end + +function [finalMatrix, whatIterationAreWeOn, eigenValues] = QRShiftLoop(matrixSize, Matrix, eigenValues, minThreshold, whatIterationAreWeOn) + while matrixSize >= 2 + flag = 0; + [Matrix, matrixSize, whatIterationAreWeOn, eigenValues] = findEigenValue(Matrix, matrixSize, whatIterationAreWeOn, eigenValues, minThreshold, flag); + [finalMatrix, matrixSize, Matrix] = deflateMatrix(Matrix, matrixSize); + end + eigenValues(size(eigenValues, 2) + 1) = Matrix(1, 1); +end + +function [Matrix, matrixSize, whatIterationAreWeOn, eigenValues] = findEigenValue(Matrix, matrixSize, whatIterationAreWeOn, eigenValues, minThreshold, flag) + while flag == 0 + eigenValueCorner = getEigenValueFromCorner(Matrix, matrixSize); + [Matrix, whatIterationAreWeOn] = shiftAndIterate(matrixSize, Matrix, eigenValueCorner, whatIterationAreWeOn); + [flag, eigenValues] = thresholdBreached(Matrix, matrixSize, minThreshold, eigenValues); + end +end + +function eigenValueCorner = getEigenValueFromCorner(Matrix, matrixSize) + corner = Matrix((matrixSize - 1) : matrixSize, (matrixSize - 1) : matrixSize); + % get 2 x 2 corner of the matrix + eigenValueCorner = solveCharactersticEquation(corner); +end + +function [Matrix, whatIterationAreWeOn] = shiftAndIterate(matrixSize, Matrix, eigenValueCorner, whatIterationAreWeOn) + % shift and iterate algorithm + identityMatrix = eye(matrixSize); + Matrix = Matrix - identityMatrix * eigenValueCorner; + [Q, R] = gramSchmidtAlgorithm(Matrix); + Matrix = R * Q + identityMatrix * eigenValueCorner; + whatIterationAreWeOn = whatIterationAreWeOn + 1; +end + +function [finalMatrix, matrixSize, Matrix] = deflateMatrix(Matrix, matrixSize) + finalMatrix(1 : matrixSize, 1 : matrixSize) = Matrix; + matrixSize = matrixSize - 1; + Matrix = Matrix(1:matrixSize, 1:matrixSize); +end + +function [flag, eigenValues] = thresholdBreached(Matrix, matrixSize, minThreshold, eigenValues) + flag = 0; + threshold = max(abs(Matrix(matrixSize, 1:(matrixSize - 1)))); + + % once we zero the row (or rather get close enough to zero) exit loop + if (threshold <= minThreshold) + eigenValues(size(eigenValues, 2) + 1) = Matrix(matrixSize, matrixSize); + flag = 1; + end +end + +% finds the eigenvalue of a 2x2 matrix that is closer to the lower right corner +function eigen = solveCharactersticEquation(Matrix) + [eigen1, eigen2] = calculateZeros(Matrix); + +end + +function eigen = valueCloserToLowerRightCorner(eigen1, eigen) + if abs(Matrix(4) - eigen1) < abs(Matrix(4) - eigen2) + eigen = eigen1; + else + eigen = eigen2; + end +end + +function [zeroOne, zeroTwo] = calculateZeros(Matrix) + b = Matrix(1) + Matrix(4); + ac = Matrix(1) * Matrix(4) - Matrix(2) * Matrix(3); + delta = (b) ^ 2 - 4 * ac; + % get delta of quadratic equation + squareRootDelta = sqrt(delta); + % square delta + zeroOne = (b - squareRootDelta) / 2; + zeroTwo = (b + squareRootDelta) / 2; +end + +% performs QR or QRdash decomposition of a matrix +function [Q, R] = gramSchmidtAlgorithm(Matrix) + [columns, Q, R, d] = initializeGramSchmid(Matrix); + [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d); + [Q, R] = normalizeColumns(columns, Q, R); +end + +function [columns, Q, R, d] = initializeGramSchmid(Matrix) + % We start with empty matrices + [rows, columns] = size(Matrix); + Q = zeros(rows, columns); + R = zeros(columns, columns); + d = zeros(1, columns); +end + +function [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d) + for i = 1 : columns + Q(:, i) = Matrix(:, i); + R(i, i) = 1; + d(i) = Q(:, i)' * Q(:, i); + for i2 = i + 1 : columns + R(i, i2) = (Q(:, i)' * Matrix(:, i2)) / d(i); + Matrix(:, i2) = Matrix(:, i2) - R(i, i2) * Q(:, i); + end + end +end + +function [Q, R] = normalizeColumns(columns, Q, R) + for i = 1 : columns + dd = norm(Q(:, i)); + Q(:, i) = Q(:, i) / dd; + R(i, i:columns) = R(i, i : columns) * dd; + end +end \ No newline at end of file diff --git a/ENUME/projectA/QRShifts.m b/ENUME/projectA/QRShifts.m new file mode 100644 index 00000000..feec8a68 --- /dev/null +++ b/ENUME/projectA/QRShifts.m @@ -0,0 +1,132 @@ +% finds the eigenValues of a matrix using QR with shifts +function [eigenValues, whatIterationAreWeOn, Matrix] = QRShifts(Matrix) + eigenFromMatlab = eig(Matrix); + [eigenValues, whatIterationAreWeOn, matrixSize, minThreshold] = initiateValues(Matrix); + [Matrix, whatIterationAreWeOn, eigenValues] = QRShiftLoop(matrixSize, Matrix, eigenValues, minThreshold, whatIterationAreWeOn); + dispResults(Matrix, whatIterationAreWeOn, eigenFromMatlab); +end + +function [eigenValues, whatIterationAreWeOn, matrixSize, minThreshold] = initiateValues(Matrix) + eigenValues = double.empty(1, 0); + whatIterationAreWeOn = 0; + matrixSize = size(Matrix, 1); + minThreshold = 1e-6; +end + +function [Matrix, whatIterationAreWeOn, eigenValues] = QRShiftLoop(matrixSize, Matrix, eigenValues, minThreshold, whatIterationAreWeOn) + while matrixSize >= 2 + flag = 0; + [Matrix, matrixSize, whatIterationAreWeOn, eigenValues] = findEigenValue(Matrix, matrixSize, whatIterationAreWeOn, eigenValues, minThreshold, flag); + [matrixSize, Matrix] = deflateMatrix(Matrix, matrixSize); + end + eigenValues(size(eigenValues, 2) + 1) = Matrix(1, 1); +end + +function [Matrix, matrixSize, whatIterationAreWeOn, eigenValues] = findEigenValue(Matrix, matrixSize, whatIterationAreWeOn, eigenValues, minThreshold, flag) + while flag == 0 + eigenValueCorner = getEigenValueFromCorner(Matrix, matrixSize); + [Matrix, whatIterationAreWeOn] = shiftAndIterate(matrixSize, Matrix, eigenValueCorner, whatIterationAreWeOn); + [flag, eigenValues] = thresholdBreached(Matrix, matrixSize, minThreshold, eigenValues); + end +end + +function eigenValueCorner = getEigenValueFromCorner(Matrix, matrixSize) + corner = Matrix((matrixSize - 1) : matrixSize, (matrixSize - 1) : matrixSize); + % get 2 x 2 corner of the matrix + eigenValueCorner = solveCharactersticEquation(corner); +end + +function [Matrix, whatIterationAreWeOn] = shiftAndIterate(matrixSize, Matrix, eigenValueCorner, whatIterationAreWeOn) + % shift and iterate algorithm + identityMatrix = eye(matrixSize); + Matrix = Matrix - identityMatrix * eigenValueCorner; + [Q, R] = gramSchmidtAlgorithm(Matrix); + Matrix = R * Q + identityMatrix * eigenValueCorner; + whatIterationAreWeOn = whatIterationAreWeOn + 1; +end + +function [matrixSize, Matrix] = deflateMatrix(Matrix, matrixSize) + matrixSize = matrixSize - 1; + Matrix = Matrix(1:matrixSize, 1:matrixSize); +end + +function [flag, eigenValues] = thresholdBreached(Matrix, matrixSize, minThreshold, eigenValues) + flag = 0; + threshold = max(abs(Matrix(matrixSize, 1:(matrixSize - 1)))); + + % once we zero the row (or rather get close enough to zero) exit loop + if (threshold <= minThreshold) + eigenValues(size(eigenValues, 2) + 1) = Matrix(matrixSize, matrixSize); + flag = 1; + end +end + +% finds the eigenvalue of a 2x2 matrix that is closer to the lower right corner +function eigen = solveCharactersticEquation(Matrix) + [eigenOne, eigenTwo] = calculateZeros(Matrix); + eigen = valueCloserToLowerRightCorner(eigenOne, eigenTwo, Matrix); +end + +function eigen = valueCloserToLowerRightCorner(eigenOne, eigenTwo, Matrix) + if abs(Matrix(4) - eigenOne) < abs(Matrix(4) - eigenTwo) + eigen = eigenOne; + else + eigen = eigenTwo; + end +end + +function [zeroOne, zeroTwo] = calculateZeros(Matrix) + b = Matrix(1) + Matrix(4); + ac = Matrix(1) * Matrix(4) - Matrix(2) * Matrix(3); + delta = (b) ^ 2 - 4 * ac; + % get delta of quadratic equation + squareRootDelta = sqrt(delta); + % square delta + zeroOne = (b - squareRootDelta) / 2; + zeroTwo = (b + squareRootDelta) / 2; +end + +% performs QR or QRdash decomposition of a matrix +function [Q, R] = gramSchmidtAlgorithm(Matrix) + [columns, Q, R, d] = initializeGramSchmid(Matrix); + [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d); + [Q, R] = normalizeColumns(columns, Q, R); +end + +function [columns, Q, R, d] = initializeGramSchmid(Matrix) + % We start with empty matrices + [rows, columns] = size(Matrix); + Q = zeros(rows, columns); + R = zeros(columns, columns); + d = zeros(1, columns); +end + +function [Q, R] = factorizeColumnsOfQ(columns, Q, Matrix, R, d) + for i = 1 : columns + Q(:, i) = Matrix(:, i); + R(i, i) = 1; + d(i) = Q(:, i)' * Q(:, i); + for i2 = i + 1 : columns + R(i, i2) = (Q(:, i)' * Matrix(:, i2)) / d(i); + Matrix(:, i2) = Matrix(:, i2) - R(i, i2) * Q(:, i); + end + end +end + +function [Q, R] = normalizeColumns(columns, Q, R) + for i = 1 : columns + dd = norm(Q(:, i)); + Q(:, i) = Q(:, i) / dd; + R(i, i:columns) = R(i, i : columns) * dd; + end +end + +function dispResults(Matrix, whatIterationAreWeOn, eigenFromMatlab) + disp("Final matrix: "); + disp(Matrix); + disp("Number of iterations: "); + disp(whatIterationAreWeOn); + disp("Eigen values from eig(Matrix):"); + disp(eigenFromMatlab); + disp("Our eigen values: "); +end \ No newline at end of file diff --git a/ENUME/projectA/matrix4.m b/ENUME/projectA/matrix4.m new file mode 100644 index 00000000..902705cf --- /dev/null +++ b/ENUME/projectA/matrix4.m @@ -0,0 +1,7 @@ +function A = matrix4() + A = [1, 1, 7, 5, 2; + 1, 8, 5, 4, 4; + 7, 5, 0, 8, 8; + 5, 4, 8, 0, 8; + 2, 4, 8, 8, 1]; +end \ No newline at end of file diff --git a/ENUME/projectA/projectA.aux b/ENUME/projectA/projectA.aux index bd049376..170f05fa 100644 --- a/ENUME/projectA/projectA.aux +++ b/ENUME/projectA/projectA.aux @@ -76,6 +76,7 @@ \@writefile{toc}{\contentsline {section}{\numberline {4.1}Problem}{34}{section.4.1}\protected@file@percent } \@writefile{toc}{\contentsline {section}{\numberline {4.2}Theoretical introduction}{34}{section.4.2}\protected@file@percent } \@writefile{toc}{\contentsline {subsection}{\numberline {4.2.1}Eigenvalues}{34}{subsection.4.2.1}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {4.2.2}QR method for finding eigenvalues}{34}{subsection.4.2.2}\protected@file@percent } \@writefile{toc}{\contentsline {section}{\numberline {4.3}Solution}{35}{section.4.3}\protected@file@percent } \@writefile{toc}{\contentsline {section}{\numberline {4.4}Discussion of the result}{35}{section.4.4}\protected@file@percent } \@writefile{toc}{\contentsline {chapter}{\numberline {5}Code appendix}{36}{chapter.5}\protected@file@percent } diff --git a/ENUME/projectA/projectA.fdb_latexmk b/ENUME/projectA/projectA.fdb_latexmk index 8481e777..4cc4e4db 100644 --- a/ENUME/projectA/projectA.fdb_latexmk +++ b/ENUME/projectA/projectA.fdb_latexmk @@ -1,5 +1,5 @@ # Fdb version 3 -["pdflatex"] 1636689518 "projectA.tex" "projectA.pdf" "projectA" 1636689519 +["pdflatex"] 1636691043 "projectA.tex" "projectA.pdf" "projectA" 1636691044 "/etc/texmf/web2c/texmf.cnf" 1635008344 475 c0e671620eb5563b2130f56340a5fde8 "" "/usr/share/texlive/texmf-dist/fonts/enc/dvips/base/8r.enc" 1165713224 4850 80dc9bab7f31fb78a000ccfed0e27cab "" "/usr/share/texlive/texmf-dist/fonts/map/fontname/texfonts.map" 1577235249 3524 cb3e574dea2d1052e39280babc910dc8 "" @@ -150,13 +150,13 @@ "errorsA.eps" 1636686866 1613352 8726309b9c94c5647644a20db10b9a1f "" "errorsB.eps" 1636686908 1614432 a5ba1015251d54e816c9a68cee4a6856 "" "iterations.eps" 1636684222 67653 3e4ba61ec0de12fb403d5a37cff1a286 "" - "projectA.aux" 1636689519 11132 20eada88a46049d6ee6d233ae73c741e "pdflatex" - "projectA.out" 1636689519 4281 2a0e019b1a585d28cf17a4b51046f2f3 "pdflatex" - "projectA.tex" 1636689517 46594 ab168ba0c7ab373ffc6e69da68a7221a "" - "projectA.toc" 1636689519 6703 dabbf52c638a1d1edbb1af66800929c8 "pdflatex" + "projectA.aux" 1636691044 11277 caaf7da803404592e338d9a317e94cde "pdflatex" + "projectA.out" 1636691044 4368 1c6cbb3ce7f06483447696a06c0d4a4f "pdflatex" + "projectA.tex" 1636691043 47801 537dfc2ebcc72961c79d33220b0035f5 "" + "projectA.toc" 1636691044 6807 31586b2bc1ebe1d5362c6d33d41b6823 "pdflatex" (generated) - "projectA.toc" - "projectA.out" "projectA.pdf" - "projectA.log" "projectA.aux" + "projectA.toc" + "projectA.log" + "projectA.out" diff --git a/ENUME/projectA/projectA.log b/ENUME/projectA/projectA.log index 15ceac86..62bece29 100644 --- a/ENUME/projectA/projectA.log +++ b/ENUME/projectA/projectA.log @@ -1,4 +1,4 @@ -This is pdfTeX, Version 3.14159265-2.6-1.40.21 (TeX Live 2020/Debian) (preloaded format=pdflatex 2021.10.23) 12 NOV 2021 04:58 +This is pdfTeX, Version 3.14159265-2.6-1.40.21 (TeX Live 2020/Debian) (preloaded format=pdflatex 2021.10.23) 12 NOV 2021 05:24 entering extended mode restricted \write18 enabled. file:line:error style messages enabled. @@ -520,7 +520,7 @@ Package epstopdf Info: Source file: (epstopdf) Command: (epstopdf) \includegraphics on input line 336. Package epstopdf Info: Output file is already uptodate. - + File: errorsA-eps-converted-to.pdf Graphic file (type pdf) Package pdftex.def Info: errorsA-eps-converted-to.pdf used on input line 336. @@ -540,7 +540,7 @@ Package epstopdf Info: Source file: (epstopdf) Command: (epstopdf) \includegraphics on input line 340. Package epstopdf Info: Output file is already uptodate. - + File: errorsB-eps-converted-to.pdf Graphic file (type pdf) Package pdftex.def Info: errorsB-eps-converted-to.pdf used on input line 340. @@ -587,7 +587,7 @@ Package epstopdf Info: Source file: (epstopdf) Command: (epstopdf) \includegraphics on input line 934. Package epstopdf Info: Output file is already uptodate. - + File: iterations-eps-converted-to.pdf Graphic file (type pdf) Package pdftex.def Info: iterations-eps-converted-to.pdf used on input line 934. @@ -598,111 +598,111 @@ Chapter 4. ] [35] Chapter 5. -LaTeX Font Info: Trying to load font information for TS1+fvm on input line 964. +LaTeX Font Info: Trying to load font information for TS1+fvm on input line 989. (/usr/share/texlive/texmf-dist/tex/latex/bera/ts1fvm.fd File: ts1fvm.fd 2004/09/07 scalable font definitions for TS1/fvm. ) LaTeX Font Info: Font shape `TS1/fvm/m/n' will be -(Font) scaled to size 10.20007pt on input line 964. +(Font) scaled to size 10.20007pt on input line 989. [36 ] [37] [38] [39] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1069. +(textcomp) Default family used instead on input line 1094. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1071. +(textcomp) Default family used instead on input line 1096. [40] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1082. +(textcomp) Default family used instead on input line 1107. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1083. +(textcomp) Default family used instead on input line 1108. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1092. +(textcomp) Default family used instead on input line 1117. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1104. +(textcomp) Default family used instead on input line 1129. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1115. +(textcomp) Default family used instead on input line 1140. [41] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1118. +(textcomp) Default family used instead on input line 1143. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1119. +(textcomp) Default family used instead on input line 1144. [42] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1150. -Package textcomp Info: Symbol \textminus not provided by -(textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1151. -Package textcomp Info: Symbol \textminus not provided by -(textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1165. +(textcomp) Default family used instead on input line 1175. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. (textcomp) Default family used instead on input line 1176. +Package textcomp Info: Symbol \textminus not provided by +(textcomp) font family fvm in TS1 encoding. +(textcomp) Default family used instead on input line 1190. +Package textcomp Info: Symbol \textminus not provided by +(textcomp) font family fvm in TS1 encoding. +(textcomp) Default family used instead on input line 1201. [43] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1198. +(textcomp) Default family used instead on input line 1223. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1209. +(textcomp) Default family used instead on input line 1234. [44] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1231. +(textcomp) Default family used instead on input line 1256. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1232. +(textcomp) Default family used instead on input line 1257. [45] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1245. +(textcomp) Default family used instead on input line 1270. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1247. +(textcomp) Default family used instead on input line 1272. [46] Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1282. +(textcomp) Default family used instead on input line 1307. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1284. +(textcomp) Default family used instead on input line 1309. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1286. +(textcomp) Default family used instead on input line 1311. Package textcomp Info: Symbol \textminus not provided by (textcomp) font family fvm in TS1 encoding. -(textcomp) Default family used instead on input line 1314. +(textcomp) Default family used instead on input line 1339. [47] [48] [49] [50 ] (./projectA.aux) Package rerunfilecheck Info: File `projectA.out' has not changed. -(rerunfilecheck) Checksum: 2A0E019B1A585D28CF17A4B51046F2F3;4281. +(rerunfilecheck) Checksum: 1C6CBB3CE7F06483447696A06C0D4A4F;4368. ) Here is how much of TeX's memory you used: - 12894 strings out of 479304 - 223409 string characters out of 5869778 - 882484 words of memory out of 5000000 - 29410 multiletter control sequences out of 15000+600000 + 12897 strings out of 479304 + 223457 string characters out of 5869778 + 882696 words of memory out of 5000000 + 29411 multiletter control sequences out of 15000+600000 423928 words of font info for 77 fonts, out of 8000000 for 9000 1141 hyphenation exceptions out of 8191 - 81i,8n,88p,715b,2253s stack positions out of 5000i,500n,10000p,200000b,80000s + 81i,17n,88p,715b,2253s stack positions out of 5000i,500n,10000p,200000b,80000s {/usr/share/texmf/fonts/enc/dvips/cm-super/cm-super-t1.enc}{/usr/share/texmf/fonts/enc/dvips/cm-super/cm-super-ts1.enc}{/usr/share/texlive/texmf-dist/fonts/enc/dvips/base/8r.enc} -Output written on projectA.pdf (51 pages, 559880 bytes). +Output written on projectA.pdf (51 pages, 562613 bytes). PDF statistics: - 1081 PDF objects out of 1200 (max. 8388607) - 986 compressed objects within 10 object streams - 456 named destinations out of 1000 (max. 500000) - 504 words of extra memory for PDF output out of 10000 (max. 10000000) + 1088 PDF objects out of 1200 (max. 8388607) + 993 compressed objects within 10 object streams + 458 named destinations out of 1000 (max. 500000) + 512 words of extra memory for PDF output out of 10000 (max. 10000000) diff --git a/ENUME/projectA/projectA.out b/ENUME/projectA/projectA.out index da6eae3c..b518c83c 100644 --- a/ENUME/projectA/projectA.out +++ b/ENUME/projectA/projectA.out @@ -29,33 +29,34 @@ \BOOKMARK [1][-]{section.4.1}{Problem}{chapter.4}% 29 \BOOKMARK [1][-]{section.4.2}{Theoretical introduction}{chapter.4}% 30 \BOOKMARK [2][-]{subsection.4.2.1}{Eigenvalues}{section.4.2}% 31 -\BOOKMARK [1][-]{section.4.3}{Solution}{chapter.4}% 32 -\BOOKMARK [1][-]{section.4.4}{Discussion of the result}{chapter.4}% 33 -\BOOKMARK [0][-]{chapter.5}{Code appendix}{}% 34 -\BOOKMARK [1][-]{section.5.1}{Task 2 Code}{chapter.5}% 35 -\BOOKMARK [2][-]{subsection.5.1.1}{Main function}{section.5.1}% 36 -\BOOKMARK [2][-]{subsection.5.1.2}{checkIfMatrixIsSquareMatrix}{section.5.1}% 37 -\BOOKMARK [2][-]{subsection.5.1.3}{gaussianEliminationWithPartialPivoting}{section.5.1}% 38 -\BOOKMARK [2][-]{subsection.5.1.4}{partialPivoting}{section.5.1}% 39 -\BOOKMARK [2][-]{subsection.5.1.5}{partialPivotingSwapOneRow}{section.5.1}% 40 -\BOOKMARK [2][-]{subsection.5.1.6}{swapRowMatrix}{section.5.1}% 41 -\BOOKMARK [2][-]{subsection.5.1.7}{swapValueVector}{section.5.1}% 42 -\BOOKMARK [2][-]{subsection.5.1.8}{gaussianElimination}{section.5.1}% 43 -\BOOKMARK [2][-]{subsection.5.1.9}{substractRows}{section.5.1}% 44 -\BOOKMARK [2][-]{subsection.5.1.10}{backSubstitutionPhase}{section.5.1}% 45 -\BOOKMARK [2][-]{subsection.5.1.11}{iterativeResidualCorrection}{section.5.1}% 46 -\BOOKMARK [2][-]{subsection.5.1.12}{improveSolution}{section.5.1}% 47 -\BOOKMARK [1][-]{section.5.2}{Task 3 code}{chapter.5}% 48 -\BOOKMARK [2][-]{subsection.5.2.1}{initializeValues}{section.5.2}% 49 -\BOOKMARK [2][-]{subsection.5.2.2}{decomposeMatrix}{section.5.2}% 50 -\BOOKMARK [2][-]{subsection.5.2.3}{jacobiLoop}{section.5.2}% 51 -\BOOKMARK [2][-]{subsection.5.2.4}{jacobiInsideLoop}{section.5.2}% 52 -\BOOKMARK [2][-]{subsection.5.2.5}{jacobiEquation}{section.5.2}% 53 -\BOOKMARK [2][-]{subsection.5.2.6}{gaussSeidelLoop}{section.5.2}% 54 -\BOOKMARK [2][-]{subsection.5.2.7}{gaussiInsideLoop}{section.5.2}% 55 -\BOOKMARK [2][-]{subsection.5.2.8}{gaussSeidelEquation}{section.5.2}% 56 -\BOOKMARK [2][-]{subsection.5.2.9}{checkError}{section.5.2}% 57 -\BOOKMARK [2][-]{subsection.5.2.10}{endOfLoop}{section.5.2}% 58 -\BOOKMARK [2][-]{subsection.5.2.11}{dispFinalResults}{section.5.2}% 59 -\BOOKMARK [2][-]{subsection.5.2.12}{plotIterations}{section.5.2}% 60 -\BOOKMARK [2][-]{subsection.5.2.13}{plotIterations}{section.5.2}% 61 +\BOOKMARK [2][-]{subsection.4.2.2}{QR method for finding eigenvalues}{section.4.2}% 32 +\BOOKMARK [1][-]{section.4.3}{Solution}{chapter.4}% 33 +\BOOKMARK [1][-]{section.4.4}{Discussion of the result}{chapter.4}% 34 +\BOOKMARK [0][-]{chapter.5}{Code appendix}{}% 35 +\BOOKMARK [1][-]{section.5.1}{Task 2 Code}{chapter.5}% 36 +\BOOKMARK [2][-]{subsection.5.1.1}{Main function}{section.5.1}% 37 +\BOOKMARK [2][-]{subsection.5.1.2}{checkIfMatrixIsSquareMatrix}{section.5.1}% 38 +\BOOKMARK [2][-]{subsection.5.1.3}{gaussianEliminationWithPartialPivoting}{section.5.1}% 39 +\BOOKMARK [2][-]{subsection.5.1.4}{partialPivoting}{section.5.1}% 40 +\BOOKMARK [2][-]{subsection.5.1.5}{partialPivotingSwapOneRow}{section.5.1}% 41 +\BOOKMARK [2][-]{subsection.5.1.6}{swapRowMatrix}{section.5.1}% 42 +\BOOKMARK [2][-]{subsection.5.1.7}{swapValueVector}{section.5.1}% 43 +\BOOKMARK [2][-]{subsection.5.1.8}{gaussianElimination}{section.5.1}% 44 +\BOOKMARK [2][-]{subsection.5.1.9}{substractRows}{section.5.1}% 45 +\BOOKMARK [2][-]{subsection.5.1.10}{backSubstitutionPhase}{section.5.1}% 46 +\BOOKMARK [2][-]{subsection.5.1.11}{iterativeResidualCorrection}{section.5.1}% 47 +\BOOKMARK [2][-]{subsection.5.1.12}{improveSolution}{section.5.1}% 48 +\BOOKMARK [1][-]{section.5.2}{Task 3 code}{chapter.5}% 49 +\BOOKMARK [2][-]{subsection.5.2.1}{initializeValues}{section.5.2}% 50 +\BOOKMARK [2][-]{subsection.5.2.2}{decomposeMatrix}{section.5.2}% 51 +\BOOKMARK [2][-]{subsection.5.2.3}{jacobiLoop}{section.5.2}% 52 +\BOOKMARK [2][-]{subsection.5.2.4}{jacobiInsideLoop}{section.5.2}% 53 +\BOOKMARK [2][-]{subsection.5.2.5}{jacobiEquation}{section.5.2}% 54 +\BOOKMARK [2][-]{subsection.5.2.6}{gaussSeidelLoop}{section.5.2}% 55 +\BOOKMARK [2][-]{subsection.5.2.7}{gaussiInsideLoop}{section.5.2}% 56 +\BOOKMARK [2][-]{subsection.5.2.8}{gaussSeidelEquation}{section.5.2}% 57 +\BOOKMARK [2][-]{subsection.5.2.9}{checkError}{section.5.2}% 58 +\BOOKMARK [2][-]{subsection.5.2.10}{endOfLoop}{section.5.2}% 59 +\BOOKMARK [2][-]{subsection.5.2.11}{dispFinalResults}{section.5.2}% 60 +\BOOKMARK [2][-]{subsection.5.2.12}{plotIterations}{section.5.2}% 61 +\BOOKMARK [2][-]{subsection.5.2.13}{plotIterations}{section.5.2}% 62 diff --git a/ENUME/projectA/projectA.pdf b/ENUME/projectA/projectA.pdf index 0d5932e1..506d95f7 100644 Binary files a/ENUME/projectA/projectA.pdf and b/ENUME/projectA/projectA.pdf differ diff --git a/ENUME/projectA/projectA.synctex.gz b/ENUME/projectA/projectA.synctex.gz index 8c4fd025..3142ee76 100644 Binary files a/ENUME/projectA/projectA.synctex.gz and b/ENUME/projectA/projectA.synctex.gz differ diff --git a/ENUME/projectA/projectA.tex b/ENUME/projectA/projectA.tex index ac773df0..120c9c17 100755 --- a/ENUME/projectA/projectA.tex +++ b/ENUME/projectA/projectA.tex @@ -947,6 +947,31 @@ This can be rewritten as: \[ (\mathbf{A} - \lambda\mathbf{I})\mathbf{v} = 0 \] And further we get characteristic equation: \[ det(\mathbf{A} - \lambda\mathbf{I}) = 0 \] + +\subsection{QR method for finding eigenvalues} +First we start by transforming matrix \textbf{A} into tridiagonal form. It increases the effectiveness of calculations. + +Single step of QR method consists of: +\[ \textbf{A}^{(k)} = \textbf{Q}^{(k)}\textbf{R}^{(k)} \] +\[ \textbf{A}^{(k+1)} = \textbf{R}^{(k)}\textbf{Q}^{(k)} \] +$\textbf{Q}^{(k)}$ being orthogonal, therefore: +\[ \textbf{R}^{(k)} = (\textbf{Q}^{(k)}))^{-1}\textbf{A}^{(k)} = \textbf{Q}^{(k)T}\textbf{A}^{(k)}\textbf{Q}^{(k)} \] +therefore: +\[ \textbf{A}^{(k+1)} = \textbf{Q}^{(k)T}\textbf{A}^{(k)}\textbf{Q}^{(k)} \] +If we have symmetric matrix \textbf{A} it converges to the diagonal matrix $diag(\lambda_i)$ + +Since this method can be slowly convergent we use shifts. Single iteration of QR method with shifts looks like this: +\begin{equation} +\begin{split} +\textbf{A}^{(k)} - p_k\textbf{I} &= \textbf{Q}^{(k)}\textbf{R}^{(k)} \\ +\textbf{A}^{(k+1)} &= \textbf{R}^{(k)}\textbf{Q}^{(k)} + p_k\textbf{I} \\ +&= \textbf{Q}^{(k)T}(\textbf{A}^{(k)} - p_k\textbf{I})\textbf{Q}^{(k)} + p_k\textbf{I} \\ +&=\textbf{Q}^{(k)T}\textbf{A}^{(k)}\textbf{Q}^{(k)} \\ +\end{split} +\end{equation} + +Wher $p_k$ sholud be chosen as a best estimate of $\lambda_{i+1}$ + \section{Solution} \section{Discussion of the result} diff --git a/ENUME/projectA/projectA.toc b/ENUME/projectA/projectA.toc index 4f0cb514..d70a850e 100644 --- a/ENUME/projectA/projectA.toc +++ b/ENUME/projectA/projectA.toc @@ -50,6 +50,7 @@ \contentsline {section}{\numberline {4.1}Problem}{34}{section.4.1}% \contentsline {section}{\numberline {4.2}Theoretical introduction}{34}{section.4.2}% \contentsline {subsection}{\numberline {4.2.1}Eigenvalues}{34}{subsection.4.2.1}% +\contentsline {subsection}{\numberline {4.2.2}QR method for finding eigenvalues}{34}{subsection.4.2.2}% \contentsline {section}{\numberline {4.3}Solution}{35}{section.4.3}% \contentsline {section}{\numberline {4.4}Discussion of the result}{35}{section.4.4}% \contentsline {chapter}{\numberline {5}Code appendix}{36}{chapter.5}%