diff --git a/w3/gaussLR.m b/w3/gaussLR.m new file mode 100644 index 0000000..b2de143 --- /dev/null +++ b/w3/gaussLR.m @@ -0,0 +1,32 @@ +function A = gaussLR(A) + % Gaussian Elimination for reforming the matrix into left-lower and right-upper triangular matrix + % A: A R^(n \times n) matrix + + % actually not used here (but it was planned to be used): bsxfun, which makes this possible: bsxfun(@times, [1,2;3,4], [1;2]) = [1,2;6,8]. It multiplies the first row by 1, and the second one by 2. + % newer (and actually (still not) used here): [1,2;3,4].*[1;2] = [1,2;6,8] + + % for convenience + len = length(A) + + for i = 1:len + prefactor = A(i,i) + if prefactor == 0 + error("LR-Zerlegung failed. Please use pivotization, or check that the matrix is regular!") + end + + % elimination (find R) + factors = zeros(len, 1); + for j = (i+1):len + factors(j) = A(j, i) / prefactor; + end + + % A(:,i:len) describes the matrix A with only the last i cols (so, for example, a 3x3 matrix with (:,2,len) with len=3) is the last 2 cols of the matrix, and with that a 2x3 matrix + + % extra step for traceability, could also be compressed to a single line + multmatrix = (-factors * A(i,i:len)) + A(:,i:len) = A(:,i:len) + multmatrix; + + % save factors for L + A(:,i) = A(:,i) + factors + end +end \ No newline at end of file diff --git a/w3/lrPivot.m b/w3/lrPivot.m new file mode 100644 index 0000000..f720590 --- /dev/null +++ b/w3/lrPivot.m @@ -0,0 +1,47 @@ +function [A, P] = lrPivot(A) + % Gaussian Elimination for reforming the matrix into left-lower and right-upper triangular matrix + % A: A R^(n \times n) matrix + % P: Permutations-Matrix + + % actually not used here (but it was planned to be used): bsxfun, which makes this possible: bsxfun(@times, [1,2;3,4], [1;2]) = [1,2;6,8]. It multiplies the first row by 1, and the second one by 2. + % newer (and actually (still not) used here): [1,2;3,4].*[1;2] = [1,2;6,8] + + % for convenience + len = length(A) + pvec = 1:len + + + % Generate LR-Zerlegung + for i = 1:len + + % pivot + [maxVal, maxValPos] = max(abs(A(pvec(i:end), i))) + + % swap + temp = pvec(i) + pvec(i) = pvec(maxValPos + i - 1) + pvec(maxValPos + i - 1) = temp + + prefactor = A(pvec(i),i) + if prefactor == 0 + error("LR-Zerlegung failed. Please check that the matrix is regular!") + end + + % elimination (find R) + factors = zeros(len, 1) + for j = (i+1):len + factors(pvec(j)) = A(pvec(j), i) / prefactor + end + + % A(:,i:len) describes the matrix A with only the last i cols (so, for example, a 3x3 matrix with (:,2,len) with len=3) is the last 2 cols of the matrix, and with that a 2x3 matrix + + % extra step for traceability, could also be compressed to a single line + multmatrix = (-factors * A(pvec(i),i:len)) + A(:,i:len) = A(:,i:len) + multmatrix + + % save factors for L + A(:,i) = A(:,i) + factors + end + + P = eye(len)(:,pvec) +end \ No newline at end of file diff --git a/w3/solveLR.m b/w3/solveLR.m new file mode 100644 index 0000000..804322b --- /dev/null +++ b/w3/solveLR.m @@ -0,0 +1,30 @@ +function x_ret = solveLR(A, b) + len = length(b) + z = zeros(len,1) + + % solve Lz = b + for i = 1:len + z(i) = b(i); + + for j = 1:(i-1) + z(i) = z(i) - z(j) * A(i, j); + end + + % no division necessary as L is unipotent + end + + x = zeros(len,1) + + % solve Rx = z + for i = len:-1:1 + x(i) = z(i); + + for j = len:-1:(i+1) + x(i) = x(i) - x(j) * A(i, j); + end + + x(i) = x(i) / A(i, i) + end + + x_ret = x +end \ No newline at end of file diff --git a/w3/solveLrPivot.m b/w3/solveLrPivot.m new file mode 100644 index 0000000..b554015 --- /dev/null +++ b/w3/solveLrPivot.m @@ -0,0 +1,30 @@ +function x = solveLrPivot(A, P, b) + len = length(b) + z = zeros(len,1) + + % solve Lz = b + for i = 1:len + z(i) = b(i) + + for j = 1:(i-1) + z(i) = z(i) - z(j) * A(P(i), j) + end + + % no division necessary as L is unipotent + end + + x = zeros(len,1) + + % solve Rx = z + for i = len:-1:1 + x(i) = z(i) + + for j = len:-1:(i+1) + x(i) = x(i) - x(j) * A(P(i), j) + end + + x(i) = x(i) / A(P(i), i) + end + + x_ret = x +end \ No newline at end of file diff --git a/w3/testGaussLR.m b/w3/testGaussLR.m new file mode 100644 index 0000000..e059d9a --- /dev/null +++ b/w3/testGaussLR.m @@ -0,0 +1,14 @@ +function testGaussLR() + success = testSingle( [ 1,1,1; 4,3,-1; 3,5,3 ], [1,1,1; 4,-1,-5; 3,-2,-10] ) + success = success & testSingle( [1,0;0,1], [1,0;0,1] ) + success = success & testSingle( [6,-4,7; -12,5,-12; 18,0,22], [6,-4,7; -2,-3,2; 3,-4,9] ) + if success + disp("It works!") + else + disp("It doesn't work :/") + end +end + +function success = testSingle(init, expected) + success = isequal(gaussLR(init), expected); +end \ No newline at end of file diff --git a/w3/testSolve.m b/w3/testSolve.m new file mode 100644 index 0000000..ffe3413 --- /dev/null +++ b/w3/testSolve.m @@ -0,0 +1,40 @@ +function testSolve + result = Test([0, 1; 1, 1], [1; 1]) + result = result & Test([11,44,1; 0.1,0.4,3; 0,1,-1], [1; 1; 1]) + result = result & Test([0.001,1,1; -1,0.004,0.004; -1000,0.004,0.000004], [1; 1; 1]) + + if result(1) % nopivot + disp("Pivot-less Solving worked!") + else + disp("Pivot-less Solving broke!") + end + + if result(2) % pivot + disp("Pivot-full Solving worked!") + else + disp("Pivot-full Solving broke!") + end +end + +function [nopivot, pivot] = Test(A, b) + try + nopivot = isequal(Solve(A, b), linsolve(A, b)) + catch + nopivot = false + end + + try + pivot = isequal(SolvePivot(A, b), linsolve(A, b)) + catch + pivot = false + end +end + +function x = Solve(A, b) + x = solveLR(gaussLR(A), b) +end + +function x = SolvePivot(A, b) + [A, T] = lrPivot(A) + x = solveLrPivot(A, T, b) +end diff --git a/w3/testSolveLR.m b/w3/testSolveLR.m new file mode 100644 index 0000000..e705aac --- /dev/null +++ b/w3/testSolveLR.m @@ -0,0 +1,14 @@ +function testSolveLR() + success = test( [1,3,2; 2,15,2; 1,3,4], [1;2;3] ) + success = success & test ( [1,1,1; 4,3,-1; 3,5,3], [4;2;0] ) + success = success & test ( [1, 0, 0; 0, 1, 0; 0, 0, 1], [4; 2; 0] ) + if success + disp("It works!") + else + disp("It broke!") + end +end + +function success = test(A, b) + success = isequal(linsolve(A, b), solveLR(gaussLR(A), b)) +end \ No newline at end of file