From 1c0131eca403b439cd5ee2517622cec5293dcbcf Mon Sep 17 00:00:00 2001 From: =?utf8?q?Nils=20Forss=C3=A9n?= Date: Fri, 16 Sep 2022 14:57:49 +0200 Subject: [PATCH] pushing --- DrawCurveFit.m | 57 +++++++++++++++++++++++++++++++++++++++++++++++ gaussSeidel.asv | 35 +++++++++++++++++++++++++++++ gaussSeidel.m | 39 ++++++++++++++++++++++++++++++++ graph_nils.m | 10 +++++++++ lufact.m | 42 ++++++++++++++++++++++++++++++++++ matrix_vec_mult.m | 12 ++++++++++ pressure_data.m | 8 +++++++ thinQR.m | 38 +++++++++++++++++++++++++++++++ 8 files changed, 241 insertions(+) create mode 100644 DrawCurveFit.m create mode 100644 gaussSeidel.asv create mode 100644 gaussSeidel.m create mode 100644 graph_nils.m create mode 100644 lufact.m create mode 100644 matrix_vec_mult.m create mode 100644 pressure_data.m create mode 100644 thinQR.m diff --git a/DrawCurveFit.m b/DrawCurveFit.m new file mode 100644 index 0000000..879873b --- /dev/null +++ b/DrawCurveFit.m @@ -0,0 +1,57 @@ +function DrawCurveFit(c, t, y) +%% +% Given the original data points (t_i, y_i) and the coefficients c_j +% for the curve fitting ansatz +% f(t) = c_1 + c_2 cos(2 pi t / T) + c_3 sin(2 pi t / T) +% + c_4 cos(4 pi t / T) + c_5 sin(4 pi t / T) +% where T = t(end) - t(1). +% This function plots the curve fitting and the original points in the +% same figure. +% +% INPUT: +% c - Vector of coefficient values +% t - Vector of original time values +% y - Vector of original data values +% +% OUTPUT: +% prints a figure to the screen + +%% +% Compute length of the time interval + + T = t(end) - t(1); + +%% +% Get an enriched set of uniform points in the interval [t(1), t(end)] +% for nicer plotting + + tt = linspace(t(1), t(end), 200); + +% Compute the curve ansatz on the enriched points + + f = c(1) + c(2) * cos(2 * pi * tt / T) + c(3) * sin(2 * pi * tt / T)... + + c(4) * cos(4 * pi * tt / T) + c(5) * sin(4 * pi * tt / T); + +%% +% Make a new figure (avoids overwriting an existing figure) + + figure + +%% +% Plot the original data points with blue dots + + plot(t, y, 'bo', 'MarkerFaceColor', 'b', 'MarkerSize', 7) + +% Plot the curve fitting in the same figure using a solid red line + + hold on + plot(tt, f, 'r-', 'LineWidth', 1.5) + +% Set x-axis limit, label axes, and increase font size of the figure + + xlim([t(1), t(end)]) + xlabel("time (s)",'interpreter','latex') + ylabel("blood pressure (mmHg)",'interpreter','latex') + set(gca, 'FontSize', 20) + +end \ No newline at end of file diff --git a/gaussSeidel.asv b/gaussSeidel.asv new file mode 100644 index 0000000..c3502ea --- /dev/null +++ b/gaussSeidel.asv @@ -0,0 +1,35 @@ +function [x,k] = gaussSeidel(A,b,x0,tol) +%% +% Gauss-Seidel iterative method to approximate the solution of a +% linear system Ax=b up to a user defined tolerance +% +% INPUT: +% A - n by n square, non-singular matrix +% b - n by 1 right hand side vector +% x0 - n by 1 vector containing that initial guess for the iteration +% tol - user set tolerance for the stopping condition in the iteration +% +% OUTPUT: +% x - n by 1 vector containing the iterative solution +% k - number of iterations +%% +% get the system size + n = length(A); + max_its = 250; + x = zeros(n, 1); + for k = 1:max_its + for i = 1:n + sigma = 0; + for j = 1:n + if j ~= i + sigma = sigma + (A(i,j)*x(j)); + x(i) = (b(i) - sigma)/a(i,i); + if norm(b - A*x, 2) << + +%% +% Gauss-Seidel iteration which overwrites the current approximate solution +% with the new approximate solution (pseudocode available in the lecture +% notes on page 46) + +% +end \ No newline at end of file diff --git a/gaussSeidel.m b/gaussSeidel.m new file mode 100644 index 0000000..bc8e15b --- /dev/null +++ b/gaussSeidel.m @@ -0,0 +1,39 @@ +function [x,k] = gaussSeidel(A,b,x0,tol) +%% +% Gauss-Seidel iterative method to approximate the solution of a +% linear system Ax=b up to a user defined tolerance +% +% INPUT: +% A - n by n square, non-singular matrix +% b - n by 1 right hand side vector +% x0 - n by 1 vector containing that initial guess for the iteration +% tol - user set tolerance for the stopping condition in the iteration +% +% OUTPUT: +% x - n by 1 vector containing the iterative solution +% k - number of iterations +%% +% get the system size + n = length(A); + max_its = 500; + k = 0; + x_k = x0; + for k_it = 1:max_its + x_k = x0; + for i = 1:n + sigma = 0; + for j = 1:n + if j ~= i + sigma = sigma + (A(i,j)*x0(j)); + end + end + x0(i) = (b(i) - sigma)/A(i,i); + end + if norm(x0 - x_k, 2) < tol * norm(x0, 2) + k = k_it; + break + end + end + x = x0; + k= k_it; +end \ No newline at end of file diff --git a/graph_nils.m b/graph_nils.m new file mode 100644 index 0000000..7ca4163 --- /dev/null +++ b/graph_nils.m @@ -0,0 +1,10 @@ +loglog(ns,ts,'.-r','MarkerSize',25,'LineWidth',1.5) % plot the data +xlabel('matrix size (n)') +ylabel('time (sec)') +hold on +set(gca,'FontSize',16); +loglog(ns,ts(end)*(ns/ns(end)).^2,'--k','LineWidth',1.5) % plot a reference line + +grid on +axis tight +legend('data','O(n^2)','location','southeast') \ No newline at end of file diff --git a/lufact.m b/lufact.m new file mode 100644 index 0000000..68dcdce --- /dev/null +++ b/lufact.m @@ -0,0 +1,42 @@ +function [L,U,P] = lufact(A) +%% +% Computes the LU factorization of the matrix A. The factorization is done +% in-place and then the L and U matrices are extracted at the end for +% output. The factorization is PA = LU +% +% INPUT: +% A - n by n square, non-singular matrix +% +% OUTPUT: +% L - lower triangular matrix with ones along the main diagonal +% U - upper triangular matrix +% P - permutation matrix for pivoting +%% +% get the size of the system + n = length(A); +%% +% initialize the pivoting matrix + P = eye(n); +%% +% Vectorized in-place LU factorization (with row pivoting) that keeps +% track of the total permutations by scrambleing the matrix P + for j = 1:n-1 +%% +% Find the index of the largest pivot element in the current column + +%% +% Swap the rows within the in-place array as well as the permutation matrix P + + +%% +% Perform the in-place elimination and save the new column of L + i = j+1:n; % indices for the "active" matrix portion + A(i,j) = A(i,j)/A(j,j); + A(i,i) = A(i,i) - A(i,j)*A(j,i); +% A, return + end +%% +% Extract L and U from the in-place form + U = triu(A); + L = eye(n) + tril(A,-1); +end diff --git a/matrix_vec_mult.m b/matrix_vec_mult.m new file mode 100644 index 0000000..b04fc9a --- /dev/null +++ b/matrix_vec_mult.m @@ -0,0 +1,12 @@ +ns = 500:500:8000; +ts = 0*ns; +for j = 1:length(ns) + A = randn(ns(j),ns(j)); % create a matrix with random number entries + x = randn(ns(j),1); + tic; % start the timer + for i = 1:15 % repeat fifteen times + A*x; + end + t = toc; % read the timer + ts(j) = t/15; % normalise the timings with the number of runs +end diff --git a/pressure_data.m b/pressure_data.m new file mode 100644 index 0000000..e4c8709 --- /dev/null +++ b/pressure_data.m @@ -0,0 +1,8 @@ +t = (linspace(0, 0.384, 49))'; +y = [60.5620 60.8120 61.8750 64.1250 67.7500 72.3120 77.0000... + 81.1880 84.5620 87.3750 89.7500 91.8120 93.5620 94.9380... + 95.8750 96.4380 96.6250 96.5000 96.0620 95.4380 94.4380... + 92.9380 90.6250 87.8120 85.3750 83.4380 81.9380 80.5000... + 79.1880 78.0000 76.8120 75.7500 74.6880 73.6880 72.6250... + 71.6250 70.6250 69.6880 68.7500 67.8750 67.0000 66.1250... + 65.2500 64.3120 63.3750 62.5620 61.8120 61.1880 60.7500]'; \ No newline at end of file diff --git a/thinQR.m b/thinQR.m new file mode 100644 index 0000000..b51ee07 --- /dev/null +++ b/thinQR.m @@ -0,0 +1,38 @@ +function [Q,R] = thinQR(A) +%% +% This computes the "thin" version of the QR factorization using modified +% Gram-Schmidt to find the orthogonal columns of the matrix Q. +% +% INPUT: +% A - m x n matrix where m >= n that is assumed to have linearly +% independent columns +% +% OUTPUT: +% Q - m x n orthogonal matrix such that Q'*Q = eye(n) +% R - n x n upper triangular matrix + +%% +% pull the size of the rectangular system and allocate memory + [m,n] = size(A); +% + v = zeros(m,n); + Q = zeros(m,n); + R = zeros(n,n); +%% +% copy the columns of the matrix A into vectors that we will use for the +% modified Gram-Schmidt. Uses extra memory but done for clarity. + for j = 1:n + v(:,j) = A(:,j); + end +%% +% now go through and determine the columns of Q and the upper +% triangular matrix R + for j = 1:n + R(j,j) = norm(v(:,j), 2); + Q(:,j) = v(:,j) / R(j,j); + for k = j+1:n + R(j,k) = Q(:,j)' * v(:,k); + v(:,k) = v(:,k) - R(j,k) * Q(:,j); + end + end +end \ No newline at end of file -- 2.30.2