From c916b65bd9829499db672e5813968bf50716e565 Mon Sep 17 00:00:00 2001 From: PolishPigeon Date: Fri, 7 Jan 2022 22:44:18 +0100 Subject: [PATCH] Beautifying task1.m --- ENUME/projectC/functionDataPoints.m | 16 ++++++ ENUME/projectC/task1.m | 84 ++++++++++++++++++++++++++++ ENUME/projectC/task1lsapprox.m | 86 +++++++++++++---------------- 3 files changed, 137 insertions(+), 49 deletions(-) create mode 100644 ENUME/projectC/functionDataPoints.m create mode 100644 ENUME/projectC/task1.m diff --git a/ENUME/projectC/functionDataPoints.m b/ENUME/projectC/functionDataPoints.m new file mode 100644 index 00000000..755d28c2 --- /dev/null +++ b/ENUME/projectC/functionDataPoints.m @@ -0,0 +1,16 @@ +function dataPoints = functionDataPoints() + dataPoints = (-5:5)'; + dataPoints(:, 2) = [ + -6.5743 + 0.9765 + 3.1026 + 1.8572 + 1.3165 + -0.6144 + 0.1032 + 0.3729 + 2.5327 + 7.3857 + 9.4892 + ]; +end \ No newline at end of file diff --git a/ENUME/projectC/task1.m b/ENUME/projectC/task1.m new file mode 100644 index 00000000..8f2aca3a --- /dev/null +++ b/ENUME/projectC/task1.m @@ -0,0 +1,84 @@ +function task1(polynomialDegree) + for currentPolynomialDegree = 0 : polynomialDegree + dataPoints = functionDataPoints(); + % obtain LSPSolutions of approximating polynomial + [LSPSolutions, error, conditionNumber] = approximate(dataPoints, currentPolynomialDegree); + displayInfo(currentPolynomialDegree, error, conditionNumber); + plotGraph(currentPolynomialDegree, dataPoints, LSPSolutions) + + end +end + +function displayInfo(currentPolynomialDegree, error, conditionNumber) + disp("Current Polynomial Degree:"); + disp(currentPolynomialDegree); + disp("Error: "); + disp(error); + disp("Condition number of Gram's Matrix: "); + disp(conditionNumber); +end + +function plotGraph(currentPolynomialDegree, dataPoints, LSPSolutions) + figure; + grid on; + hold on; + plotDataPoints(currentPolynomialDegree, dataPoints); + plotApproximation(dataPoints, LSPSolutions); + hold off; +end + +function plotDataPoints(currentPolynomialDegree, dataPoints) + title(['Polynomial approximation of the funtion using: ', num2str(currentPolynomialDegree), ' degree']); + scatter(dataPoints(:, 1), dataPoints(:, 2), 72, 'x'); +end + +function plotApproximation(dataPoints, LSPSolutions) + x = dataPoints(1):0.05:dataPoints(end, 1); + y = valueApproximationAtx(LSPSolutions, x); + plot(x, y); +end + + +% find the approximating polynomial of the given degree +function [LSPSolutions, error, conditionNumber] = approximate(dataPoints, polydeg) + % initialize A matrix + A = zeros(size(dataPoints, 1), polydeg + 1); + + A = calculateACells(dataPoints, A); + LSPSolutions = solveLSP(A, dataPoints); + + % Calculate error of solution + error = norm(dataPoints(:, 2) - A * LSPSolutions); + % Calculate condition number of Gram's matrix + GramsMatrix = A' * A; + conditionNumber = cond(GramsMatrix); +end + +function Matrix = calculateACells(dataPoints, Matrix) + for row = 1 : size(Matrix, 1) + for column = 1 : size(Matrix, 2) + Matrix(row, column) = dataPoints(row, 1) ^ (column - 1); + end + end +end + +function solutions = solveLSP(Matrix, dataPoints) + [Q, eqsys, invqtq] = QRDecomposition(Matrix); + eqsys(:, end + 1) = invqtq * Q' * dataPoints(:, 2); + solutions = backSubstitution(eqsys); +end + +% evaluate the value of an approximation at the given x +function arrayOfValues = valueApproximationAtx(LSPSolutions, arrayOfArguments) + arrayOfValues = zeros(1, size(arrayOfArguments, 2)); + arrayOfValues = calculateArrayofValues(arrayOfArguments, arrayOfValues, LSPSolutions); + +end + +function arrayOfValues = calculateArrayofValues(arrayOfArguments, arrayOfValues, LSPSolutions) + for x = 1 : size(arrayOfArguments, 2) + for i = 1 : size(LSPSolutions, 1) + arrayOfValues(x) = arrayOfValues(x) + LSPSolutions(i) * arrayOfArguments(x) ^ (i - 1); + end + end +end diff --git a/ENUME/projectC/task1lsapprox.m b/ENUME/projectC/task1lsapprox.m index fe97601d..bcc2ee3d 100644 --- a/ENUME/projectC/task1lsapprox.m +++ b/ENUME/projectC/task1lsapprox.m @@ -1,68 +1,56 @@ -% define function data points -taskfunc = (-5:5)'; -taskfunc(:, 2) = [ - -6.5743 - 0.9765 - 3.1026 - 1.8572 - 1.3165 - -0.6144 - 0.1032 - 0.3729 - 2.5327 - 7.3857 - 9.4892 -]; - -% perform the task -for polydeg = 0:3 - % obtain factors of approximating polynomial - [factors, error, gramcond] = approximate(taskfunc, polydeg); - - % print error and condition number - disp(['Approximation degree ', num2str(polydeg), ':']); - disp(['Error: ', num2str(error)]); - disp(['Condition number: ', num2str(gramcond)]); - - % plot data points - figure; - grid on; - hold on; - title(['Polynomial approximation of the function (degree ', ... - num2str(polydeg), ')']); - scatter(taskfunc(:, 1), taskfunc(:, 2)); - - % plot approximation - x = taskfunc(1):0.05:taskfunc(end, 1); - y = evalapprox(factors, x); - plot(x, y); - - % finish and print graph - hold off; - set(gcf, 'PaperPosition', [0 0 6 4]); - set(gcf, 'PaperSize', [6 4]); - %print(['report/approx', num2str(polydeg)], '-dpdf'); +function task1() + for polydeg = 0:3 + dataPoints = functionDataPoints(); + % obtain factors of approximating polynomial + [factors, error, gramcond] = approximate(dataPoints, polydeg); + + % print error and condition number + disp(['Approximation degree ', num2str(polydeg), ':']); + disp(['Error: ', num2str(error)]); + disp(['Condition number: ', num2str(gramcond)]); + + % plot data points + figure; + grid on; + hold on; + title(['Polynomial approximation of the function (degree ', ... + num2str(polydeg), ')']); + scatter(dataPoints(:, 1), dataPoints(:, 2)); + + % plot approximation + x = dataPoints(1):0.05:dataPoints(end, 1); + y = evalapprox(factors, x); + plot(x, y); + + % finish and print graph + hold off; + set(gcf, 'PaperPosition', [0 0 6 4]); + set(gcf, 'PaperSize', [6 4]); + %print(['report/approx', num2str(polydeg)], '-dpdf'); + end end + + % find the approximating polynomial of the given degree -function [factors, error, gramcond] = approximate(func, polydeg) +function [factors, error, gramcond] = approximate(dataPoints, polydeg) % define the A matrix used for solving the error minimization problem - A = zeros(size(func, 1), polydeg + 1); + A = zeros(size(dataPoints, 1), polydeg + 1); % calculate cells of A using natural basis for row = 1:size(A, 1) for col = 1:size(A, 2) - A(row, col) = func(row, 1) ^ (col - 1); + A(row, col) = dataPoints(row, 1) ^ (col - 1); end end % solve least-square problem using QRdash decomposition [Q, eqsys, invqtq] = QRDecomposition(A); - eqsys(:, end + 1) = invqtq * Q' * func(:, 2); + eqsys(:, end + 1) = invqtq * Q' * dataPoints(:, 2); factors = backSubstitution(eqsys); % calculate error and condition number of Gram's matrix - error = norm(func(:, 2) - A * factors); + error = norm(dataPoints(:, 2) - A * factors); gramcond = cond(A' * A); end