--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
projects / TANA21.git / commitdiff