Adding references

ENUME/references/ENUME_Project1/README.md

@@ -0,0 +1,30 @@
+# ENUME_Project1
+A Numerical Method project cointaining 4 tasks.
+## 1. 
+Write a program finding macheps in the MATLAB environment on a lab computer or your
+computer. 
+## 2. 
+Write a general program solving a system of n linear equations Ax = b using the indicated
+method Apply the program to solve the system of linear equations for given matrix A and vector b, for increasing numbers
+of equations n = 10,20,40,80,160,… until the solution time becomes prohibitive (or the
+method fails), for:
+### a) 
+Aij { for i = j -> 13, for i = j-1 or i = j+1 -> 4, other -> 0 \
+Bi = 2.4 + 0.6i, \
+i,j = 1..n;
+### b)
+Aij = 4/[5(i + j – 1)], Bi = 1/(2 i), i – odd; Bi = 0, i – even, i, j = 1,…,n;
+
+## 3
+Write a general program for solving the system of n linear equations Ax = b using the
+Gauss-Seidel and Jacobi iterative algorithms. Apply it for the system: \
+13x1 + 2x2 – 8x3 + x4 =16 \
+x1 + 10x2 + 5x3 – 2x4 = 24 \
+6x1 + 2x2 – 23x3 + 15x4 = 184 \
+x1 + 2x2 – x3 + 13x4 = 82
+
+## 4
+Write a program of the QR method for finding eigenvalues of 5×5 matrices: \
+a) without shifts; \
+b) with shifts calculated on the basis of an eigenvalue of the 2×2 right-lower-corner 
+submatrix. 

ENUME/references/ENUME_Project1/ResidualCorrection.m

@@ -0,0 +1,30 @@
+function [TableOfErrors] = ResidualCorrection(n, subpoint)
+    TableOfErrors = zeros(n, 1);
+    TableOfErrorsAC = zeros(n, 1);
+    TableOfErrorsAC2 = zeros(n, 1);
+    TableOfNValues = zeros(n, 1);
+    for i = 1:n
+        TableOfNValues(i) = 10*2^(i-1);
+        disp(i);
+        x = TASK2(10*2^(i-1));
+        if subpoint == 'A'
+            x = TaskA(x);
+        elseif subpoint == 'B'
+            x = TaskB(x);
+        end
+        TableOfErrors(i) = norm(x.errors);
+        x = ResidualCorrection(x);
+        x= GetErrors(x);
+        TableOfErrorsAC(i) = norm(x.errors);
+        x = ResidualCorrection(x);
+        x= GetErrors(x);
+        TableOfErrorsAC2(i) = norm(x.errors);
+    end
+    plot(TableOfNValues, TableOfErrors, '-o');
+    hold on;
+    ylabel('Maximal absolute value of an error') 
+    xlabel('N')
+    plot(TableOfNValues, TableOfErrorsAC, '-or');
+    plot(TableOfNValues, TableOfErrorsAC2, '-og');
+    hold off;
+end

ENUME/references/ENUME_Project1/TASK1.m

@@ -0,0 +1,8 @@
+function macheps = Task1()
+    macheps = 1;
+    while 1 + (macheps/2)>1
+        macheps = macheps/2;
+    end
+end
+
+

ENUME/references/ENUME_Project1/TASK2.m

@@ -0,0 +1,161 @@
+classdef TASK2
+    properties
+        %main array size n by n
+        A;
+        %second array size n
+        B;
+        %answer
+        x;
+        %errors
+        errors;
+        %size of a matricx
+        n;
+    end
+    methods (Access = 'protected')
+        function [obj, locx] = gaussWithPartialPiv(obj)
+            locx = zeros(obj.n, 1);
+            for j = 1:obj.n
+                %finding the maximal value of a column
+                index = j;
+                locMax = abs(obj.A(j, j));
+                for i = j+1:obj.n
+                    if abs(obj.A(i, j)) > locMax
+                        locMax = abs(obj.A(i, j));
+                        index = i;
+                    end
+                end
+                %checking neccessary condition for the further execution
+                if locMax == 0
+                    disp(locMax);
+                    disp('Wrong Matrix');
+                    return
+                end
+                if index ~= j
+                    %swaping the rows in neccessary (highest value found
+                    %not in j row)
+                    obj.A([index, j],:) = obj.A([j, index],:);
+                    obj.B([index, j]) = obj.B([j, index]);
+                end
+                for i = j+1:obj.n
+                    r = obj.A(i, j)/obj.A(j,j);
+                    %if r == 0 we dont need to continue
+                    if r ~= 0
+                        %changing the main array A
+                        for locj = j+1:obj.n
+                            obj.A(i, locj) = obj.A(i, locj) - r*obj.A(j, locj);
+                        end
+                        %and the B array aswell
+                        obj.B(i) = obj.B(i) - r*obj.B(j);
+                    end
+                end
+            end
+            %now finally we can obtain the results
+            %Xn value is rather obvious
+            locx(obj.n) = obj.B(obj.n)/obj.A(obj.n, obj.n);
+            for i = obj.n-1:-1:1
+                buffor = 0;
+                for j = i+1:obj.n
+                    buffor = buffor + obj.A(i, j)*locx(j);
+                end
+                locx(i) = (obj.B(i) - buffor)/obj.A(i,i);
+            end
+        end
+        function obj = TaskAArray(obj)
+            [row, columns] = size(obj.A);
+            if row ~= size(obj.B)
+                disp('A and B array size is different!');
+                return 
+            end
+            for i = 1:row
+                %setting the A array
+                for j = 1:columns
+                    if i == j
+                        obj.A(i, j) = 13; 
+                    elseif i == j-1 || i == j+1
+                        obj.A(i, j) = 4;
+                    else
+                        obj.A(i, j) = 0;
+                    end
+                end
+                %setting he B array
+                obj.B(i) = 2.4 + 0.6*i;
+            end
+        end
+        function obj = TaskBArray(obj)
+            [row, columns] = size(obj.A);
+            if row ~= size(obj.B)
+                disp('A and B array size is differn!');
+                return 
+            end
+            for i = 1:row
+                %setting the A array
+                for j = 1:columns
+                    obj.A(i, j) = 4/(5*(i + j - 1));
+                end
+                %setting he B array
+                if mod(i, 2) == 0
+                    obj.B(i) = 1/(2*i);
+                else
+                    obj.B(i) = 0;
+                end
+            end
+        end
+        function obj = SetSize(obj, n)        
+            obj.A = zeros(n);
+            obj.B = zeros(n, 1);
+            obj.x = zeros(n,1);
+            obj.errors = zeros(n,1);
+            obj.n = n;
+        end
+    end
+    %public methods
+    methods
+        function obj = TASK2(n)
+            obj.n = n;
+        end
+        function obj = GetErrors(obj)  
+            %colculating the error following the formula r = Ax - b
+            for i = 1:obj.n
+                result = 0;
+                for j = i:obj.n
+                    result = result + obj.A(i,j) * obj.x(j);
+                end
+                obj.errors(i) = result - obj.B(i);
+            end
+        end
+        function obj = ResidualCorrection(obj)
+            newObj = TASK2(obj.n);
+            newObj = obj;
+            newObj.B = newObj.errors;
+            [newObj, newObj.x] = gaussWithPartialPiv(newObj);
+            for i = 1:obj.n
+                obj.x(i) = obj.x(i) - newObj.x(i);
+            end
+        end
+        function obj = TaskA(obj)
+            obj = SetSize(obj, obj.n);
+            obj = TaskAArray(obj);
+            [obj, obj.x] = gaussWithPartialPiv(obj);
+            obj = GetErrors(obj);
+        end
+        function obj = TaskB(obj)
+            obj = SetSize(obj, obj.n);
+            obj = TaskBArray(obj);
+            [obj, obj.x] = gaussWithPartialPiv(obj);
+            obj = GetErrors(obj);
+        end
+        function DispSolutionAndError(obj)
+            format long
+            disp('x for n ='); disp(obj.n);
+            obj.x
+            disp('errors for n ='); disp(obj.n);
+            obj.errors
+        end
+         function DispObjects(obj)
+            obj.A
+            obj.B
+            obj.x
+            obj.errors
+         end
+    end
+end

ENUME/references/ENUME_Project1/TASK3.m

@@ -0,0 +1,155 @@
+classdef TASK3 < TASK2
+    properties
+        acc = 1e-10;
+    end
+    methods (Access = 'private')
+        function [obj, array]= gaussSeidel(obj, type)
+            array = [];
+            err = inf;
+            k = 1;
+            m = 10000;
+            while  k <= m && err > obj.acc
+                localMax = 0;
+                for i = 1 : obj.n 
+                    s = 0;
+                    for j = 1 : obj.n 
+                        s = s+(-obj.A(i,j)*obj.x(j));
+                    end
+                    s = (s+obj.B(i))/obj.A(i,i);
+                    if abs(s) > localMax 
+                        localMax  = abs(s);
+                    end
+                    obj.x(i) = obj.x(i) + s;
+                end
+                if err > localMax
+                    err = localMax ;
+                end
+                k = k+1;
+                obj = GetErrors(obj);
+                if(type == 'N')%norm
+                    array(end+1) = norm(obj.errors);
+                elseif(type == 'E')%maximal value of s
+                    array(end+1) = localMax;
+                end
+            end
+            disp(k - 1);
+        end
+        function [obj, array]= jacobi(obj, type)
+            array = [];
+            err = inf;
+            k = 1;
+            m = 10000;
+            x1 = zeros(obj.n, 1);
+            while  k <= m && err > obj.acc
+                localMax = 0;
+                for i = 1 : obj.n 
+                    s = 0;
+                    for j = 1 : obj.n 
+                        s = s+(-obj.A(i,j)*x1(j));
+                    end
+                    x1 = obj.x;
+                    s = (s+obj.B(i))/obj.A(i,i);
+                    obj.x(i) = obj.x(i) + s;
+                    if abs(s) > localMax 
+                        localMax  = abs(s);
+                    end
+                end
+                if err > localMax
+                    err = localMax ;
+                end
+                k = k+1;
+                obj = GetErrors(obj);
+                if(type == 'N')%norm
+                    array(end+1) = norm(obj.errors);
+                elseif(type == 'E')%maximal value of s
+                    array(end+1) = localMax;
+                end
+            end
+            disp(k - 1);
+        end
+	end
+    methods
+        function obj = TASK3()
+            obj = obj@TASK2(4);
+        end
+        function obj = SetTask3(obj)
+            obj = SetSize(obj, 4);
+            obj.A = [13 2 -8 1; 1 10 5 -2; 6 2 -23 15; 1 22 -1 13];
+            obj.B = [16 24 184 82];
+        end
+        function obj = Example(obj)
+            obj = SetSize(obj, 4);
+            obj.A = [5 -2 3 0; -3 9 1 -2; 2 -1 -7 1; 4 3 -5 7];
+            obj.B = [-1 2 3 0.5];
+            [obj, array] = gaussSeidel(obj, 'N');
+            plot(array, '-o');
+        end
+        %type in 'G' for Gauss Seidel and 'J' for Jacobi method.
+        function obj = Task3System(obj, alg)
+            obj = SetSize(obj, 4);
+            obj.A = [13 2 -8 1; 1 10 5 -2; 6 2 -23 15; 1 22 -1 13];
+            obj.B = [16 24 184 82];
+            if alg == 'G'
+                [obj, array] = gaussSeidel(obj);
+                plot(array, '-o');
+            elseif alg == 'J'
+                [obj, array] = jacobi(obj);
+                plot(array, '-or');
+            end
+        end
+        function obj = Task3a(obj, type)
+            if type ~= 'N' && type ~= 'E'
+                disp('Wrong type parameter!');
+                return
+            end
+            obj = SetTask3(obj);
+            [obj, array1] = gaussSeidel(obj, type);
+            obj = SetTask3(obj);
+            [obj, array2] = jacobi(obj, type);
+            plot(array1, 'ob-');
+            hold on
+            plot(array2, 'or-');
+            hold off
+        end
+        function obj = Task3With2a(obj, input, type)
+            obj = SetSize(obj, 10);
+            obj = TaskAArray(obj);
+            if input == 'G'
+                [obj, table] = gaussSeidel(obj, type);
+                plot(table, '-o');
+            elseif input == 'J'
+                [obj, table] = jacobi(obj,type);
+                plot(table, '-o');
+            elseif input == 'B'
+                [obj, table] = gaussSeidel(obj,type);
+                obj = SetSize(obj, 10);
+                obj = TaskAArray(obj);
+                [obj, table2] = jacobi(obj,type);
+                plot(table, '-o');
+                hold on
+                plot(table2, '-or');
+                hold off
+            end
+        end
+        function obj = Task3With2b(obj, input, type)
+            obj = SetSize(obj, 10);
+            obj = TaskBArray(obj);
+            if input == 'G'
+                [obj, table] = gaussSeidel(obj, type);
+                plot(table, '-o');
+            elseif input == 'J'
+                [obj, table] = jacobi(obj,type);
+                plot(table, '-o');
+            elseif input == 'B'
+                [obj, table] = gaussSeidel(obj,type);
+                obj = SetSize(obj, 10);
+                obj = TaskBArray(obj);
+                [obj, table2] = jacobi(obj,type);
+                plot(table, '-o');
+                hold on
+                plot(table2, '-o');
+                hold off
+            end
+        end
+    end
+end

ENUME/references/ENUME_Project1/TASK4.m

@@ -0,0 +1,95 @@
+classdef TASK4 < TASK2
+    properties
+        tol = 1e-6;
+        imax = 10000;
+    end
+    methods
+        function obj = TASK4()
+            obj = obj@TASK2(5);
+        end
+        function [obj, Q, R]=QR(obj)
+            Q=zeros(obj.n,obj.n);
+            R=zeros(obj.n,obj.n);
+            d=zeros(1,obj.n);
+            for i=1:obj.n
+                Q(:,i)=obj.A(:,i);
+                R(i,i)=1;
+                d(i)=Q(:,i)'*Q(:,i);
+                for j=i+1:obj.n
+                    R(i,j)=(Q(:,i)'*obj.A(:,j))/d(i);
+                    obj.A(:,j)=obj.A(:,j)-R(i,j)*Q(:,i);
+                end
+            end
+            for i=1:obj.n
+                dd=norm(Q(:,i));
+                Q(:,i)=Q(:,i)/dd;
+                R(i,i:obj.n)=R(i,i:obj.n)*dd;
+            end
+        end
+        function [obj, eig, i] = EigvalQRshifts(obj)
+            obj.n=size(obj.A,1);
+            eig=diag(ones(obj.n));
+            INITIALsubmatrix=obj.A; 
+            for k=obj.n:-1:2
+                DK=INITIALsubmatrix;
+                i=0;
+                while i<=obj.imax && max(abs(DK(k,1:k-1)))>obj.tol
+                    DD=DK(k-1:k,k-1:k);
+                    [ev1,ev2]=roots(1,-(DD(1,1)+DD(2,2)),DD(2,2)*DD(1,1)-DD(2,1)*DD(1,2));
+                    if abs(ev1-DD(2,2)) < abs(ev2-DD(2,2))
+                        shift=ev1;
+                    else
+                        shift=ev2;
+                    end
+                    DP=DK-eye(k)*shift;
+                    [Q1,R1]=qr(DP);
+                    DK=R1*Q1+eye(k)*shift;
+                    i=i+1;
+                end
+                if i > obj.imax
+                    error("Too many iterations!");
+                end
+                eig(k)=DK(k,k);
+                if k > 2
+                    INITIALsubmatrix=DK(1:k-1,1:k-1);
+                else
+                    eig(1)=DK(1,1);
+                end
+                DK
+            end
+        end
+        function [obj, eigenvalues, i] = EigvalQRNoShift(obj)
+            i=1;
+            while i <= obj.imax && max(max(obj.A-diag(diag(obj.A)))) > obj.tol
+                [obj, Q1,R1] = QR(obj);
+                obj.A=R1*Q1;
+                i=i+1;
+            end
+            if i > obj.imax
+                error("Too many iterations!");
+            end
+            eigenvalues=diag(obj.A);
+        end
+        function obj = SetExampleShifts(obj)
+            obj = SetSize(obj, 5);
+            obj.A = [2 33 8 -3 4; 33 1 -6 5 -3; 8 -6 -5 -6 8; -3 5 -6 3 2; 4 -3 8 2 45];
+            [obj, eig, k] = EigvalQRshifts(obj);
+            eig
+            k
+        end
+        function obj = SetExample(obj)
+            obj = SetSize(obj, 5);
+            obj.A = [2 33 8 -3 4; 33 1 -6 5 -3; 8 -6 -5 -6 8; -3 5 -6 3 2; 4 -3 8 2 45];
+            [obj, eig, k] = EigvalQRNoShift(obj);
+            eig
+            k
+        end
+    end
+end
+function [x1, x2] = roots(a, b, c)
+    first = -b + sqrt(b * b - 4 * a * c);
+	second = -b - sqrt(b * b - 4 * a * c);
+	li = max(abs(first), abs(second));
+	x1 = li/(2*a);
+	x2 = ((-b)/a);
+end

ENUME/references/ENUME_ProjectA_No31/task 1/macheps.m

@@ -0,0 +1,8 @@
+function [macheps] = macheps()
+%MACHEPS function calculates machine epsilon
+%   finds smallest number E that 1+E>1
+macheps = 1;
+while( 1 + ( macheps / 2 ) > 1 )
+    macheps = macheps/2;
+end
+end

ENUME/references/ENUME_ProjectA_No31/task 2/euclideanNorm.m

@@ -0,0 +1,21 @@
+function [r] = euclideanNorm(rMatrix)
+%EUCLIDEANNORM calculates Euclidean Norm
+%   yes
+sr = size(rMatrix);
+if( size(sr) > 2 )
+    return;
+end
+if( sr(1) == 1 )
+    n = sr(2);
+elseif( sr(2) == 1 );
+    n = sr(1);
+else
+    return;
+end
+sum = 0;
+for i = 1:n
+    sum = sum + rMatrix(i)*rMatrix(i);
+end
+r = sqrt(sum);
+end
+

ENUME/references/ENUME_ProjectA_No31/task 2/matrixGen2a.m

@@ -0,0 +1,21 @@
+function [A, b] = matrixGen2a(n)
+%matrixGen2a Generates matrices for the task 2a
+%   generates matrices for a system Ax = B of n linear equations
+
+A = zeros(n, n);
+b = zeros(n, 1);
+
+for i = 1:n
+    b(i) = 0.9*i;
+   for j = 1:n
+       if(i == j)
+           A(i, j) = 11;
+       elseif (i == j-1)
+           A(i, j) = 5;
+       elseif (i == j+1)
+           A(i, j) = 5;
+       end
+   end
+end
+end
+

ENUME/references/ENUME_ProjectA_No31/task 2/matrixGen2b.m

@@ -0,0 +1,18 @@
+function [A,b] = matrixGen2b(n)
+%matrixGen2b Generates matrices for the task 2b
+%   generates matrices for a system Ax = B of n linear equations
+
+A = zeros(n, n);
+b = zeros(n, 1);
+
+for i = 1:n
+   if( mod(i,2) == 0 )
+       b(i) = 2/(3*i);
+   else
+       b(i) = 0;
+   end
+   for j = 1:n
+       A(i,j) = 7/(8*(i+j+1));
+   end 
+end
+end

ENUME/references/ENUME_ProjectA_No31/task 2/residualCorrection.m

@@ -0,0 +1,17 @@
+function [Rnew,Xnew] = residualCorrection(r,x,AT,b,A)
+%RESIDUALCORRECTION calculates an iteration of residual correction
+%   currently not working properly
+n = length(x);
+dx = zeros(n, 1);
+for o = 1:n
+    k = n-o+1;
+    sum = 0;
+    for j = k+1:n
+        sum = sum + (AT(k,j) * dx(j));
+    end
+    dx(k) = (r(k) - sum)/AT(k,k);
+end
+Xnew = x - dx;
+Rnew = A*Xnew - b;
+end
+

ENUME/references/ENUME_ProjectA_No31/task 2/solveIndicated.m

@@ -0,0 +1,46 @@
+function [x, A, b] = solveIndicated(A, b)
+%SOLVEINDICATED solves system of linear equations Ax=b
+%   uses the indicated method (Gaussain elimination with partial pivoting)
+%   returns x - the solution and A,b - the system after transformation
+
+sa = size(A);
+sb = size(b);
+if( sa(1) ~= sa(2) )
+    return;
+elseif ( sa(1) ~= sb(1) )
+    return;
+end
+
+n = sa(1);
+
+for k = 1:n
+%   Partial pivoting
+    i = k;
+    for j = k+1:n
+        if ( abs(A(j, k)) > abs(A(i, k)) )
+            i = j;
+        end
+    end
+    if( k~=i )
+        A([k i], :) = A([i k], :);
+        b([k i]) = b([i k]);
+    end
+%   Gauss transform
+    for j = k+1:n
+        l = A(j,k) / A(k,k);
+        A(j, :) = A(j, :) - A(k, :) * l;
+        b(j) = b(j) - b(k) * l;
+    end
+end
+%Backwards substitution
+x = zeros(n, 1);
+for o = 1:n
+    k = n-o+1;
+    sum = 0;
+    for j = k+1:n
+        sum = sum + (A(k,j) * x(j));
+    end
+    x(k) = (b(k) - sum)/A(k,k);
+end
+end
+

ENUME/references/ENUME_ProjectA_No31/task 2/task2.m

@@ -0,0 +1,130 @@
+clear all;
+
+repeatSmall = 4;
+repeatBig = 9;
+repeatEqua = 10;
+%-----TASK 2A
+%big graph
+repeats = zeros(1,repeatBig);
+errors = zeros(1,repeatBig);
+n = repeatEqua;
+for i = 1:1:repeatBig
+    [A, b] = matrixGen2a( n );
+    [x] = solveIndicated( A, b );
+    repeats(i) = n;
+    errors(i) = euclideanNorm( A*x - b );
+    n = n*2;
+end
+figure(1)
+plot(repeats, errors, 'o');
+title(sprintf("Errors in task 2a, up to %d equations", n/2));
+xlabel('Number of equations');
+ylabel('Euclidean norm of errors');
+grid on;
+box off;
+saveas(1, "./plots/2aBIG.fig");
+saveas(1, "./plots/2aBIG.png");
+%small
+repeats = zeros(1,repeatSmall);
+errors = zeros(1,repeatSmall);
+n = repeatEqua;
+for i = 1:1:repeatSmall
+    [A, b] = matrixGen2a( n );
+    [x] = solveIndicated( A, b );
+    repeats(i) = n;
+    errors(i) = euclideanNorm( A*x - b );
+    n = n*2;
+end
+figure(2)
+plot(repeats, errors, 'o');
+title(sprintf("Errors in task 2a, up to %d equations", n/2));
+xlabel('Number of equations');
+ylabel('Euclidean norm of errors');
+grid on;
+box off;
+saveas(2, "./plots/2aSMALL.fig");
+saveas(2, "./plots/2aSMALL.png");
+[Aa, ba] = matrixGen2a( 10 );
+[xa, AaT, baT] = solveIndicated( Aa, ba );
+Ra = Aa*xa - ba;
+%TODO: RESIDUAL
+repeatResidual = 10;
+residualA = zeros(1,repeatResidual+1);
+residualA(1) = euclideanNorm(Ra);
+xaR = xa;
+RaR = Ra;
+for i = 1:1:repeatResidual
+    [RaR, xaR] = residualCorrection(RaR,xaR,AaT,ba,Aa);
+    residualA(i+1) = euclideanNorm(RaR);
+end
+figure(5)
+plot(0:1:repeatResidual,residualA,'o');
+title("Results of residual correction");
+xlabel('Iteration');
+ylabel('Euclidean norm of errors');
+grid on;
+box off;
+saveas(5, "./plots/residualA.fig");
+saveas(5, "./plots/residualA.png");
+%-----TASK 2B
+%big graph
+repeats = zeros(1,repeatBig);
+errors = zeros(1,repeatBig);
+n = repeatEqua;
+for i = 1:1:repeatBig
+    [A, b] = matrixGen2b( n );
+    [x] = solveIndicated( A, b );
+    repeats(i) = n;
+    errors(i) = euclideanNorm( A*x - b );
+    n = n*2;
+end
+figure(3)
+plot(repeats, errors, 'o');
+title(sprintf("Errors in task 2b, up to %d equations", n/2));
+xlabel('Number of equations');
+ylabel('Euclidean norm of errors');
+grid on;
+box off;
+saveas(3, "./plots/2bBIG.fig");
+saveas(3, "./plots/2bBIG.png");
+%small
+repeats = zeros(1,repeatSmall);
+errors = zeros(1,repeatSmall);
+n = repeatEqua;
+for i = 1:1:repeatSmall
+    [A, b] = matrixGen2b( n );
+    [x] = solveIndicated( A, b );
+    repeats(i) = n;
+    errors(i) = euclideanNorm( A*x - b );
+    n = n*2;
+end
+figure(4)
+plot(repeats, errors, 'o');
+title(sprintf("Errors in task 2b, up to %d equations", n/2));
+xlabel('Number of equations');
+ylabel('Euclidean norm of errors');
+grid on;
+box off;
+saveas(4, "./plots/2bSMALL.fig");
+saveas(4, "./plots/2bSMALL.png");
+[Ab, bb] = matrixGen2b( 10 );
+[xb, AbT, bbT] = solveIndicated( Ab, bb );
+Rb = Ab*xb - bb;
+%TODO: RESIDUAL
+residualB = zeros(1,repeatResidual+1);
+residualB(1) = euclideanNorm(Rb);
+xbR = xb;
+RbR = Rb;
+for i = 1:1:repeatResidual
+    [RbR, xbR] = residualCorrection(Rb,xb,AbT,bb,Ab);
+    residualB(i+1) = euclideanNorm(RbR);
+end
+figure(6)
+plot(0:1:repeatResidual,residualB,'o');
+title("Results of residual correction");
+xlabel('Iteration');
+ylabel('Euclidean norm of errors');
+grid on;
+box off;
+saveas(6, "./plots/residualB.fig");
+saveas(6, "./plots/residualB.png");

ENUME/references/ENUME_ProjectA_No31/task 3/GaussSeidelMethod.m

@@ -0,0 +1,40 @@
+function [x, errors] = GaussSeidelMethod(A, b)
+%GAUSSSEIDELMETHOD solves a system Ax = b using the Gauss-Seidel iterative method
+%   The accuracy target is 10e-10
+%   returns x - the solution and errors - vector of errors for each
+%   iteration
+
+[L, D, U] = decomposeLDU(A);
+
+%Checking covergence conditions
+if ~(rowDominant(A) || columnDominant(A))
+    sr = max(abs(eig(-inv(D)*(L+U))));
+    if sr >= 1
+        x = sr;
+        return
+    end
+end
+
+%Initial guess
+x = zeros(length(A), 1);
+
+n = length(A);
+iteration = 1;
+errors(iteration) = vecnorm(A*x - b);
+
+while errors(iteration) > 10^-10
+    w = U*x - b;
+    for i = 1:1:n
+        sum = 0;
+        for j = 1:1:i-1
+            sum = sum - L(i,j) * x(j);
+        end
+        sum = sum - w(i);
+        x(i) = sum / D(i,i);
+    end
+    
+    iteration = iteration + 1;
+    
+    errors(iteration) = vecnorm(A*x - b);
+end
+

ENUME/references/ENUME_ProjectA_No31/task 3/JacobiMethod.m

@@ -0,0 +1,32 @@
+function [x, errors] = JacobiMethod(A, b)
+%JACOBIMETHOD solves a system Ax = b using the Jacobi iterative method
+%   The accuracy target is 10e-10
+%   returns x - the solution and errors - vector of errors for each
+%   iteration
+
+[L, D, U] = decomposeLDU(A);
+
+%Checking covergence conditions
+if ~(rowDominant(A) || columnDominant(A))
+    sr = max(abs(eig(-inv(D)*(L+U))));
+    if sr >= 1
+        x = sr;
+        return
+    end
+end
+
+%Initial guess
+x = zeros(length(A), 1);
+
+Di = inv(D);
+iteration = 1;
+errors(iteration) = vecnorm(A*x - b);
+
+while errors(iteration) > 10^-10
+    x = -Di * (L+U) * x + Di*b;
+    
+    iteration = iteration + 1;
+    
+    errors(iteration) = vecnorm(A*x - b);
+end
+

ENUME/references/ENUME_ProjectA_No31/task 3/columnDominant.m

@@ -0,0 +1,20 @@
+function [out] = columnDominant(A)
+%COLUMNDOMINANT checks column dominance of matrix A
+%   returns true when dominant, false otherwise
+n = size(A);
+
+for j = 1:1:n
+    sum = 0;
+    for i = 1:1:n
+        if i == j
+            continue
+        end
+        sum = sum + abs(A(i,j));
+    end
+    if abs(A(j,j)) <= sum
+        out = false;
+        return
+    end
+end
+out = true;
+end

ENUME/references/ENUME_ProjectA_No31/task 3/decomposeLDU.m

@@ -0,0 +1,27 @@
+function [L,D,U] = decomposeLDU(A)
+%DECOMPOSELDU decomposes the matrix A into L+D+U
+%   returns L - lowerdiagonal, D - diagonal, U - upperdiagonal
+
+n = size(A);
+L = zeros(n);
+D = zeros(n);
+U = zeros(n);
+
+for i = 2:1:n
+    for j = 1:1:i-1
+        L(i,j) = A(i,j);
+    end
+end
+
+for i = 1:1:n
+    D(i,i) = A(i,i);
+end
+
+for i = 1:1:n
+    for j = i+1:1:n
+        U(i,j) = A(i,j);
+    end
+end
+
+end
+

ENUME/references/ENUME_ProjectA_No31/task 3/matrixGen2a.m

@@ -0,0 +1,21 @@
+function [A, b] = matrixGen2a(n)
+%matrixGen2a Generates matrices for the task 2a
+%   generates matrices for a system Ax = B of n linear equations
+
+A = zeros(n, n);
+b = zeros(n, 1);
+
+for i = 1:n
+    b(i) = 0.9*i;
+   for j = 1:n
+       if(i == j)
+           A(i, j) = 11;
+       elseif (i == j-1)
+           A(i, j) = 5;
+       elseif (i == j+1)
+           A(i, j) = 5;
+       end
+   end
+end
+end
+

ENUME/references/ENUME_ProjectA_No31/task 3/matrixGen2b.m

@@ -0,0 +1,18 @@
+function [A,b] = matrixGen2b(n)
+%matrixGen2b Generates matrices for the task 2b
+%   generates matrices for a system Ax = B of n linear equations
+
+A = zeros(n, n);
+b = zeros(n, 1);
+
+for i = 1:n
+   if( mod(i,2) == 0 )
+       b(i) = 2/(3*i);
+   else
+       b(i) = 0;
+   end
+   for j = 1:n
+       A(i,j) = 7/(8*(i+j+1));
+   end 
+end
+end

ENUME/references/ENUME_ProjectA_No31/task 3/rowDominant.m

@@ -0,0 +1,21 @@
+function [out] = rowDominant(A)
+%ROWDOMINANT checks row dominance of matrix A
+%   returns true when dominant, false otherwise
+n = size(A);
+
+for i = 1:1:n
+    sum = 0;
+    for j = 1:1:n
+        if i == j
+            continue
+        end
+        sum = sum + abs(A(i,j));
+    end
+    if abs(A(i,i)) <= sum
+        out = false;
+        return
+    end
+end
+out = true;
+end
+

ENUME/references/ENUME_ProjectA_No31/task 3/task3.m

@@ -0,0 +1,56 @@
+clear all;
+%First case
+A = [20 10 8 1;
+    1 10 6 2;
+    2 1 8 4;
+    1 1 1 4;];
+b = [100; 80; 60; 40];
+%From task 2
+[Aa, ba] = matrixGen2a( 10 );
+[Ab, bb] = matrixGen2b( 10 );
+
+[xJ, errorsJ] = JacobiMethod(A, b);
+[xGS, errorsGS] = GaussSeidelMethod(A, b);
+[xJa, errorsJa] = JacobiMethod(Aa, ba);
+[xGSa, errorsGSa] = GaussSeidelMethod(Aa, ba);
+[sr] = JacobiMethod(Ab, bb);
+
+%Plots
+figure(1);
+plot(1:1:length(errorsJ), errorsJ, '.')
+title("Jacobi method, matrices 3");
+xlabel("Iteration"); ylabel("Euclidean norm of errors");
+saveas(1, "./plots/Jacobi3.fig");
+saveas(1, "./plots/Jacobi3.png");
+
+figure(2);
+plot(1:1:length(errorsJa), errorsJa, '.')
+title("Jacobi method, matrices 2a");
+xlabel("Iteration"); ylabel("Euclidean norm of errors");
+saveas(2, "./plots/Jacobi2a.fig");
+saveas(2, "./plots/Jacobi2a.png");
+
+figure(3);
+plot(1:1:length(errorsGS), errorsGS, '.')
+title("Gauss-Seidel method, matrices 3");
+xlabel("Iteration"); ylabel("Euclidean norm of errors");
+saveas(3, "./plots/GaussSeidel3.fig");
+saveas(3, "./plots/GaussSeidel3.png");
+
+figure(4)
+plot(1:1:length(errorsGSa), errorsGSa, '.')
+title("Gauss-Seidel method, matrices 2a");
+xlabel("Iteration"); ylabel("Euclidean norm of errors");
+saveas(4, "./plots/GaussSeidel2a.fig");
+saveas(4, "./plots/GaussSeidel2a.png");
+
+figure(5);
+hold on;
+plot(1:1:length(errorsJ), errorsJ, 'o')
+plot(1:1:length(errorsGS), errorsGS, 'o')
+legend('Jacobi', 'Gauss-Seidel');
+title("Jacobi and Gauss-Seidel comparison");
+xlabel("Iteration"); ylabel("Euclidean norm of errors");
+xlim([0 40]);
+saveas(5, "./plots/compare.fig");
+saveas(5, "./plots/compare.png");

ENUME/references/ENUME_ProjectA_No31/task 4/QRfactorize.m

@@ -0,0 +1,26 @@
+function [Q,R] = QRfactorize(A)
+%QRFACTORIZE factorizes the matrix A using QR factorization
+%   returns the matrices Q and R
+[m, n] = size(A);
+Q = zeros(m, n);
+R = zeros(m, n);
+d = zeros(m, n);
+%Factorization
+for i = 1:1:n
+    Q(:, i) = A(:, i);
+    R(i, i) = 1;
+    d(i) = Q(:, i)' * Q(:, i);
+    
+    for j = i+1:1:n
+        R(i, j) = (Q(: ,i)' * A(:, j)) / d(i);
+        A(:, j) = A(:, j) - R(i, j) * Q(:, i);
+    end
+end
+%Normalization
+for i = 1:1:n
+    dd = norm(Q(:, i));
+    Q(:, i) = Q(:, i) / dd;
+    R(i, i:n) = R(i, i:n) * dd;
+end
+end
+

ENUME/references/ENUME_ProjectA_No31/task 4/eigenvalueQRnoshift.m

@@ -0,0 +1,17 @@
+function [eigenvalues,iteration,Afinal] = eigenvalueQRnoshift(A)
+%EIGENVALUEQRNOSHIFT calculates eigenvalues using the QR method with no
+%shifts
+%   returns eigenvalues, the number of iterations and the transformed
+%   matrix A
+iteration = 1;
+while max(max(A-diag(diag(A)))) > 10^-6
+    [Q, R] = QRfactorize(A);
+    A = R * Q;
+    iteration = iteration+1;
+end
+
+eigenvalues = diag(A);
+Afinal = A;
+
+end
+

ENUME/references/ENUME_ProjectA_No31/task 4/eigenvalueQRshift.m

@@ -0,0 +1,37 @@
+function [eigenvalues,iteration,Afinal] = eigenvalueQRshift(A)
+%EIGENVALUEQRSHIFT calculates eigenvalues using the QR method with shifts
+%   returns eigenvalues, the number of iterations and the transformed
+%   matrix A
+n = size(A, 1);
+eigenvalues = diag(ones(n));
+    
+initialSub = A;
+iteration = 0;
+    for k = n:-1:2
+        DK = initialSub;
+        
+        while max(abs(DK(k, 1:k-1))) > 10^-6
+            DD = DK(k-1:k, k-1:k);
+            [ev1, ev2] = quadpolynroots(1, -(DD(1,1) + DD(2,2)), DD(2,2) * DD(1,1) - DD(2,1) * DD(1,2));
+            
+            if abs(ev1 - DD(2, 2)) < abs(ev2 - DD(2, 2))
+                shift = ev1;
+            else
+                shift = ev2;
+            end
+            DP = DK - eye(k) * shift;
+            [Q, R] = QRfactorize(DP);
+            DK = R * Q + eye(k) * shift;
+            iteration = iteration + 1;
+        end
+        eigenvalues(k) = DK(k, k);
+        A(1:k, 1:k) = DK(1:k, 1:k);
+        if k > 2
+            initialSub = DK(1:k-1, 1:k-1);
+        else
+            eigenvalues(1) = DK(1, 1);
+        end
+    end
+    Afinal = A; 
+end
+

ENUME/references/ENUME_ProjectA_No31/task 4/quadpolynroots.m

@@ -0,0 +1,18 @@
+function [x1, x2] = quadpolynroots(a,b,c)
+%UNTITLED4 Summary of this function goes here
+%   Detailed explanation goes here
+
+l1 = -b + sqrt(b*b - 4*a*c);
+l2 = -b - sqrt(b*b - 4*a*c);
+
+if abs(l1) > abs(l2)
+    ctr = l1;
+else
+    ctr = l2;
+end
+
+x1 = ctr/(2 * a);
+x2 = ((-b) / a) - x1;
+
+end
+

ENUME/references/ENUME_ProjectA_No31/task 4/task4.m

@@ -0,0 +1,11 @@
+clear all;
+
+A = [23 12 3 5 10;
+    12 14 8 5 22;
+    3 8 9 13 11;
+    5 5 13 10 17;
+    10 22 11 17 25];
+
+[eigNS, iteNS, finNS] = eigenvalueQRnoshift(A);
+[eigS, iteS, finS] = eigenvalueQRshift(A);
+eigE = eig(A);

ENUME/references/MNUM/Projekt1/Zadanie1/Zadanie1.m

@@ -0,0 +1,10 @@
+function [] = Zadanie1()
+    x = 1.5; %poczštkowe zainicjalizowanie zmiennej
+    g = 1.0;
+    while( x > 1 ) %przechodzenie przez pętle i dzielenie epsilona 
+        g = g/2;   %tak długo aż dodanie go nie wpłynie na wynik
+        x = 1.0 + g;
+    end
+    g = g*2; %jeden przebieg pętli w tył
+    fprintf('Wyznaczony epsilon maszynowy: %d \n',g);
+end

ENUME/references/MNUM/Projekt1/Zadanie2/Zadanie2a.m

@@ -0,0 +1,16 @@
+function [] = Zadanie2a(n)
+    A = a_genA(n);          %generator macierzy wejściowej A
+    b = a_genB(n);          %generator macierzy wejściowej b
+    b=b';
+    
+    tic                     %rozpoczęcie pomiaru czasu
+    L1 = cholesky(A,n);     %rozkład metodą Cholesky'ego-Banachiewicza
+    x = solveEq(L1,b);      %rozwiązanie równania
+    t = toc;                %koniec pomiaru czasu
+    
+    r = b - A * x;          %obliczenie residuum
+    br = norm(r);           %obliczenie błędu rozwiązania jako normy z residuum
+    
+    fprintf("Zmierzony czas rozwiązania: %d \n",t);
+    fprintf("Liczba równań i błąd rozwiązania: \n %d %d \n",n,br);
+end

ENUME/references/MNUM/Projekt1/Zadanie2/Zadanie2b.m

@@ -0,0 +1,16 @@
+function [] = Zadanie2b(n)
+    A = b_genA(n);
+    b = b_genB(n);
+    b=b';
+    
+    tic                     %rozpoczęcie pomiaru czasu
+    L1 = cholesky(A,n);     %rozkład metodš Cholesky'ego-Banachiewicza
+    x = solveEq(L1,b);      %rozwišzanie równania
+    t = toc;                %koniec pomiaru czasu
+    
+    r = b - A * x;          %obliczenie residuum
+    br = norm(r);           %obliczenie błędu rozwišzania jako normy z residuum
+    
+    fprintf("Zmierzony czas rozwišzania: %d \n",t);
+    fprintf("Liczba równań i błšd rozwišzania: \n %d %d \n",n,br);
+end

ENUME/references/MNUM/Projekt1/Zadanie2/Zadanie2c.m

@@ -0,0 +1,16 @@
+function [] = Zadanie2c(n)
+    A = c_genA(n);
+    b = c_genB(n);
+    b=b';
+    
+    tic                     %rozpoczęcie pomiaru czasu
+    L1 = cholesky(A,n);     %rozkład metodš Cholesky'ego-Banachiewicza
+    x = solveEq(L1,b);      %rozwišzanie równania
+    t = toc;                %koniec pomiaru czasu
+    
+    r = b - A * x;          %obliczenie residuum
+    br = norm(r);           %obliczenie błędu rozwišzania jako normy z residuum
+    
+    fprintf("Zmierzony czas rozwišzania: %d \n",t);
+    fprintf("Liczba równań i błšd rozwišzania: \n %d %d \n",n,br);
+end

ENUME/references/MNUM/Projekt1/Zadanie2/a/a_genA.m

@@ -0,0 +1,5 @@
+function [A] = a_genA(n)
+    v1 = ones(1,n)*10;
+    v2 = ones(1,n-1)*4;
+    A = diag(v1) + diag(v2,1) + diag(v2,-1);
+end

ENUME/references/MNUM/Projekt1/Zadanie2/a/a_genB.m

@@ -0,0 +1,4 @@
+function [b] = a_genB(n)
+    bottom = 2.5 - 0.5*n;
+    b = 2:-0.5:bottom;
+end

ENUME/references/MNUM/Projekt1/Zadanie2/b/b_genA.m

@@ -0,0 +1,12 @@
+function [A] = b_genA(n)
+    v1 = 4*n^2-1:-1:4*n^2-n;
+    A = diag(v1);
+        for i = 1:n-1
+            j=i+1;
+            while(j <= n)
+                A(i,j) = 2*(i+j)+1; %dodanie wartości do odpowiedniego indeksu macierzy
+                j=j+1;
+            end
+        end
+    A = triu(A)+triu(A,1)'; %odbicie symetrzyczne wobec przekątnej
+end

ENUME/references/MNUM/Projekt1/Zadanie2/b/b_genB.m

@@ -0,0 +1,4 @@
+function [b] = b_genB(n)
+    top = 2.5 + 0.6*n;
+    b = 3.1:0.6:top;
+end

ENUME/references/MNUM/Projekt1/Zadanie2/backwardSubstitution.m

@@ -0,0 +1,7 @@
+function [x] = backwardSubstitution(A,b)
+    m = length(b);                                      %wyłuskanie ostatniego indeksu
+    x(m,1) = b(m)/A(m,m);                               %wyliczenie ostatniego elementu
+    for i = m-1:-1:1                                    %iteracja po wszystkiech równaniach poczšwszy od przedostatniego wiersza
+        x(i,1)=(b(i)-A(i,i+1:m)*x(i+1:m,1))./A(i,i);    %wyznaczenie niewiadomych
+    end
+end

ENUME/references/MNUM/Projekt1/Zadanie2/c/c_genA.m

@@ -0,0 +1,12 @@
+function [A] = c_genA(n)
+    v1 = 0.2*n+0.3:0.3:0.2*n+0.3*n;
+    A = diag(v1);
+        for i = 1:n-1
+            j=i+1;
+            while(j <= n)
+                A(i,j) = 1/(4*(i+j+1)); %dodanie wartości do odpowiedniego indeksu macierzy
+                j=j+1;
+            end
+        end
+    A = triu(A)+triu(A,1)'; %odbicie symetryczne wobec przekątnej
+end

ENUME/references/MNUM/Projekt1/Zadanie2/c/c_genB.m

@@ -0,0 +1,4 @@
+function [b] = c_genB(n)
+    b = 1:1:n;
+    b = (1./b).*(5/3); %odwrócenie i przemnożenie każdego elementu
+end

ENUME/references/MNUM/Projekt1/Zadanie2/cholesky.m

@@ -0,0 +1,12 @@
+function[L1] = cholesky(A,n)    
+    L1 = zeros(n,n);
+    for i = 1:n                                                 %przechodzenie po kolumnach
+        for j = i:n                                             %przechodzenie wierszy wewnątrz kolumny
+            if(i == j)
+                L1(i,i) = sqrt(A(i,i)-sumDiag(L1,i));           %przypisanie wartości na przekątnej macierzy
+            else
+                L1(j,i) = (A(j,i)-sumRest(L1,i,j))/L1(i,i);     %przypisanie wartości w reszcie przypadków
+            end
+        end
+    end
+end

ENUME/references/MNUM/Projekt1/Zadanie2/forwardSubstitution.m

@@ -0,0 +1,7 @@
+function [x] = forwardSubstitution(A,b)
+    m = length(b);                                      %wyłuskanie ostatniego indeksu
+    x(1,1) = b(1)/A(1,1);                               %wyliczenie pierwszego elementu
+    for i = 2:m                                         %iteracja po wszystkich równaniach poczšwszy od drugiego wiersza
+        x(i,1)=(b(i)-A(i,1:i-1)*x(1:i-1,1))./A(i,i);    %wyznaczanie niewiadomych
+    end
+end

ENUME/references/MNUM/Projekt1/Zadanie2/solveEq.m

@@ -0,0 +1,5 @@
+function [x] = solveEq(L1,b)
+    b = b';                             %transponowanie macierzy
+    y = forwardSubstitution(L1,b);      
+    x = backwardSubstitution(L1',y);
+end

ENUME/references/MNUM/Projekt1/Zadanie2/sumDiag.m

@@ -0,0 +1,6 @@
+function[sum] = sumDiag(L1,i)
+    sum=0;
+    for k = 1:1:i-1
+       sum=sum+L1(i,k)^2;
+    end
+end

ENUME/references/MNUM/Projekt1/Zadanie2/sumRest.m

@@ -0,0 +1,6 @@
+function[sum] = sumRest(L1,i,j)
+    sum=0;
+    for k = 1:1:i-1
+        sum=sum+(L1(j,k)*L1(i,k));
+    end
+end

ENUME/references/MNUM/Projekt1/Zadanie3/GaussSeidel.m

@@ -0,0 +1,21 @@
+function[x,iter,od] = GaussSeidel(A,b,n,e) 
+    x = zeros(n,1);
+    
+    %Stworzenie macierzy naddiagonalne, poddiagonalnej i diagonalnej
+    [L,D,U] = rozkladGaussSeidel(A,n);
+    b=b';
+    r = 1;
+    iter = 1;                                  
+    while(r>e || norm(A*x-b)>e)                                   %kolejne iteracje
+        y = x;                                  %zapamietujemy wektor x z poprzedniej iteracji
+        w = U*x - b;
+        for i = 1:n
+            x(i) = (-L(i,:)*x- w(i))/D(i,i);
+        end
+    r = norm(x-y);                                %liczmy blad z nomrmy euklidesowej, po kazdej iteracji
+    iter = iter + 1;                              %zliczamy ilosc iteracji
+    
+    end
+    %fprintf('norm(A*x-b):%d\n',norm(A*x-b));
+    od=norm(A*x-b);
+end

ENUME/references/MNUM/Projekt1/Zadanie3/Zadanie3.m

@@ -0,0 +1,24 @@
+function[x] = Zadanie3(n,e2)
+    i=1;
+    e=1;
+    A = b_genA(n);
+    b = b_genB(n);
+    
+    %Sprawdzenie warunku dostatecznego zbieżności
+    if(warDost(A,n) == 0)
+        disp('Warunek silnej dominacji diagonalnej nie jest spelniony');
+        return
+    end
+    
+    while e>e2 % 0.0000001
+        [x,iter(i),od] = GaussSeidel(A,b,n,e);
+        e = e/10;
+        r(i) = od;
+        i = i+1;
+    end
+    plot(iter, r);
+    disp(norm(A*x-b));
+    title('Zaleznosc bledu wyniku od ilosci iteracji')
+    xlabel('Ilosc iteracji');
+    ylabel('Blad wyniku');
+end

ENUME/references/MNUM/Projekt1/Zadanie3/rozkladGaussSeidel.m

@@ -0,0 +1,15 @@
+function[L,D,U] = rozkladGaussSeidel(A,n)
+    L = zeros(n,n); %macierz poddiagonalna
+    U = zeros(n,n); %macierz naddiagonalna
+    D = diag(diag(A)); %macierz diagonalna
+    
+    for i = 1:n
+        for j = 1:n
+            if(i<j)
+                U(i,j) = A(i,j);
+            elseif(i>j)
+                L(i,j) = A(i,j);
+            end
+        end
+    end
+end

ENUME/references/MNUM/Projekt1/Zadanie3/warDost.m

@@ -0,0 +1,17 @@
+function[bool] = warDost(A,n)
+    for i = 1:n
+        %Sumujemy wszystkie elementy w wierszu oprócz głównej przekštnej
+        sum = 0;
+        for j = 1:n
+            if i~=j
+                sum = sum + abs(A(i,j));
+            end
+        end
+        %Sprawdzamy warunek silnej dominancji dla jednego wiersza 
+        if sum > abs(A(i,i))
+            bool = 0;
+            return;
+        end
+    end
+    bool = 1;
+end

ENUME/references/MNUM/Projekt2/Zadanie1/Zadanie1.m

@@ -0,0 +1,73 @@
+function [] = Zadanie1(czyPrzesun,czySym)
+    % czyPrzesun - wartosc 1 jesli wersja z przesunieciami
+    % czySym - wartosc 1 jestli macierz symetryczna
+
+    iterQR = 0; %liczba wartosci wlasnej QR bez przesuniec
+    iterQRS = 0; %liczba wartosci wlasnej QR z przesunieciami
+    
+    timeQR = 0; %czas wykonania alg QR bez przesuniec
+    timeQRS = 0; %czas wykonania alg QR z przesunieciami
+    timeE = 0; %czas wykonania alg eig() wbudowanym w MATLAB
+    
+    SIZE = 20; %rozmiar maciery
+    matrixNumberR = 30; %ilosc losowych macierzy
+    
+    ILEMACQR = 0;
+    ILEMACQRS = 0;
+
+    %Generowanie danych
+    for i=1:matrixNumberR
+        
+        A = genA(SIZE,czySym);
+        tolerance = 0.0000001;
+        imax = 200;
+        
+        start = tic;
+        [~,D] = eig(A);
+        timeEig = toc(start);
+        timeE = timeE + timeEig;
+        
+        if czyPrzesun == 1 %QR z przesunieciami (eignes - wartWlasne)
+            [eigens, iteracje, time, ok] = zPrzesun(A, tolerance, imax);
+            if ok == 1
+                ILEMACQRS = ILEMACQRS + 1;
+                iterQRS = iterQRS + iteracje;
+                timeQRS = timeQRS + time;
+            end
+        else %QR bez przesuniec  (eigens - wartWlasne)
+            [eigens, iteracje, time, ok] = bezPrzesun(A, tolerance, imax);
+            if ok == 1
+                ILEMACQR = ILEMACQR + 1;
+                iterQR = iterQR + iteracje;
+                timeQR = timeQR + time;
+            end 
+        end
+    end
+    
+    fprintf('Ilosc macierzy: %d\n',matrixNumberR);
+    fprintf('Wielkosc macierzy: %d\n',SIZE);
+    SREDNIAczasEig = timeE / matrixNumberR;
+    
+    if czyPrzesun == 1
+        %Wyniki z przesunieciami:
+        SREDNIAQRSI = iterQRS / ILEMACQRS;
+        SREDNIAczasQRS = timeQRS / ILEMACQRS;
+        fprintf('Z przesunieciami:\n');
+        fprintf('Ilosc zkonczonych sukcesem: %d\n', ILEMACQRS);
+        fprintf('Srednia ilosci iteracji %d\n',SREDNIAQRSI);
+        fprintf('Sredni czas obliczen %d\n',SREDNIAczasQRS);
+    else
+        %Wyniki bez przesuniec:
+        SREDNIAQRI = iterQR / ILEMACQR;
+        SREDNIAczasQR = timeQR / ILEMACQR;
+        fprintf('Bez przesuniec:\n');
+        fprintf('Ilosc zkonczonych sukcesem: %d\n', ILEMACQR);
+        fprintf('Srednia ilosci iteracji %d\n',SREDNIAQRI);
+        fprintf('Sredni czas obliczen %d\n',SREDNIAczasQR);
+    end
+    
+    fprintf('Sredni czas obliczen eig %d\n',SREDNIAczasEig);
+    d = diag(D);
+    disp(sort(d));
+    disp(sort(eigens));
+end

ENUME/references/MNUM/Projekt2/Zadanie1/bezPrzesun.m

@@ -0,0 +1,20 @@
+function [wartWlasne, iteracje, time, ok] = bezPrzesun(D, tolerance, imax)
+    %imax - maksymalna liczba iteracji
+    % ok - czy funkcja wykonala sie z przekroczeniem imax
+    start = tic;
+    ok = 1; 
+    i = 1; 
+    
+    while i <= imax && max(max(D-diag(diag(D)))) > tolerance
+        [Q1, R1] = qrZmodGS(D);
+        D = R1*Q1; %macierz po przekształceniu
+        i = i + 1;
+    end
+    
+    if i > imax 
+        ok = 0;
+    end
+    iteracje = i;
+    wartWlasne = diag(D); %wykstraktowanie wektora wartości własnych
+    time = toc(start);
+end

ENUME/references/MNUM/Projekt2/Zadanie1/genA.m

@@ -0,0 +1,10 @@
+function [A] = genA(SIZE,czySym)
+        A = rand(SIZE);
+        while rank(A)~= SIZE %powstanie macierzy o pelnym rzedzie
+            A = rand(SIZE);
+        end
+        
+        if czySym == 1
+            A = A'+A; %dla symetrycznych
+        end
+end

ENUME/references/MNUM/Projekt2/Zadanie1/qrZmodGS.m

@@ -0,0 +1,25 @@
+%rozkład QR (waski) macierzy zmodyfikowanym alg. Grama-Schmidta dla 
+%macierzy mxn (m>=n) o pelnym rzedzie, rzeczywistej lub zespolonej 
+%na podstawie ksišzki prof. Tatjewskiego
+function [Q,R] = qrZmodGS(A)
+    [m n] = size(A);
+    Q = zeros(m,n); 
+    R = zeros(n,n);
+    d = zeros(1,n);
+    %rozkład A z kolumnami Q ortogonalnymi
+    for i=1:n
+        R(i,i)=1;
+        Q(:,i)=A(:,i);
+        d(i)=Q(:,i)'*Q(:,i);
+        for j=i+1:n
+            R(i,j)=(Q(:,i)'*A(:,j))/d(i);
+            A(:,j)=A(:,j)-R(i,j)*Q(:,i);
+        end
+    end
+    %normowanie rozkladu (kolumny Q ortogonalne)
+    for i=1:n
+        dd = norm(Q(:,i));
+        Q(:,i) = Q(:,i)/dd;
+        R(i,i:n) = R(i,i:n)*dd;
+    end
+end

ENUME/references/MNUM/Projekt2/Zadanie1/zPrzesun.m

@@ -0,0 +1,44 @@
+function [wartWlasne, iteracje, time, ok] = zPrzesun(A,tolerance,imax)
+    % ok - czy funkcja wykonala sie z przekroczeniem imax
+    % tolerance - toleracnaj jako górna granica wartości elementów zerowanych
+    start = tic;
+    ok = 1; 
+    n = size(A,1);
+    wartWlasne = diag(zeros(n));
+    iteracje = 0;
+    INITIALsubmatrix = A;
+    for k=n:-1:2
+        DK = INITIALsubmatrix(1:k, 1:k); %macierz potrzebna do wyznaczenia wartosci w1
+        i = 0;
+       
+        while i <= imax && max(abs(DK(k,1:k-1))) > tolerance
+            
+            ev = roots([1, -(DK(k-1,k-1)+DK(k,k)), DK(k,k)*DK(k-1,k-1)-DK(k,k-1)*DK(k-1,k)]);
+            if abs(ev(1)-DK(k,k)) < abs(ev(2)-DK(k,k))
+                shift = ev(1); % nasze przesuniecie jako najbliższa DK(k,k) wartosc 
+                               % wlasna analizowanej macierzy 2x2
+            else
+                shift = ev(2);
+            end
+            DK = DK - eye(k)*shift; %nasza macierz przesunięta
+            [Q1, R1] = qrZmodGS(DK); %faktoryzacja QR
+            DK = R1*Q1 + eye(k)*shift; %macierz przekształcona
+            i = i+1;
+            
+            iteracje = iteracje + 1;
+        end
+        
+        if i > imax
+            ok = 0;
+            break;
+        end
+        wartWlasne(k) = DK(k,k);
+        
+        if k>2
+            INITIALsubmatrix = DK(1:k-1,1:k-1); %definicja macierzy
+        else
+            wartWlasne(1) = DK(1,1); %ostatnia wartosc wlasna
+        end
+    end
+    time = toc(start);
+end

ENUME/references/MNUM/Projekt2/Zadanie2/Zadanie2.m

@@ -0,0 +1,25 @@
+function [] = Zadanie2(STOPIEN)
+    x = -5:1:5;
+    x2 = -5:0.1:5;
+    y = [-7.7743 -0.2235 1.9026 0.6572 0.1165 -1.8144 -1.0968 -0.8261 1.3327 6.1857 8.2891];
+    x=x';
+    y=y';
+    a = najmnKwadratow(STOPIEN,x,y);
+  
+    eukNorm = euklidesNorm(a,y,x);
+    czebNorm = czebyszewNorm(a,y,x);
+    
+    fprintf('Błšd aproksymacji w normie Czebyszewa dla równania rzędu %d wynosi %d\n',STOPIEN,czebNorm);
+    fprintf('Błšd aproksymacji w normie Euklidesowej dla równania rzędu %d wynosi %d\n',STOPIEN,eukNorm);
+    
+    %Rysowanie wykresu
+    scatter(x,y,'filled')
+    hold on
+    p = polyval(a,x2);
+    plot(x2,p);
+    grid;
+    title('Aproksymacja przy stopniu równym 10 - równania normalne');
+    xlabel('x');
+    ylabel('y');
+    legend('Punkty oryginalne','funkcja aproksymujšca')
+end

ENUME/references/MNUM/Projekt2/Zadanie2/czebyszewNorm.m

@@ -0,0 +1,13 @@
+function [maxError] = czebyszewNorm(a,f,x)
+    p = polyval(a,x);
+    [s,~] = size(f);
+    i = 1;
+    maxError = 0;
+    while i<=s
+        err = abs(p(i)-f(i));
+        if err > maxError 
+            maxError = err;
+        end
+        i = i + 1;
+    end
+end

ENUME/references/MNUM/Projekt2/Zadanie2/euklidesNorm.m

@@ -0,0 +1,11 @@
+function [err] = euklidesNorm(a,f,x)
+    p = polyval(a,x);
+    [s,~] = size(f);
+    i = 1;
+    sum=0;
+    while i<=s
+        sum = sum + (p(i)-f(i))^2;
+        i = i + 1;
+    end
+    err = sqrt((1/s)*sum);
+end

ENUME/references/MNUM/Projekt2/Zadanie2/najmnKwadratow.m

@@ -0,0 +1,24 @@
+function [a] = najmnKwadratow(n, x, y)
+    % n - zadany stopien wielomianu
+    % x - wektor argumentow
+    % y - wektor wartosci
+    
+    % a - wektor wyznaczonych wspolczynnikow
+    
+    
+    [m, ~] = size(x); %pobranie rozmiarow
+    A = zeros(m, n+1);
+    
+    for i=1:m %wiersze
+        for j=0:n %kolumny
+            A(i,n+1-j) = x(i)^(j); %wypełniamy odpowiednimi wartooeciami potegi x
+        end
+    end
+    
+    %Rozkład QR
+    [Q, R] = qrZmodGS(A); %rozklad qr metoda g-s z poprzedniego zadania
+    a = R \ (Q'*y);
+    
+    %Uklad równań normalnych
+    %a = (A'*A)\(A'*y); %rozwiazanie ukaadu równan  
+end