diff --git a/ENUME/projectB/Code/task2MM1.asv b/ENUME/projectB/Code/task2MM1.asv deleted file mode 100644 index c7b6a7fa..00000000 --- a/ENUME/projectB/Code/task2MM1.asv +++ /dev/null @@ -1,110 +0,0 @@ -interval = [-5, 10]; -rootBrackets = rootBracketing(@polynomial, interval(1), interval(2)); - -printGraph(@polynomial, 'MM1', @mm1, interval, rootBrackets, 'Approximate zeros of function for method of '); - -printComplexGraph(@polynomial, 'MM1', @mm1, [-1 + i, 0], 'Aproximate complex roots of polynomial'); - -function y = polynomial(x) - y = -2 * x^4 + 12 * x^3 + 4* x^2 + 1 * x + 3; -end - -function [approximation, iterations] = mm1(polynomial, a, b, tolerance) - - [approximation, approximationValue, iterations] = initialize(a, b, polynomial); - [approximation, iterations] = mm1Loop(approximation, tolerance, approximationValue, iterations, polynomial); -end - -function [approximation, approximationValue, iterations] = initialize(a, b, polynomial) - approximation = [a, b, (a + b) / 2]; - approximationValue = arrayfun(polynomial, approximation); - iterations = [approximation(3); polynomial(approximation(3))]; -end - -function [approximation, iterations] = mm1Loop(approximation, tolerance, approximationValue, iterations, polynomial) - while abs(polynomial(approximation(3))) > tolerance - [approximation, approximationValue, iterations] = insideLoop(approximation, approximationValue, polynomial, iterations); - end -end - -function [approximation, approximationValue, iterations] = insideLoop(approximation, approximationValue, polynomial, iterations) - eqsys = createEquationSystem(approximation, approximationValue); - [zplus, zminus] = rootsOfQuadraticFormula(eqsys, approximationValue); - [approximation, approximationValue, iterations] = updateApproximations(zplus, zminus, approximation, iterations, polynomial); -end - -function eqsys = createEquationSystem(approximation, approximationValue) - [z0, z1, diff0, diff1] = initializeEquationSystem(approximation, approximationValue); - eqsys = solveEquationSystem(z0, diff0, z1, diff1); -end - -function [zplus, zminus] = rootsOfQuadraticFormula(eqsys, approximationValue) - [a, b, c] = createApproximatedQuadraticFormula(eqsys, approximationValue); - [zplus, zminus] = findRootsOfQuadraticFormula(a, b, c); -end - -function [approximation, approximationValue, iterations] = updateApproximations(zplus, zminus, approximation, iterations, polynomial) - newapprx = chooseNewRoot(zplus, zminus, approximation); - iterations = addZeroToIterationVector(newapprx, iterations, polynomial); - worstapprxindex = getWorstApproximationIndex(approximation, newapprx); - [approximation, approximationValue] = deleteWorstApproximation(worstapprxindex, approximation, polynomial, newapprx); -end - -function [z0, z1, diff0, diff1] = initializeEquationSystem(approximation, approximationValue) - z0 = approximation(1) - approximation(3); - z1 = approximation(2) - approximation(3); - diff0 = approximationValue(1) - approximationValue(3); - diff1 = approximationValue(2) - approximationValue(3); -end - -function eqsys = solveEquationSystem(z0, diff0, z1, diff1) - 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); -end - -function [a, b, c] = createApproximatedQuadraticFormula(eqsys, approximationValue) - a = eqsys(1, 3); - b = eqsys(2, 3); - c = approximationValue(3); -end - -function [zplus, zminus] = findRootsOfQuadraticFormula(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 newapprx = chooseNewRoot(zplus, zminus, approximation) - if abs(zplus) < abs(zminus) - newapprx = approximation(3) + zplus; - else - newapprx = approximation(3) + zminus; - end -end - -function iterations = addZeroToIterationVector(newapprx, iterations, polynomial) - zero = newapprx; - iterations(:, size(iterations, 2) + 1) = [zero, polynomial(zero)]; -end - -function worstapprxindex = getWorstApproximationIndex(approximation, newapprx) - worstapprxindex = -1; - worstapprxdiff = 0; - for i = 1:size(approximation, 2) - diff = abs(approximation(i) - newapprx); - if diff > worstapprxdiff - worstapprxindex = i; - worstapprxdiff = diff; - end - end -end - -function [approximation, approximationValue] = deleteWorstApproximation(worstapprxindex, approximation, polynomial, newapprx) - approximation(worstapprxindex) = []; - approximation(3) = newapprx; - approximationValue = arrayfun(polynomial, approximation); -end \ No newline at end of file diff --git a/ENUME/projectB/Code/task2MM2.asv b/ENUME/projectB/Code/task2MM2.asv deleted file mode 100644 index 67c6c69d..00000000 --- a/ENUME/projectB/Code/task2MM2.asv +++ /dev/null @@ -1,70 +0,0 @@ -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/task3.asv b/ENUME/projectB/Code/task3.asv new file mode 100644 index 00000000..c4ccf633 --- /dev/null +++ b/ENUME/projectB/Code/task3.asv @@ -0,0 +1,72 @@ +interval = [-5, 10]; +rootBrackets = rootBracketing(@polynomial, interval(1), interval(2)); + +printGraph(@polynomial, 'Laguerre', @laguerre, interval, rootBrackets, 'Approximate zeros of function for method of '); + +printComplexGraph(@polynomial, 'Laguerre', @laguerre, [-1 + i, 0], 'Aproximate complex roots of polynomial'); + +function y = polynomial(x) + y = -2 * x^4 + 12 * x^3 + 4* x^2 + 1 * x + 3; +end + +function [zero, iterations] = laguerre(polynomial, a, b, tolerance) + [degree, zero, iterations] = initialize(a, b, polynomial); + [zero, iterations] = laguerreLoop(polynomial, zero, tolerance, iterations, degree); + +end + +function [degree, zero, iterations] = initialize(a, b, polynomial) + degree = 4; + zero = (a + b) / 2; + iterations = [zero; polynomial(zero)]; +end + +function [zero, iterations] = laguerreLoop(polynomial, zero, tolerance, iterations, degree) + while abs(polynomial(zero)) > tolerance + [iterations, zero] = insideLoop(polynomial, zero, degree, iterations); + end +end + +function [iterations, zero] = insideLoop(polynomial, zero, degree, iterations) + [derrivative0, derrivative1, derrivative2] = calculateDerrivatives(polynomial, zero); + [zPlus, zMinus] = calculateZ(degree, derrivative0, derrivative1, derrivative2); + + + % update answer + zero = newzero; + iterations(:, size(iterations, 2) + 1) = [zero, polynomial(zero)]; +end + +function [derrivative0, derrivative1, derrivative2] = calculateDerrivatives(polynomial, zero) + % calculate derivatives + derrivative0 = polynomial(zero); + derrivative1 = derivative(polynomial, zero, 1); + derrivative2 = derivative(polynomial, zero, 2); +end + +function [zPlus, zMinus] = calculateZ(degree, derrivative0, derrivative1, derrivative2) + expressionUnderSquareRoot = (degree - 1) * ((degree - 1) * derrivative1 ^ 2 - degree * derrivative0 * derrivative2); + lagsqrt = sqrt(expressionUnderSquareRoot); + + zPlus = degree * derrivative0 / (derrivative1 + lagsqrt); + zMinus = degree * derrivative0 / (derrivative1 - lagsqrt); +end + +function [] = chooseNewZero() + % choose value closer to current approximation + if abs(zPlus) < abs(zMinus) + newzero = zero - zPlus; + else + newzero = zero - zMinus; + end + +% calculate the nth derivative of polynomial 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/task3.m b/ENUME/projectB/Code/task3.m new file mode 100644 index 00000000..d837642a --- /dev/null +++ b/ENUME/projectB/Code/task3.m @@ -0,0 +1,72 @@ +interval = [-5, 10]; +rootBrackets = rootBracketing(@polynomial, interval(1), interval(2)); + +printGraph(@polynomial, 'Laguerre', @laguerre, interval, rootBrackets, 'Approximate zeros of function for method of '); + +printComplexGraph(@polynomial, 'Laguerre', @laguerre, [-1 + i, 0], 'Aproximate complex roots of polynomial'); + +function y = polynomial(x) + y = -2 * x^4 + 12 * x^3 + 4* x^2 + 1 * x + 3; +end + +function [zero, iterations] = laguerre(polynomial, a, b, tolerance) + [degree, zero, iterations] = initialize(a, b, polynomial); + [zero, iterations] = laguerreLoop(polynomial, zero, tolerance, iterations, degree); + +end + +function [degree, zero, iterations] = initialize(a, b, polynomial) + degree = 4; + zero = (a + b) / 2; + iterations = [zero; polynomial(zero)]; +end + +function [zero, iterations] = laguerreLoop(polynomial, zero, tolerance, iterations, degree) + while abs(polynomial(zero)) > tolerance + [iterations, zero] = insideLoop(polynomial, zero, degree, iterations); + end +end + +function [iterations, zero] = insideLoop(polynomial, zero, degree, iterations) + [derrivative0, derrivative1, derrivative2] = calculateDerrivatives(polynomial, zero); + [zPlus, zMinus] = calculateZ(degree, derrivative0, derrivative1, derrivative2); + newZero = chooseNewZero(zPlus, zMinus, zero); + [zero, iterations] = updateZeros(newZero, iterations, polynomial); +end + +function [derrivative0, derrivative1, derrivative2] = calculateDerrivatives(polynomial, zero) + derrivative0 = polynomial(zero); + derrivative1 = derivative(polynomial, zero, 1); + derrivative2 = derivative(polynomial, zero, 2); +end + +function [zPlus, zMinus] = calculateZ(degree, derrivative0, derrivative1, derrivative2) + expressionUnderSquareRoot = (degree - 1) * ((degree - 1) * derrivative1 ^ 2 - degree * derrivative0 * derrivative2); + lagsqrt = sqrt(expressionUnderSquareRoot); + + zPlus = degree * derrivative0 / (derrivative1 + lagsqrt); + zMinus = degree * derrivative0 / (derrivative1 - lagsqrt); +end + +function newZero = chooseNewZero(zPlus, zMinus, zero) + if abs(zPlus) < abs(zMinus) + newZero = zero - zPlus; + else + newZero = zero - zMinus; + end +end + +function [zero, iterations] = updateZeros(newZero, iterations, polynomial) + zero = newZero; + iterations(:, size(iterations, 2) + 1) = [zero, polynomial(zero)]; +end + +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