Add Task 1 code

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

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

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

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

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

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

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
+

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