Add task 3 code

This commit is contained in:
PolishPigeon 2021-12-03 04:11:21 +01:00
parent 207d97d4cc
commit 828a019165
4 changed files with 144 additions and 180 deletions

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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