kuhyx/WUT_Computer_Science
1
0
Fork 0
mirror of https://github.com/kuhyx/WUT_Computer_Science.git synced 2026-11-11 12:27:40 +01:00
Code Issues Actions Packages Projects Releases Wiki Activity

Adding task 2 M2

Browse Source
This commit is contained in:
Krzysztof Rudnicki 2021-12-03 03:28:09 +01:00
parent 55874207e7
commit 207d97d4cc
6 changed files with 125 additions and 201 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
ENUME/projectB/Code/deriv.m
View File

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

4
ENUME/projectB/Code/polynomial.m
View File

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

39
ENUME/projectB/Code/printcomplex.m
View File

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

70
ENUME/projectB/Code/task2MM2.asv Normal file
View File

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

87
ENUME/projectB/Code/task2MM2.m
View File

@ -3,7 +3,7 @@ rootBrackets = rootBracketing(@polynomial, interval(1), interval(2));
printGraph(@polynomial, 'MM2', @mm2, interval, rootBrackets, 'Approximate zeros of function for method of '); 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 % the polynomial function for task 2
function y = polynomial(x) function y = polynomial(x)
@ -11,37 +11,60 @@ function y = polynomial(x)
end end
% find roots of polynomial using MM2 % find roots of polynomial using MM2
function [approx, steps] = mm2(func, a, b, tolerance) function [approximation, iterations] = mm2(polynomial, a, b, tolerance)
% define current and (dummy) previous approximation point [approximation, iterations] = initialize(a, b, polynomial);
approx = (a + b) / 2; [approximation, iterations] = mm2Loop(approximation, iterations, polynomial, tolerance);
end
% 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
func(approx) function [approximation, iterations] = initialize(a, b, polynomial)
while abs(func(approx)) > tolerance approximation = (a + b) / 2;
% calculate approximating parabola using first and second derivative iterations = [approximation; polynomial(approximation)];
c = func(approx); end
b = deriv(func, approx, 1);
a = deriv(func, approx, 2) / 2; function [approximation, iterations] = mm2Loop(approximation, iterations, polynomial, tolerance)
while abs(polynomial(approximation)) > tolerance
% find roots of parabola [approximation, iterations] = insideLoop(approximation, polynomial, iterations);
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
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

114
ENUME/projectB/Code/task2muller.m
View File

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