From e461df7573c2b7b7e26c965d8cf2d8e175d67378 Mon Sep 17 00:00:00 2001 From: Nguyễn Gia Phong Date: Mon, 16 Dec 2019 21:31:18 +0700 Subject: [usth/MATH2.2] Numerical Methods The future starts now. --- usth/MATH2.2/labwork/3/GaussPivot.m | 25 +++++++++++++++++++++++++ usth/MATH2.2/labwork/3/LU.m | 10 ++++++++++ usth/MATH2.2/labwork/3/LU_pivot.m | 20 ++++++++++++++++++++ usth/MATH2.2/labwork/3/gauss.m | 12 ++++++++++++ usth/MATH2.2/labwork/3/labwork3.pdf | Bin 0 -> 72465 bytes usth/MATH2.2/labwork/3/subst.m | 17 +++++++++++++++++ 6 files changed, 84 insertions(+) create mode 100644 usth/MATH2.2/labwork/3/GaussPivot.m create mode 100644 usth/MATH2.2/labwork/3/LU.m create mode 100644 usth/MATH2.2/labwork/3/LU_pivot.m create mode 100644 usth/MATH2.2/labwork/3/gauss.m create mode 100644 usth/MATH2.2/labwork/3/labwork3.pdf create mode 100644 usth/MATH2.2/labwork/3/subst.m (limited to 'usth/MATH2.2/labwork/3') diff --git a/usth/MATH2.2/labwork/3/GaussPivot.m b/usth/MATH2.2/labwork/3/GaussPivot.m new file mode 100644 index 0000000..7764b3a --- /dev/null +++ b/usth/MATH2.2/labwork/3/GaussPivot.m @@ -0,0 +1,25 @@ +function x = GaussPivot (A, b) + [m n] = size (A); + m ~= n && error ('Matrix A must be square'); + nb = n + 1; + Aug = [A b]; + % forward elimination + for k = 1 : n - 1 + % partial pivoting + [big i] = max (abs (Aug(k:n, k))); + ipr = i + k - 1; + if ipr ~= k + Aug([k ipr], :) = Aug([ipr k], :); + end + for i = k + 1 : n + coeff = Aug(i, k) / Aug(k, k); + Aug(i, k:nb) = Aug(i, k:nb) - coeff * Aug(k, k:nb); + end + end + % back substitution + x = zeros (n, 1); + x(n) = Aug(n, nb) / Aug(n, n); + for i = n - 1 : -1 : 1 + x(i) = (Aug(i, nb) - Aug(i, i+1:n) * x(i+1:n)) / Aug(i, i); + end +end diff --git a/usth/MATH2.2/labwork/3/LU.m b/usth/MATH2.2/labwork/3/LU.m new file mode 100644 index 0000000..860486d --- /dev/null +++ b/usth/MATH2.2/labwork/3/LU.m @@ -0,0 +1,10 @@ +function [L U] = LU (A) + [n n1] = size (A); + [L U] = deal (eye (n), A); + for k = 1:n + for i = k + 1 : n + L(i, k) = U(i, k) / U(k, k); + U(i, :) = U(i, :) - L(i, k)*U(k, :); + end + end +end diff --git a/usth/MATH2.2/labwork/3/LU_pivot.m b/usth/MATH2.2/labwork/3/LU_pivot.m new file mode 100644 index 0000000..c506697 --- /dev/null +++ b/usth/MATH2.2/labwork/3/LU_pivot.m @@ -0,0 +1,20 @@ +function [L U P] = LU_pivot (A) + [n _] = size (A); + [L P U] = deal (eye (n), eye (n), A); + for k = 1:n + [pivot m] = max (abs (U(k:n, k))); + m = m + k - 1; + if (m ~= k) + U([m k], :) = U([k m], :); % interchange rows m and k in U + P([m k], :) = P([k m], :); % interchange rows m and k in P + if k >= 2; % very important point + % interchange rows m and k in columns 1:k-1 of L + L([m k], 1:k-1) = L([k m], 1:k-1); + end + end + for i = k + 1 : n + L(i, k) = U(i, k) / U(k, k); + U(i, :) -= L(i, k)*U(k, :); + end + end +end diff --git a/usth/MATH2.2/labwork/3/gauss.m b/usth/MATH2.2/labwork/3/gauss.m new file mode 100644 index 0000000..09128dc --- /dev/null +++ b/usth/MATH2.2/labwork/3/gauss.m @@ -0,0 +1,12 @@ +function x = gauss (A, b) + issquare (A) || error ('Matrix A must be square'); + Aug = [A b]; + [m n] = size (Aug); + % forward elimination + for k = 1 : m - 1 + for l = k + 1 : m + Aug(l, k:n) = Aug(l, k:n) - Aug(k, k:n)*Aug(l, k)/Aug(k, k); + end + end + x = ubst (Aug); +end diff --git a/usth/MATH2.2/labwork/3/labwork3.pdf b/usth/MATH2.2/labwork/3/labwork3.pdf new file mode 100644 index 0000000..4e0c7e0 Binary files /dev/null and b/usth/MATH2.2/labwork/3/labwork3.pdf differ diff --git a/usth/MATH2.2/labwork/3/subst.m b/usth/MATH2.2/labwork/3/subst.m new file mode 100644 index 0000000..5c6e783 --- /dev/null +++ b/usth/MATH2.2/labwork/3/subst.m @@ -0,0 +1,17 @@ +function x = subst (aug) + [m n] = size (aug); + x = zeros (m, 1); + if (istril (aug)) + x(m) = aug(1, n) / aug(1, m); + for k = 2 : m + x(k) = (aug(k, n) - aug(k, k+1:m)*x(1:k-1)) / aug(k, k); + end + elseif (istriu (aug)) + x(m) = aug(m, n) / aug(m, m); + for k = m - 1 : -1 : 1 + x(k) = (aug(k, n) - aug(k, k+1:m)*x(k+1:m)) / aug(k, k); + end + else + error ('aug must be a triangular matrix'); + end +end -- cgit 1.4.1