From 207d97d4cc83ee0f785d06ee9b83b14c348d0264 Mon Sep 17 00:00:00 2001 From: PolishPigeon Date: Fri, 3 Dec 2021 03:28:09 +0100 Subject: [PATCH] Adding task 2 M2 --- ENUME/projectB/Code/deriv.m | 12 --- ENUME/projectB/Code/polynomial.m | 4 - ENUME/projectB/Code/printcomplex.m | 39 ---------- ENUME/projectB/Code/task2MM2.asv | 70 ++++++++++++++++++ ENUME/projectB/Code/task2MM2.m | 87 ++++++++++++++-------- ENUME/projectB/Code/task2muller.m | 114 ----------------------------- 6 files changed, 125 insertions(+), 201 deletions(-) delete mode 100644 ENUME/projectB/Code/deriv.m delete mode 100644 ENUME/projectB/Code/polynomial.m delete mode 100644 ENUME/projectB/Code/printcomplex.m create mode 100644 ENUME/projectB/Code/task2MM2.asv delete mode 100644 ENUME/projectB/Code/task2muller.m diff --git a/ENUME/projectB/Code/deriv.m b/ENUME/projectB/Code/deriv.m deleted file mode 100644 index 3f932383..00000000 --- a/ENUME/projectB/Code/deriv.m +++ /dev/null @@ -1,12 +0,0 @@ -% calculate the nth derivative of func at x -function y = deriv(func, x, deg) - % base case: zeroth derivative - if deg == 0 - y = func(x); - return - end - - % recurse to find the nth derivative - step = sqrt(eps); - y = (deriv(func, x + step, deg - 1) - deriv(func, x - step, deg - 1)) / (2 * step); -end diff --git a/ENUME/projectB/Code/polynomial.m b/ENUME/projectB/Code/polynomial.m deleted file mode 100644 index 2b7856af..00000000 --- a/ENUME/projectB/Code/polynomial.m +++ /dev/null @@ -1,4 +0,0 @@ -% the polynomial function for task 2 -function y = polynomial(x) - y = -2 * x^4 + 12 * x^3 + 4* x^2 + 1 * x + 3; -end diff --git a/ENUME/projectB/Code/printcomplex.m b/ENUME/projectB/Code/printcomplex.m deleted file mode 100644 index c7155e18..00000000 --- a/ENUME/projectB/Code/printcomplex.m +++ /dev/null @@ -1,39 +0,0 @@ -% graph the complex roots of a function -function printComplexGraph(func, algorithms, bracket, plottitle, outputsuffix) - % perform task for all available algorithms - for alg = 1:size(algorithms, 1) - [algname, algfunc] = algorithms{alg, :}; - - % prepare plot - figure; - grid on; - hold on; - title([plottitle, ' (', algname, ')']); - xlabel("Real part"); - ylabel("Imaginary part"); - set(gca, 'XAxisLocation', 'origin'); - - % find all zeros within the bracket using the given algorithm - [zero, steps] = algfunc(func, bracket(1), bracket(2), 1e-15); - - % plot steps on graph - plot(real(steps(1, :)), imag(steps(1, :)), '-o'); - text(real(steps(1, 1)), imag(steps(1, 1)), 'start', ... - 'HorizontalAlignment', 'center', ... - 'VerticalAlignment', 'top'); - text(real(steps(1, end)), imag(steps(1, end)), 'end', ... - 'HorizontalAlignment', 'center', ... - 'VerticalAlignment', 'top'); - - % finish and print graph - hold off; - set(gcf, 'PaperPosition', [0 0 6 4]); - set(gcf, 'PaperSize', [6 4]); - printComplexGraph(['report/', func2str(algfunc), outputsuffix], '-dpdf'); - - % print root table - disp([plottitle, ' (', algname, ')']); - columns = {'step', 'root', 'abs value at root'}; - disp(table([1:size(steps, 2)]', steps(1, :)', abs(steps(2, :))', 'VariableNames', columns)); - end -end diff --git a/ENUME/projectB/Code/task2MM2.asv b/ENUME/projectB/Code/task2MM2.asv new file mode 100644 index 00000000..67c6c69d --- /dev/null +++ b/ENUME/projectB/Code/task2MM2.asv @@ -0,0 +1,70 @@ +interval = [-5, 10]; +rootBrackets = rootBracketing(@polynomial, interval(1), interval(2)); + +printGraph(@polynomial, 'MM2', @mm2, interval, rootBrackets, 'Approximate zeros of function for method of '); + +printComplexGraph(@polynomial, 'MM2', @mm2, [-1 + i, 0], 'Aproximate complex roots of polynomial'); + +% the polynomial function for task 2 +function y = polynomial(x) + y = -2 * x^4 + 12 * x^3 + 4* x^2 + 1 * x + 3; +end + +% find roots of polynomial using MM2 +function [approximation, iterations] = mm2(polynomial, a, b, tolerance) + [approximation, iterations] = initialize(a, b, polynomial); + [approximation, iterations] = mm2Loop(approximation, iterations, polynomial, tolerance); +end + +function [approximation, iterations] = initialize(a, b, polynomial) + approximation = (a + b) / 2; + iterations = [approximation; polynomial(approximation)]; +end + +function [approximation, iterations] = mm2Loop(approximation, iterations, polynomial, tolerance) + while abs(polynomial(approximation)) > tolerance + [approximation, iterations] = insideLoop(approximation, polynomial, iterations); + end +end + +function [approximation, iterations] = insideLoop(approximation, polynomial, iterations) + [a, b, c] = getABC(approximation, polynomial); + [zPlus, zMinus] = findRoots(a, b, c); + newApproximation = chooseNewApproximation(zPlus, zMinus, approximation); + + approximation = newApproximation; + iterations(:, size(iterations, 2) + 1) = [approximation, polynomial(approximation)]; +end + +function [a, b, c] = getABC(approximation, polynomial) + c = polynomial(approximation); + b = derivative(polynomial, approximation, 1); + a = derivative(polynomial, approximation, 2) / 2; +end + +function [zPlus, zMinus] = findRoots(a, b, c) + zPlus = -2 * c / (b + sqrt(b ^ 2 - 4 * a * c)); + zMinus = -2 * c / (b - sqrt(b ^ 2 - 4 * a * c)); +end + +function newApproximation = chooseNewApproximation(zPlus, zMinus, approximation) + if abs(zPlus) < abs(zMinus) + newApproximation = approximation + zPlus; + else + newApproximation = approximation + zMinus; + end +end + +function [approximation, iterations] = updateApproximations(newApproximation, iterations) + + +% calculate the nth derivative of func at x +function y = derivative(function_, x, degree) + if degree == 0 + y = function_(x); + return + end + + step = sqrt(eps); + y = (derivative(function_, x + step, degree - 1) - derivative(function_, x - step, degree - 1)) / (2 * step); +end diff --git a/ENUME/projectB/Code/task2MM2.m b/ENUME/projectB/Code/task2MM2.m index 449c9ad9..1a0b1038 100644 --- a/ENUME/projectB/Code/task2MM2.m +++ b/ENUME/projectB/Code/task2MM2.m @@ -3,7 +3,7 @@ rootBrackets = rootBracketing(@polynomial, interval(1), interval(2)); printGraph(@polynomial, 'MM2', @mm2, interval, rootBrackets, 'Approximate zeros of function for method of '); -printComplexGraph(@polynomial, 'MM2', @mm2, [-1+i, 0], 'Aproximate complex roots of polynomial'); +printComplexGraph(@polynomial, 'MM2', @mm2, [-1 + i, 0], 'Aproximate complex roots of polynomial'); % the polynomial function for task 2 function y = polynomial(x) @@ -11,37 +11,60 @@ function y = polynomial(x) end % find roots of polynomial using MM2 -function [approx, steps] = mm2(func, a, b, tolerance) - % define current and (dummy) previous approximation point - approx = (a + b) / 2; - - % initialize output - steps = [approx; func(approx)]; - - % iterate algorithm until the error is within tolerance - % the error is defined as the diff between the prev and the current approx +function [approximation, iterations] = mm2(polynomial, a, b, tolerance) + [approximation, iterations] = initialize(a, b, polynomial); + [approximation, iterations] = mm2Loop(approximation, iterations, polynomial, tolerance); +end - func(approx) - while abs(func(approx)) > tolerance - % calculate approximating parabola using first and second derivative - c = func(approx); - b = deriv(func, approx, 1); - a = deriv(func, approx, 2) / 2; - - % find roots of parabola - zplus = -2 * c / (b + sqrt(b ^ 2 - 4 * a * c)); - zminus = -2 * c / (b - sqrt(b ^ 2 - 4 * a * c)); - - % choose root closer to current approximation - if abs(zplus) < abs(zminus) - newapprox = approx + zplus; - else - newapprox = approx + zminus; - end - - % update answer and prev approx - prevapprox = approx; - approx = newapprox; - steps(:, size(steps, 2) + 1) = [approx, func(approx)]; +function [approximation, iterations] = initialize(a, b, polynomial) + approximation = (a + b) / 2; + iterations = [approximation; polynomial(approximation)]; +end + +function [approximation, iterations] = mm2Loop(approximation, iterations, polynomial, tolerance) + while abs(polynomial(approximation)) > tolerance + [approximation, iterations] = insideLoop(approximation, polynomial, iterations); end end + +function [approximation, iterations] = insideLoop(approximation, polynomial, iterations) + [a, b, c] = getABC(approximation, polynomial); + [zPlus, zMinus] = findRoots(a, b, c); + newApproximation = chooseNewApproximation(zPlus, zMinus, approximation); + [approximation, iterations] = updateApproximations(newApproximation, iterations, polynomial); +end + +function [a, b, c] = getABC(approximation, polynomial) + c = polynomial(approximation); + b = derivative(polynomial, approximation, 1); + a = derivative(polynomial, approximation, 2) / 2; +end + +function [zPlus, zMinus] = findRoots(a, b, c) + zPlus = -2 * c / (b + sqrt(b ^ 2 - 4 * a * c)); + zMinus = -2 * c / (b - sqrt(b ^ 2 - 4 * a * c)); +end + +function newApproximation = chooseNewApproximation(zPlus, zMinus, approximation) + if abs(zPlus) < abs(zMinus) + newApproximation = approximation + zPlus; + else + newApproximation = approximation + zMinus; + end +end + +function [approximation, iterations] = updateApproximations(newApproximation, iterations, polynomial) + approximation = newApproximation; + iterations(:, size(iterations, 2) + 1) = [approximation, polynomial(approximation)]; +end + +% calculate the nth derivative of func at x +function y = derivative(function_, x, degree) + if degree == 0 + y = function_(x); + return + end + + step = sqrt(eps); + y = (derivative(function_, x + step, degree - 1) - derivative(function_, x - step, degree - 1)) / (2 * step); +end diff --git a/ENUME/projectB/Code/task2muller.m b/ENUME/projectB/Code/task2muller.m deleted file mode 100644 index c72bfc1b..00000000 --- a/ENUME/projectB/Code/task2muller.m +++ /dev/null @@ -1,114 +0,0 @@ -% define available algorithms -algorithms = { - 'MM1', @mm1; - 'MM2', @mm2 -}; - -% find all real root brackets -interval = [1, 7]; -brackets = rootbrac(@polynomial, interval(1), interval(2)); - -% find and graph real roots using both algorithms -printroots(@polynomial, algorithms, interval, brackets, ... - 'Approximate real roots of polynomial', 'realroots'); -printcomplex(@polynomial, algorithms, [-1+i, 0], ... - 'Approximate complex roots of polynomial', 'complexroots'); - -% find roots of polynomial using MM1 -function [zero, steps] = mm1(func, a, b, tolerance) - % define the three approximation points - apprx = [a, b, (a + b) / 2]; - apprxval = arrayfun(func, apprx); - - % initialize output - steps = [apprx(3); func(apprx(3))]; - - % iterate algorithm until the error is within tolerance - while abs(apprx(3) - apprx(2)) > tolerance - % prepare linear equation system to find parabola - z0 = apprx(1) - apprx(3); - z1 = apprx(2) - apprx(3); - diff0 = apprxval(1) - apprxval(3); - diff1 = apprxval(2) - apprxval(3); - - % solve equation system using Gaussian elimination (spaghetti code but fast) - eqsys = [z0 ^ 2, z0, diff0; z1 ^ 2, z1, diff1]; - reductor = eqsys(2, 1) / eqsys(1, 1); - eqsys(2, :) = eqsys(2, :) - reductor * eqsys(1, :); - eqsys(2, 1) = 0; - eqsys(2, :) = eqsys(2, :) ./ eqsys(2, 2); - eqsys(1, :) = eqsys(1, :) - eqsys(1, 2) * eqsys(2, :); - eqsys(1, :) = eqsys(1, :) ./ eqsys(1, 1); - - % define approximation parabola - a = eqsys(1, 3); - b = eqsys(2, 3); - c = apprxval(3); - - % find roots of parabola - zplus = -2 * c / (b + sqrt(b ^ 2 - 4 * a * c)); - zminus = -2 * c / (b - sqrt(b ^ 2 - 4 * a * c)); - - % choose root closer to current approximation - if abs(zplus) < abs(zminus) - newapprx = apprx(3) + zplus; - else - newapprx = apprx(3) + zminus; - end - - % update answer - zero = newapprx; - steps(:, size(steps, 2) + 1) = [zero, func(zero)]; - - % eliminate the most distant of the three approximations - worstapprxindex = -1; - worstapprxdiff = 0; - for i = 1:size(apprx, 2) - diff = abs(apprx(i) - newapprx); - if diff > worstapprxdiff - worstapprxindex = i; - worstapprxdiff = diff; - end - end - - % delete old approximation and append new one - apprx(worstapprxindex) = []; - apprx(3) = newapprx; - apprxval = arrayfun(func, apprx); - end -end - -% find roots of polynomial using MM2 -function [approx, steps] = mm2(func, a, b, tolerance) - % define current and (dummy) previous approximation point - approx = (a + b) / 2; - prevapprox = approx + b - a; - - % initialize output - steps = [approx; func(approx)]; - - % iterate algorithm until the error is within tolerance - % the error is defined as the diff between the prev and the current approx - while abs(approx - prevapprox) > tolerance - % calculate approximating parabola using first and second derivative - c = func(approx); - b = deriv(func, approx, 1); - a = deriv(func, approx, 2) / 2; - - % find roots of parabola - zplus = -2 * c / (b + sqrt(b ^ 2 - 4 * a * c)); - zminus = -2 * c / (b - sqrt(b ^ 2 - 4 * a * c)); - - % choose root closer to current approximation - if abs(zplus) < abs(zminus) - newapprox = approx + zplus; - else - newapprox = approx + zminus; - end - - % update answer and prev approx - prevapprox = approx; - approx = newapprox; - steps(:, size(steps, 2) + 1) = [approx, func(approx)]; - end -end