diff --git a/ENUME/projectB/Code/deriv.m b/ENUME/projectB/Code/deriv.m new file mode 100644 index 00000000..3f932383 --- /dev/null +++ b/ENUME/projectB/Code/deriv.m @@ -0,0 +1,12 @@ +% 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 new file mode 100644 index 00000000..2b7856af --- /dev/null +++ b/ENUME/projectB/Code/polynomial.m @@ -0,0 +1,4 @@ +% 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/printGraph.m b/ENUME/projectB/Code/printGraph.m new file mode 100644 index 00000000..f3bc99ba --- /dev/null +++ b/ENUME/projectB/Code/printGraph.m @@ -0,0 +1,34 @@ +% graph the real roots of a function +function printGraph(taskFunction, algorithmName, algorithm, interval, rootBrackets, plotTitle) + + grid on; % Get y values lines + hold on; % Retain current plot when adding new plots + title([plotTitle, algorithmName]); + set(gca, 'XAxisLocation', 'origin'); % Set properties of current axis + x = interval(1):0.01:interval(2); + % x is a vector of values between left and right interval with every value being higher by 0.01 + y = arrayfun(taskFunction, x); % We sketch the function from task for each x + plot(x, y); + + % iterate over rootBrackets and add them to the plot + for rootBracket = rootBrackets + % find all zeros within the bracket using the given algorithm + % Get all steps from the algorithm we use + [~, steps] = algorithm(taskFunction, rootBracket(1), rootBracket(2), 1e-10); + + firstStepColor = [1 0 0]; % Red + otherStepsColor = [0 1 0]; % Green + % plot first steps + scatter(steps(1, 1), steps(2, 1), [], firstStepColor); + text(steps(1, 1), steps(2, 1), 'firstStep', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top'); % It makes text appear neatly + % plot other steps + scatter(steps(1, 2:end), steps(2, 2:end), [], otherStepsColor); + + + + % print root table + disp([plotTitle, ' (', algorithmName, ')']); + columns = {'step', 'root', 'value at root'}; + disp(table((1:size(steps, 2))', steps(1, :)', steps(2, :)', 'VariableNames', columns)); + end +end diff --git a/ENUME/projectB/Code/printcomplex.m b/ENUME/projectB/Code/printcomplex.m new file mode 100644 index 00000000..e3843aeb --- /dev/null +++ b/ENUME/projectB/Code/printcomplex.m @@ -0,0 +1,39 @@ +% graph the complex roots of a function +function printcomplex(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]); + print(['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/rootBracketing.m b/ENUME/projectB/Code/rootBracketing.m new file mode 100644 index 00000000..76e7262d --- /dev/null +++ b/ENUME/projectB/Code/rootBracketing.m @@ -0,0 +1,34 @@ +% find the root brackets of a function within the given range +function rootBrackets = rootBracketing(givenFunction, intervalLeft, intervalRight) + [a, b, rootBrackets, resolution] = initializeValues(intervalLeft, intervalRight); + rootBrackets = bracketingLoop(a, b, rootBrackets, intervalRight, resolution, givenFunction); +end + +function [a, b, rootBrackets, resolution] = initializeValues(intervalLeft, intervalRight) + % define search resolution + resolution = (intervalRight - intervalLeft) / 6; + % The higher the value of denominator the less iterations will it take + % to reach the roots, however in order to have nice graph showing those + % brackets I will choose relatively small denominator - I have choosen + % the smallest natural number that still generates brackets on a graph + + % start search at the start of the range + a = intervalLeft; + b = intervalLeft + resolution; + rootBrackets = double.empty(2, 0); % initialize empty vector of size 2 +end + +function rootBrackets = bracketingLoop(a, b, rootBrackets, intervalRight, resolution, givenFunction) + while b ~= intervalRight % if the bracket can't be expanded end loop + % if the function changes sign inside the interval that means that we passed through a root that means that a bracket has been found + if sign(givenFunction(a)) ~= sign(givenFunction(b)) + % save bracket + rootBrackets(:, size(rootBrackets, 2) + 1) = [a, b]; % Add the new bracket to existing ones + end + % check next bracket + a = b; + b = min(a + resolution, intervalRight); + % Once a + resolution > intervalRight, then we will know that we + % reached beyond the interval and we must stop + end +end diff --git a/ENUME/projectB/Code/rootbrac.m b/ENUME/projectB/Code/rootbrac.m new file mode 100644 index 00000000..d84513b8 --- /dev/null +++ b/ENUME/projectB/Code/rootbrac.m @@ -0,0 +1,26 @@ +% finds the root brackets of a function within the given range +function brackets = rootBracketing(givenFunction, intervalLeft, intervalRight) + % define search resolution + resolution = (intervalRight - intervalLeft) / 10; + + % start search at the start of the range + a = intervalLeft; + b = intervalLeft + resolution; + brackets = double.empty(2, 0); + + % keep moving the interval until the range is exceeded + while 1 + % if the function changes sign inside the interval, a bracket has been found + if sign(givenFunction(a)) ~= sign(givenFunction(b)) + % save bracket + brackets(:, size(brackets, 2) + 1) = [a, b]; + end + + % if the bracket can't be expanded, return + if b == intervalRight; return; end + + % check next bracket + a = b; + b = min(a + resolution, intervalRight); + end +end diff --git a/ENUME/projectB/Code/task1Bisection.m b/ENUME/projectB/Code/task1Bisection.m new file mode 100644 index 00000000..c1f3f074 --- /dev/null +++ b/ENUME/projectB/Code/task1Bisection.m @@ -0,0 +1,75 @@ +interval = [-5, 10]; +rootBrackets = rootBracketing(@taskFunction, interval(1), interval(2)); + +printGraph(@taskFunction, 'False Position', @falsePosition, interval, rootBrackets, 'Approximate zeros of function for method of '); + +function y = taskFunction(x) + y = -2.1 + 0.3*x - x*exp(1)^(-x); +end + +function [zero, iterations] = falsePosition(taskFunction, a, b, tolerance) + + [iterations, lastTwoA, lastTwoB, i] = initialize(); + [zero, iterations, a, b, lastTwoA, lastTwoB] = firstTwoIterations(a, b, taskFunction, iterations, lastTwoA, lastTwoB); + [zero, iterations] = falsePositionLoop(taskFunction, zero, tolerance, lastTwoA, lastTwoB, i, a, b, iterations); + +end + +function [iterations, lastTwoA, lastTwoB, i] = initialize() + iterations = double.empty(2, 0); + lastTwoA = double.empty(2, 0); + lastTwoB = double.empty(2, 0); + i = 0; +end + +function [zero, iterations, a, b, lastTwoA, lastTwoB] = firstTwoIterations(a, b, taskFunction, iterations, lastTwoA, lastTwoB) + for j = 1 : 2 + zero = (a*taskFunction(b) - b * taskFunction(a)) / (taskFunction(b) - taskFunction(a)); + iterations(:, size(iterations, 2) + 1) = [zero, taskFunction(zero)]; + if sign(taskFunction(a)) ~= sign(taskFunction(zero)) + b = zero; + else + a = zero; + end + lastTwoA(j) = a; + lastTwoB(j) = b; + end +end + +function [zero, iterations] = falsePositionLoop(taskFunction, zero, tolerance, lastTwoA, lastTwoB, i, a, b, iterations) + while abs(taskFunction(zero)) > tolerance + [lastTwoA, i, a, lastTwoB, b, tolerance, zero, iterations] = insideLoop(lastTwoA, i, a, lastTwoB, b, tolerance, taskFunction, iterations); + end +end + +function [lastTwoA, i, a, lastTwoB, b, tolerance, zero, iterations] = insideLoop(lastTwoA, i, a, lastTwoB, b, tolerance, taskFunction, iterations) + [lastTwoA, lastTwoB] = changeLastTwoAB(lastTwoA, lastTwoB, i, a, b); + zero = calculateZero(lastTwoB, tolerance, a, b, lastTwoA, taskFunction); + iterations(:, size(iterations, 2) + 1) = [zero, taskFunction(zero)]; + [a, b] = newSubInterval(taskFunction, a, b, zero); + +end + +function [lastTwoA, lastTwoB] = changeLastTwoAB(lastTwoA, lastTwoB, i, a, b) + lastTwoA(mod(i, 2) + 1) = a; + lastTwoB(mod(i, 2) + 1) = b; +end + +function [zero] = calculateZero(lastTwoB, tolerance, a, b, lastTwoA, taskFunction) + if(abs(lastTwoB(1) - lastTwoB(2)) < tolerance) + zero = (a*(taskFunction(b) / 2) - b * taskFunction(a)) / (taskFunction(b) / 2 - taskFunction(a)); + elseif (abs(lastTwoA(1) - lastTwoA(2)) < tolerance) + zero = (a*taskFunction(b) - b * (taskFunction(a) / 2)) / (taskFunction(b) - (taskFunction(a) / 2)); + else + zero = (a*taskFunction(b) - b * taskFunction(a)) / (taskFunction(b) - taskFunction(a)); + end +end + +function [a, b] = newSubInterval(taskFunction, a, b, zero) + if sign(taskFunction(a)) ~= sign(taskFunction(zero)) + b = zero; + else + a = zero; + end +end + diff --git a/ENUME/projectB/Code/task1Newton.m b/ENUME/projectB/Code/task1Newton.m new file mode 100644 index 00000000..472be1cd --- /dev/null +++ b/ENUME/projectB/Code/task1Newton.m @@ -0,0 +1,43 @@ +interval = [-5, 10]; +rootBrackets = rootBracketing(@taskFunction, interval(1), interval(2)); + +printGraph(@taskFunction, 'Newton', @newtonMethod, interval, rootBrackets, 'Approximate zeros of function for method of '); + +function y = taskFunction(x) + y = -2.1 + 0.3*x - x*exp(1)^(-x); +end + +function [zero, iterations] = newtonMethod(taskFunction, a, b, tolerance) + [iterations, iteration, zero] = initialize(a, b); + [zero, iterations] = newtonLoop(iterations, iteration, zero, a, b, tolerance, taskFunction); +end + +function [iterations, step, zero] = initialize(a, b) + iterations = double.empty(2, 0); + step = sqrt(eps); + zero = (a + b) / 2; + iterations(:, size(iterations, 2) + 1) = [zero, taskFunction(zero)]; +end + +function [zero, iterations] = newtonLoop(iterations, iteration, zero, a, b, tolerance, taskFunction) + while abs(taskFunction(zero)) > tolerance + [zero, iterations] = insideLoop(taskFunction, zero, iteration, iterations, a, b); + end +end + +function [zero, iterations] = insideLoop(taskFunction, zero, iteration, iterations, a, b) + [zero, iterations] = calculateZeroIterations(taskFunction, zero, iteration, iterations); + checkForDivergence(zero, a, b); +end + +function [zero, iterations] = calculateZeroIterations(taskFunction, zero, iteration, iterations) + derivative = (taskFunction(zero + iteration) - taskFunction(zero - iteration)) / (2 * iteration); + zero = zero - taskFunction(zero) / derivative; + iterations(:, size(iterations, 2) + 1) = [zero, taskFunction(zero)]; +end + +function checkForDivergence(zero, a, b) + if zero < a || zero > b + error('Divergent iteration'); + end +end diff --git a/ENUME/projectB/report.aux b/ENUME/projectB/Report/report.aux similarity index 100% rename from ENUME/projectB/report.aux rename to ENUME/projectB/Report/report.aux diff --git a/ENUME/projectB/report.fdb_latexmk b/ENUME/projectB/Report/report.fdb_latexmk similarity index 100% rename from ENUME/projectB/report.fdb_latexmk rename to ENUME/projectB/Report/report.fdb_latexmk diff --git a/ENUME/projectB/report.fls b/ENUME/projectB/Report/report.fls similarity index 100% rename from ENUME/projectB/report.fls rename to ENUME/projectB/Report/report.fls diff --git a/ENUME/projectB/report.log b/ENUME/projectB/Report/report.log similarity index 100% rename from ENUME/projectB/report.log rename to ENUME/projectB/Report/report.log diff --git a/ENUME/projectB/report.out b/ENUME/projectB/Report/report.out similarity index 100% rename from ENUME/projectB/report.out rename to ENUME/projectB/Report/report.out diff --git a/ENUME/projectB/report.pdf b/ENUME/projectB/Report/report.pdf similarity index 100% rename from ENUME/projectB/report.pdf rename to ENUME/projectB/Report/report.pdf diff --git a/ENUME/projectB/report.synctex.gz b/ENUME/projectB/Report/report.synctex.gz similarity index 100% rename from ENUME/projectB/report.synctex.gz rename to ENUME/projectB/Report/report.synctex.gz diff --git a/ENUME/projectB/report.tex b/ENUME/projectB/Report/report.tex similarity index 100% rename from ENUME/projectB/report.tex rename to ENUME/projectB/Report/report.tex diff --git a/ENUME/projectB/report.toc b/ENUME/projectB/Report/report.toc similarity index 100% rename from ENUME/projectB/report.toc rename to ENUME/projectB/Report/report.toc diff --git a/ENUME/projectB/task-USOS_1153296.pdf b/ENUME/projectB/Task/task-USOS_1153296.pdf similarity index 100% rename from ENUME/projectB/task-USOS_1153296.pdf rename to ENUME/projectB/Task/task-USOS_1153296.pdf