feat: Adding code

ENUME/projectC/QRDecomposition.m

@@ -0,0 +1,26 @@
+function [Q, R, invqtq] = QRDecomposition(A)
+    % initialize empty matrices
+    Q = zeros(size(A));
+    R = eye(size(A, 2));
+    invqtq = zeros(size(A, 2));
+    
+    % modified Gram-Schmidt, use each column to orthogonalize the ones in front of it
+    for col = 1:size(A, 2)
+        % by the time we've reached this column, it's already been orthogonalized
+        Q(:, col) = A(:, col);
+        
+        % calculate current column dot product for R
+        coldot = dot(Q(:, col), Q(:, col));
+        invqtq(col, col) = 1 / coldot;
+        
+        % orthogonalize further columns
+        for next = (col + 1):size(A, 2)
+            % calculate R cell for this column pair
+            R(col, next) = dot(Q(:, col), A(:, next)) / coldot;
+            
+            % orthogonalize column
+            A(:, next) = A(:, next) - R(col, next) * Q(:, col);
+        end
+    end
+   
+end

ENUME/projectC/backSubstitution.m

@@ -0,0 +1,16 @@
+% solves system with triangular matrix
+function result = backSubstitution(eqsys)
+    for col = size(eqsys, 1):-1:1
+        % normalize diagonal coefficients to 1
+        eqsys(col, :) = eqsys(col, :) / eqsys(col, col);
+        
+        % eliminate factor from other rows
+        for row = (col - 1):-1:1
+            reductor = eqsys(row, col) / eqsys(col, col);
+            eqsys(row, :) = eqsys(row, :) - eqsys(col, :) * reductor;
+        end
+    end
+    
+    % rightmost column is now the result
+    result = eqsys(:, size(eqsys, 2));
+end

ENUME/projectC/rk4.m

@@ -0,0 +1,31 @@
+% solve ODE system using RK4 with constant step size
+function [x, derivs] = rk4(functs, init, interval, stepsize, maxsteps)
+    % set initial values as start points of output
+    x = init;
+    
+    % build derivatives table
+    derivs = zeros(size(x));
+    for eqnum = 1:size(functs, 1)
+        derivs(eqnum, 1) = functs{eqnum}(x(:, 1));
+    end
+    
+    % build output based on preceding values
+    stepcount = ceil((interval(2) - interval(1)) / stepsize);
+    if nargin >= 5; stepcount = min(stepcount, maxsteps - 1); end
+    for step = 1:stepcount
+        % obtain the preceding function values
+        stepval = x(:, step);
+        
+        for eqnum = 1:size(functs, 1)
+            % generic single-step iteration
+            phi = rk4phi(functs{eqnum}, stepval, stepsize);
+            x(eqnum, step + 1) = x(eqnum, step) + stepsize * phi;
+            
+            % update derivatives table
+            derivs(eqnum, step + 1) = functs{eqnum}(x(:, step + 1));
+        end
+    end
+    
+    % append arguments to output
+    x = [interval(1):stepsize:(stepcount * stepsize); x];
+end

ENUME/projectC/rk4auto.m

@@ -0,0 +1,61 @@
+% automatic step size variant of RK4
+function [x, sizes, errors] = rk4auto(functs, init, interval, initstep, eps_rel, eps_abs)
+    % set start points of output
+    args = interval(1);
+    x = init;
+    
+    % initialize output plots
+    sizes = double.empty();
+    errors = double.empty();
+    
+    % integrate function until end of interval reached
+    stepsize = initstep;
+    step = 0;
+    while 1
+        % obtain the preceding function values
+        step = step + 1;
+        stepval = x(:, step);
+        
+        % advance output function
+        for eqnum = 1:size(functs, 1)
+            % generic single-step iteration
+            phi = rk4phi(functs{eqnum}, stepval, stepsize);
+            x(eqnum, step + 1) = x(eqnum, step) + stepsize * phi;
+        end
+        
+        % stop algorithm if function integrated over the whole interval
+        args(step + 1) = args(step) + stepsize;
+        if args(end) >= interval(2); break; end
+        
+        % also calculate next step using two half-steps
+        for substep = 1:2
+            for eqnum = 1:size(functs, 1)
+                phi = rk4phi(functs{eqnum}, stepval, stepsize / 2);
+                stepval(eqnum) = stepval(eqnum) + (stepsize / 2) * phi;
+            end
+        end
+        
+        % calculate step correction factor
+        alpha = Inf;
+        for eqnum = 1:size(functs, 1)
+            % calculate approximation error
+            delta = abs(stepval(eqnum) - x(eqnum, step + 1)) / 15;
+            errors(step) = delta;
+            
+            % calculate equation-specific alpha
+            epsilon = abs(stepval(eqnum)) * eps_rel + eps_abs;
+            eqalpha = epsilon / delta;
+            
+            % minimum alpha wins
+            if eqalpha < alpha; alpha = eqalpha; end
+        end
+        alpha = alpha ^ (1/5);
+        
+        % correct step size with safety factor
+        stepsize = 0.9 * alpha * stepsize;
+        sizes(step) = stepsize;
+    end
+    
+    % append arguments to output
+    x = [args; x];
+end

ENUME/projectC/rk4phi.m

@@ -0,0 +1,8 @@
+% calculate the phi for RK4 algorithms
+function phi = rk4phi(fun, stepval, stepsize)
+    k1 = fun(stepval);
+    k2 = fun(stepval + 0.5 * stepsize * k1);
+    k3 = fun(stepval + 0.5 * stepsize * k2);
+    k4 = fun(stepval + stepsize * k3);
+    phi = (k1 + 2 * k2 + 2 * k3 + k4) / 6;
+end

ENUME/projectC/task1lsapprox.m

@@ -0,0 +1,78 @@
+% define function data points
+taskfunc = (-5:5)';
+taskfunc(:, 2) = [
+    -6.5743
+    0.9765
+    3.1026
+    1.8572
+    1.3165
+    -0.6144
+    0.1032
+    0.3729
+    2.5327
+    7.3857
+    9.4892
+];
+
+% perform the task
+for polydeg = 0:3
+    % obtain factors of approximating polynomial
+    [factors, error, gramcond] = approximate(taskfunc, polydeg);
+    
+    % print error and condition number
+    disp(['Approximation degree ', num2str(polydeg), ':']);
+    disp(['Error: ', num2str(error)]);
+    disp(['Condition number: ', num2str(gramcond)]);
+    
+    % plot data points
+    figure;
+    grid on;
+    hold on;
+    title(['Polynomial approximation of the function (degree ', ...
+        num2str(polydeg), ')']);
+    scatter(taskfunc(:, 1), taskfunc(:, 2));
+    
+    % plot approximation
+    x = taskfunc(1):0.05:taskfunc(end, 1);
+    y = evalapprox(factors, x);
+    plot(x, y);
+    
+    % finish and print graph
+    hold off;
+    set(gcf, 'PaperPosition', [0 0 6 4]);
+    set(gcf, 'PaperSize', [6 4]);
+    %print(['report/approx', num2str(polydeg)], '-dpdf');
+end
+
+% find the approximating polynomial of the given degree
+function [factors, error, gramcond] = approximate(func, polydeg)
+    % define the A matrix used for solving the error minimization problem
+    A = zeros(size(func, 1), polydeg + 1);
+    
+    % calculate cells of A using natural basis
+    for row = 1:size(A, 1)
+        for col = 1:size(A, 2)            
+            A(row, col) = func(row, 1) ^ (col - 1);
+        end
+    end
+    
+    % solve least-square problem using QRdash decomposition
+    [Q, eqsys, invqtq] = QRDecomposition(A);
+    eqsys(:, end + 1) = invqtq * Q' * func(:, 2);
+    factors = backSubstitution(eqsys);
+    
+    % calculate error and condition number of Gram's matrix
+    error = norm(func(:, 2) - A * factors);
+    gramcond = cond(A' * A);
+end
+
+% evaluate the value of an approximation at the given x
+function y = evalapprox(factors, xarray)
+    y = zeros(1, size(xarray, 2));
+    
+    for xi = 1:size(xarray, 2)
+        for i = 1:size(factors, 1)
+            y(xi) = y(xi) + factors(i) * xarray(xi) ^ (i - 1);
+        end
+    end
+end

ENUME/projectC/task2ode.asv

@@ -0,0 +1,104 @@
+% functions of the ODE system and initial values
+sysfuncts = {
+    @(x) x(2) + x(1) * (0.5 - x(1)^2 - x(2)^2);
+    @(x) -x(1) + x(2) * (0.5 - x(1)^2 - x(2)^2)
+};
+initvalues = [0; 12];
+interval = [0; 15];
+
+% define available algorithms
+algorithms = {
+    'RK4', @rk4, [0.01, 0.013408];
+    'Adams PC', @adamspc, [0.002, 0.01305]
+};
+
+% solve ODE using different algorithms and step sizes
+for alg = 1:size(algorithms, 1)
+    [algname, algfunc, stepsizes] = algorithms{alg, :};
+    
+    % solve using the given algorithm for each step size
+    stepresults = cell(size(stepsizes, 2), 3);
+    stepnames = {'optimal step', 'larger step'};
+    for stepno = 1:size(stepsizes, 2)
+        result = algfunc(sysfuncts, initvalues, interval, stepsizes(stepno));
+        stepresults(stepno, :) = { ...
+            stepsizes(stepno), ...
+            stepnames{stepno}, ...
+            result};
+    end
+    
+    % plot each component against time
+    for eqnum = 1:size(sysfuncts, 1)
+        % begin plot
+        figure; grid on; hold on;
+        title([algname, ', x_', num2str(eqnum), ' against time']);
+        set(gcf, 'PaperPosition', [0 0 6 4]);
+        set(gcf, 'PaperSize', [6 4]);
+        
+        % plot component for each step size
+        for stepresult = stepresults'
+            plot(stepresult{3}(1, :), stepresult{3}(eqnum + 1, :));
+        end
+        
+        % finish plot
+        hold off;
+        legend(stepresults{:, 2});
+%        %print(['report/', func2str(algfunc), 'x', num2str(eqnum)], '-dpdf');
+    end
+    
+    % plot first two components against each other
+    figure; grid on; hold on;
+    title([algname, ' trajectory plot (x_2 against x_1)']);
+    set(gcf, 'PaperPosition', [0 0 6 4]);
+    set(gcf, 'PaperSize', [6 4]);
+    
+    % plot for each step size
+    for stepresult = stepresults'
+        plot(stepresult{3}(2, :), stepresult{3}(3, :));
+    end
+    
+    % finish plot
+    hold off;
+    legend(stepresults{:, 2});
+%    %print(['report/', func2str(algfunc), 'traj'], '-dpdf');
+end
+
+% solve ODE using RK4 with automatic step size
+[result, sizes, errors] = rk4auto(sysfuncts, initvalues, interval, ...
+    1e-5, 10e-10, 10e-10);
+
+% plot trajectory
+figure;
+plot(result(2, :), result(3, :));
+grid on;
+title('RK4 auto-step trajectory plot (x_2 against x_1)');
+set(gcf, 'PaperPosition', [0 0 6 4]);
+set(gcf, 'PaperSize', [6 4]);
+%%print(['report/', 'rk4autotraj'], '-dpdf');
+
+% plot statistics
+stats = {
+    "RK4 auto-step step size", "rk4sizes", sizes;
+    "RK4 auto-step approximation error", "rk4errors", errors
+};
+for stat = stats'
+    figure;
+    plot(result(1, 2:(end - 1)), stat{3});
+    grid on;
+    title(stat{1});
+    set(gcf, 'PaperPosition', [0 0 6 4]);
+    set(gcf, 'PaperSize', [6 4]);
+%    %print(strcat("report/", stat{2}), '-dpdf');
+end
+
+% compare results with ODE45
+odefun = @(t, x) [ sysfuncts{1}(x); sysfuncts{2}(x) ];
+odeoptions = odeset('RelTol', 10e-10, 'AbsTol', 10e-10);
+[t, x] = ode45(odefun, interval, initvalues, odeoptions);
+figure;
+plot(x(:, 1), x(:, 2));
+grid on;
+title('ODE45 trajectory plot (x_2 against x_1)');
+set(gcf, 'PaperPosition', [0 0 6 4]);
+set(gcf, 'PaperSize', [6 4]);
+%print('report/ode45', '-dpdf');

ENUME/projectC/task2ode.m

@@ -0,0 +1,104 @@
+% functions of the ODE system and initial values
+sysfuncts = {
+    @(x) x(2) + x(1) * (0.5 - x(1)^2 - x(2)^2);
+    @(x) -x(1) + x(2) * (0.5 - x(1)^2 - x(2)^2)
+};
+initvalues = [0; 12];
+interval = [0; 15];
+
+% define available algorithms
+algorithms = {
+    'RK4', @rk4, [0.01, 0.013408];
+    'Adams PC', @adamspc, [0.002, 0.01305]
+};
+
+% solve ODE using different algorithms and step sizes
+for alg = 1:size(algorithms, 1)
+    [algname, algfunc, stepsizes] = algorithms{alg, :};
+    
+    % solve using the given algorithm for each step size
+    stepresults = cell(size(stepsizes, 2), 3);
+    stepnames = {'optimal step', 'larger step'};
+    for stepno = 1:size(stepsizes, 2)
+        result = algfunc(sysfuncts, initvalues, interval, stepsizes(stepno));
+        stepresults(stepno, :) = { ...
+            stepsizes(stepno), ...
+            stepnames{stepno}, ...
+            result};
+    end
+    
+    % plot each component against time
+    for eqnum = 1:size(sysfuncts, 1)
+        % begin plot
+        figure; grid on; hold on;
+        title([algname, ', x_', num2str(eqnum), ' against time']);
+        set(gcf, 'PaperPosition', [0 0 6 4]);
+        set(gcf, 'PaperSize', [6 4]);
+        
+        % plot component for each step size
+        for stepresult = stepresults'
+            plot(stepresult{3}(1, :), stepresult{3}(eqnum + 1, :));
+        end
+        
+        % finish plot
+        hold off;
+        legend(stepresults{:, 2});
+%        %print(['report/', func2str(algfunc), 'x', num2str(eqnum)], '-dpdf');
+    end
+    
+    % plot first two components against each other
+    figure; grid on; hold on;
+    title([algname, ' trajectory plot (x_2 against x_1)']);
+    set(gcf, 'PaperPosition', [0 0 6 4]);
+    set(gcf, 'PaperSize', [6 4]);
+    
+    % plot for each step size
+    for stepresult = stepresults'
+        plot(stepresult{3}(2, :), stepresult{3}(3, :));
+    end
+    
+    % finish plot
+    hold off;
+    legend(stepresults{:, 2});
+%    %print(['report/', func2str(algfunc), 'traj'], '-dpdf');
+end
+
+% solve ODE using RK4 with automatic step size
+[result, sizes, errors] = rk4auto(sysfuncts, initvalues, interval, ...
+    1e-5, 10e-10, 10e-10);
+
+% plot trajectory
+figure;
+plot(result(2, :), result(3, :));
+grid on;
+title('RK4 auto-step trajectory plot (x_2 against x_1)');
+set(gcf, 'PaperPosition', [0 0 6 4]);
+set(gcf, 'PaperSize', [6 4]);
+%%print(['report/', 'rk4autotraj'], '-dpdf');
+
+% plot statistics
+stats = {
+    "RK4 auto-step step size", "rk4sizes", sizes;
+    "RK4 auto-step approximation error", "rk4errors", errors
+};
+for stat = stats'
+    figure;
+    plot(result(1, 2:(end - 1)), stat{3});
+    grid on;
+    title(stat{1});
+    set(gcf, 'PaperPosition', [0 0 6 4]);
+    set(gcf, 'PaperSize', [6 4]);
+%    %print(strcat("report/", stat{2}), '-dpdf');
+end
+
+% compare results with ODE45
+odefun = @(t, x) [ sysfuncts{1}(x); sysfuncts{2}(x) ];
+odeoptions = odeset('RelTol', 10e-10, 'AbsTol', 10e-10);
+[t, x] = ode45(odefun, interval, initvalues, odeoptions);
+figure;
+plot(x(:, 1), x(:, 2));
+grid on;
+title('ODE45 trajectory plot (x_2 against x_1)');
+set(gcf, 'PaperPosition', [0 0 6 4]);
+set(gcf, 'PaperSize', [6 4]);
+%print('report/ode45', '-dpdf');