From a7d9a9da180526ef2a67c33310b6a875699472ef Mon Sep 17 00:00:00 2001 From: PolishPigeon Date: Thu, 11 Nov 2021 19:54:02 +0100 Subject: [PATCH] Task 2 more functions --- ENUME/projectA/jacobiMethod.asv | 74 ++++++++++++++++++++++----------- ENUME/projectA/jacobiMethod.m | 70 +++++++++++++++++++------------ 2 files changed, 93 insertions(+), 51 deletions(-) diff --git a/ENUME/projectA/jacobiMethod.asv b/ENUME/projectA/jacobiMethod.asv index 1077178a..5bf901a1 100644 --- a/ENUME/projectA/jacobiMethod.asv +++ b/ENUME/projectA/jacobiMethod.asv @@ -1,43 +1,69 @@ function x = jacobiMethod(Matrix, Vector) - [Rows,~] = size(Matrix); + [L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance] = initializeValues(Matrix); + [x, whichIterationAreWeOn, demandedTolerance] = jacobiLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance, Vector); + dispFinalResults(demandedTolerance, whichIterationAreWeOn, Matrix, Vector); +end + +function [L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance] = initializeValues(Matrix) + [Rows, ~] = size(Matrix); + [L, D, U] = decomposeMatrix(Matrix); + initial_x = ones(Rows, 1); + whichIterationAreWeOn = 0; + currentError = inf; % We set it to inf so that the algorithm will always start + % (See condition below) + demandedTolerance = 1e-10; +end + +function [L, D, U] = decomposeMatrix(Matrix) D = diag(diag(Matrix)); - inverseD = inv(D); U = triu(Matrix, 1); % Generates upper triangular part of matrix % where the second variable denotes on which diagonal of matrix should we % start L = tril(Matrix, -1); % Generates lower triangular part of matrix % where the second variable denotes on which diagonal of matrix should we % start - initial_x = ones(Rows, 1); - initial_x = - inverseD * ( ( L + U ) * initial_x) + inverseD * initial_x; % As per formula - % We will be using D \ initial_x and D \ () since it is faster and more - % accurate according to matlab - whichIterationAreWeOn = 0; - currentError = inf; % We set it to inf so that it the algorithm will always start - % (See condition below) - demandedTolerance = 1e-10; +end + +function [x, whichIterationAreWeOn, demandedTolerance] = jacobiLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance, Vector) while currentError >= demandedTolerance - - x = - inverseD * ( ( L + U ) * initial_x) + inverseD * initial_x; % As per formula - - currentError = norm(x - initial_x); - disp(currentError); - if currentError <= demandedTolerance - currentError = norm(Matrix*x-Vector); - disp(currentError); - if currentError <= demandedTolerance - break; - else - demandedTolerance = demandedTolerance * 2; - end + x = jacobiEquation(D, L, U, initial_x, Vector); + [flag, demandedTolerance] = checkError(x, initial_x, demandedTolerance, Matrix, Vector); + if flag == 1 + break end initial_x = x; whichIterationAreWeOn = whichIterationAreWeOn + 1; end +end + +function x = jacobiEquation(D, L, U, initial_x, Vector) + x = - D \ ( L + U ) * initial_x + D \ Vector; % As per formula + % We will be using D \ Vector and D \ ( ) instead of inverseD since + % this is faster according to matlab +end + +function [flag, demandedTolerance] = checkError(x, initial_x, demandedTolerance, Matrix, Vector) + flag = 0; + currentError = norm(x - initial_x); + if currentError <= demandedTolerance + currentError = norm(Matrix*x-Vector); + if currentError <= demandedTolerance + flag = 1; + else + demandedTolerance = demandedTolerance * 2; + end + end +end + +function [initial_x, ] + + + +function dispFinalResults(demandedTolerance, whichIterationAreWeOn, Matrix, Vector) disp("Final demandedTolerance"); disp(demandedTolerance); disp("Final Iteration: "); disp(whichIterationAreWeOn); disp("A\b matlab:"); - disp(Matrix / Vector); + disp(Matrix \ Vector); end \ No newline at end of file diff --git a/ENUME/projectA/jacobiMethod.m b/ENUME/projectA/jacobiMethod.m index be8f1848..3e0bfde8 100644 --- a/ENUME/projectA/jacobiMethod.m +++ b/ENUME/projectA/jacobiMethod.m @@ -1,22 +1,16 @@ function x = jacobiMethod(Matrix, Vector) - [L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance] = initializeValues(Matrix) - x = jacobiLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance, Vector); - disp("Final demandedTolerance"); - disp(demandedTolerance); - disp("Final Iteration: "); - disp(whichIterationAreWeOn); - disp("A\b matlab:"); - disp(Matrix \ Vector); + [L, D, U, initial_x, whichIterationAreWeOn, demandedTolerance, flag] = initializeValues(Matrix); + [x, whichIterationAreWeOn, demandedTolerance] = jacobiLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, demandedTolerance, Vector, flag); + dispFinalResults(demandedTolerance, whichIterationAreWeOn, Matrix, Vector); end -function [L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance] = initializeValues(Matrix) +function [L, D, U, initial_x, whichIterationAreWeOn, demandedTolerance, flag] = initializeValues(Matrix) [Rows, ~] = size(Matrix); [L, D, U] = decomposeMatrix(Matrix); initial_x = ones(Rows, 1); whichIterationAreWeOn = 0; - currentError = inf; % We set it to inf so that the algorithm will always start - % (See condition below) demandedTolerance = 1e-10; + flag = 0; end function [L, D, U] = decomposeMatrix(Matrix) @@ -29,28 +23,50 @@ function [L, D, U] = decomposeMatrix(Matrix) % start end -function x = jacobiLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, currentError, demandedTolerance, Vector) - while currentError >= demandedTolerance - x = jacobiEquation(D, L, U, initial_x, Vector); - currentError = norm(x - initial_x); - %disp(currentError); - if currentError <= demandedTolerance - currentError = norm(Matrix*x-Vector); - %disp(currentError); - if currentError <= demandedTolerance - break; - else - demandedTolerance = demandedTolerance * 2; - end - end - initial_x = x; - whichIterationAreWeOn = whichIterationAreWeOn + 1; +function [x, whichIterationAreWeOn, demandedTolerance] = jacobiLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, demandedTolerance, Vector, flag) + while flag ~= 1 + [x, whichIterationAreWeOn, demandedTolerance, flag, initial_x] = jacobiInsideLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, demandedTolerance, Vector); end end +function [x, whichIterationAreWeOn, demandedTolerance, flag, initial_x] = jacobiInsideLoop(Matrix, L, D, U, initial_x, whichIterationAreWeOn, demandedTolerance, Vector) + x = jacobiEquation(D, L, U, initial_x, Vector); + [flag, demandedTolerance] = checkError(x, initial_x, demandedTolerance, Matrix, Vector); + [initial_x, whichIterationAreWeOn] = endOfLoop(x, whichIterationAreWeOn); +end function x = jacobiEquation(D, L, U, initial_x, Vector) x = - D \ ( L + U ) * initial_x + D \ Vector; % As per formula % We will be using D \ Vector and D \ ( ) instead of inverseD since % this is faster according to matlab +end + +function [flag, demandedTolerance] = checkError(x, initial_x, demandedTolerance, Matrix, Vector) + flag = 0; + currentError = norm(x - initial_x); + if currentError <= demandedTolerance + currentError = norm(Matrix*x-Vector); + if currentError <= demandedTolerance + flag = 1; + else + demandedTolerance = demandedTolerance * 2; + end + end +end + +function [initial_x, whichIterationAreWeOn, flag] = endOfLoop(x, whichIterationAreWeOn) + initial_x = x; + whichIterationAreWeOn = whichIterationAreWeOn + 1; + flag = 0; +end + + + +function dispFinalResults(demandedTolerance, whichIterationAreWeOn, Matrix, Vector) + disp("Final demandedTolerance"); + disp(demandedTolerance); + disp("Final Iteration: "); + disp(whichIterationAreWeOn); + disp("A\b matlab:"); + disp(Matrix \ Vector); end \ No newline at end of file