feat: completed w3 contents
This commit is contained in:
parent
abeb814912
commit
abb34c7834
32
w3/gaussLR.m
Normal file
32
w3/gaussLR.m
Normal file
@ -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
|
47
w3/lrPivot.m
Normal file
47
w3/lrPivot.m
Normal file
@ -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
|
30
w3/solveLR.m
Normal file
30
w3/solveLR.m
Normal file
@ -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
|
30
w3/solveLrPivot.m
Normal file
30
w3/solveLrPivot.m
Normal file
@ -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
|
14
w3/testGaussLR.m
Normal file
14
w3/testGaussLR.m
Normal file
@ -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
|
40
w3/testSolve.m
Normal file
40
w3/testSolve.m
Normal file
@ -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
|
14
w3/testSolveLR.m
Normal file
14
w3/testSolveLR.m
Normal file
@ -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
|
Loading…
Reference in New Issue
Block a user