From cb42f561f5b7ad27aba88b283e33dbfb81667e85 Mon Sep 17 00:00:00 2001 From: PolishPigeon Date: Sat, 8 Jan 2022 00:52:14 +0100 Subject: [PATCH] feat: beautified RK4Phi and rk4.m functions --- ENUME/projectC/RK4Phi.m | 8 +++ ENUME/projectC/adamspc.m | 47 +++++++++++++++++ ENUME/projectC/auxplot.m | 41 +++++++++++++++ ENUME/projectC/rk4.asv | 47 +++++++++++++++++ ENUME/projectC/rk4.m | 70 ++++++++++++++++---------- ENUME/projectC/rk4auto.m | 4 +- ENUME/projectC/rk4phi.m | 8 --- ENUME/projectC/task1lsapprox.m | 66 ------------------------ ENUME/projectC/{task2ode.m => task2.m} | 16 +----- 9 files changed, 190 insertions(+), 117 deletions(-) create mode 100644 ENUME/projectC/RK4Phi.m create mode 100644 ENUME/projectC/adamspc.m create mode 100644 ENUME/projectC/auxplot.m create mode 100644 ENUME/projectC/rk4.asv delete mode 100644 ENUME/projectC/rk4phi.m delete mode 100644 ENUME/projectC/task1lsapprox.m rename ENUME/projectC/{task2ode.m => task2.m} (80%) diff --git a/ENUME/projectC/RK4Phi.m b/ENUME/projectC/RK4Phi.m new file mode 100644 index 00000000..afdc8d38 --- /dev/null +++ b/ENUME/projectC/RK4Phi.m @@ -0,0 +1,8 @@ +function Phi = RK4Phi(algorithm, stepValue, stepSize) + k1 = algorithm(stepValue); + k2 = algorithm(stepValue + 0.5 * stepSize * k1); + k3 = algorithm(stepValue + 0.5 * stepSize * k2); + k4 = algorithm(stepValue + stepSize * k3); + Phi = (k1 + 2 * k2 + 2 * k3 + k4) / 6; +end +% calculates the Phi used in RK4 algorithms diff --git a/ENUME/projectC/adamspc.m b/ENUME/projectC/adamspc.m new file mode 100644 index 00000000..83018aa1 --- /dev/null +++ b/ENUME/projectC/adamspc.m @@ -0,0 +1,47 @@ +function x = adamspc(functs, init, interval, stepsize) + % obtain first five steps from RK4 + [x, derivs] = rk4(functs, init, interval, stepsize, 5); + x = x(2:end, :); + + % define coefficient tables + beta = [1901, -2774, 2616, -1274, 251] / 720; + betastar = [475, 1427, -798, 482, -173, 27] / 1440; + + % build output based on preceding values + stepcount = ceil((interval(2) - interval(1)) / stepsize); + for step = 6:(stepcount + 1) + % predict + xpred = x(:, step - 1); + for eqnum = 1:size(functs, 1) + for prev = 1:5 + xpred(eqnum) = xpred(eqnum) ... + + stepsize * beta(prev) * derivs(eqnum, step - prev); + end + end + + % evaluate + derivpred = zeros(size(functs, 1), 1); + for eqnum = 1:size(functs, 1) + derivpred(eqnum) = functs{eqnum}(xpred); + end + + % correct + x(:, step) = x(:, step - 1); + for eqnum = 1:size(functs, 1) + for prev = 1:5 + x(eqnum, step) = x(eqnum, step) ... + + stepsize * betastar(prev + 1) * derivs(eqnum, step - prev); + end + x(eqnum, step) = x(eqnum, step) ... + + stepsize * betastar(1) * derivpred(eqnum); + end + + % evaluate + for eqnum = 1:size(functs, 1) + derivs(eqnum, step) = functs{eqnum}(x(:, step)); + end + end + + % append arguments to output + x = [interval(1):stepsize:(stepcount * stepsize); x]; +end diff --git a/ENUME/projectC/auxplot.m b/ENUME/projectC/auxplot.m new file mode 100644 index 00000000..b8012870 --- /dev/null +++ b/ENUME/projectC/auxplot.m @@ -0,0 +1,41 @@ +% ENUME MICHAŁ SZOPIŃSKI +% PROJECT C NUMBER 60 +% https://github.com/Lachcim/szopinski-enume + +% coordinate arguments +argsx = -1:0.1:1; +argsy = -2:0.5:12; + +% available functions +functs = { + "x_1", -67, 50, @(x, y) y + x * (0.5 - x^2 - y^2); + "x_2", -115, 30, @(x, y) -x + y * (0.5 - x^2 - y^2) +}; + +% plot each function +for funct = 1:size(functs, 1) + [functname, cama, camb, functh] = functs{funct, :}; + + % evaluate values + vals = zeros(length(argsx), length(argsy)); + for i = 1:size(argsx, 2) + for j = 1:size(argsy, 2) + vals(i, j) = functh(argsx(i), argsy(j)); + end + end + + % plot + figure; + surf(argsx, argsy, vals'); + title(strcat("Gradient of ", functname)); + xlabel('x_1'); + ylabel('x_2'); + zlabel(strcat("d", functname, "/dt")); + xlim([argsx(1), argsx(end)]); + ylim([argsy(1), argsy(end)]); + view(cama, camb); + grid on; + set(gcf, 'PaperPosition', [0 0 6 4]); + set(gcf, 'PaperSize', [6 4]); + print(strcat('report/', functname, 'gradient'), '-dpdf'); +end diff --git a/ENUME/projectC/rk4.asv b/ENUME/projectC/rk4.asv new file mode 100644 index 00000000..ba1e8c93 --- /dev/null +++ b/ENUME/projectC/rk4.asv @@ -0,0 +1,47 @@ +% solve ODE system using RK4 with constant step size +function [x, derivativesTable] = rk4(functions, initialValues, interval, stepSize, maxSteps) + x = initialValues; + + derivativesTable = buildDerivatiesTable(x, functions); + + % Calculate stepCount + stepCount = ceil((interval(2) - interval(1)) / stepSize); + if nargin == 5 + stepCount = min(stepCount, maxSteps - 1); + end % IF we include max steps in our function input + % (nargin is number of arguments in input) + % then choose smaller number between maxSteps and stepCount and choose + % it for stepCount + + [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, functions, derivativesTable); + + + % append arguments to output + x = [interval(1):stepSize:(stepCount * stepSize); x]; +end + +function derivativesTable = buildDerivatiesTable(x, functions) + derivativesTable = zeros(size(x)); + for eqnum = 1:size(functions, 1) + derivativesTable(eqnum, 1) = functions{eqnum}(x(:, 1)); + end +end + +function [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, functions, derivativesTable) + for step = 1 : stepCount + [x, derivativesTable] = rk4stepLoop(x, step, functions, stepSize, derivativesTable); + end +end + +function [x, derivativesTable] = rk4stepLoop(x, step, functions, stepSize, derivativesTable) + stepValue = x(:, step); + + for equationNumber = 1 : 2 + % generic single-step iteration + phi = RK4Phi(functions{equationNumber}, stepValue, stepSize); + x(equationNumber, step + 1) = x(equationNumber, step) + stepSize * phi; + + % update derivatives table + derivativesTable(equationNumber, step + 1) = functions{equationNumber}(x(:, step + 1)); + end +end \ No newline at end of file diff --git a/ENUME/projectC/rk4.m b/ENUME/projectC/rk4.m index 80fb5be6..0f880480 100644 --- a/ENUME/projectC/rk4.m +++ b/ENUME/projectC/rk4.m @@ -1,31 +1,49 @@ % 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; +function [x, derivativesTable] = rk4(equations, initialValues, interval, stepSize, maxSteps) + x = initialValues; + + derivativesTable = buildDerivatiesTable(x, equations); - % 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 + % Calculate stepCount + stepCount = ceil((interval(2) - interval(1)) / stepSize); + if nargin == 5 + stepCount = min(stepCount, maxSteps - 1); + end % IF we include max steps in our function input + % (nargin is number of arguments in input) + % then choose smaller number between maxSteps and stepCount and choose + % it for stepCount + + [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, equations, derivativesTable); + % append arguments to output - x = [interval(1):stepsize:(stepcount * stepsize); x]; + x = [interval(1):stepSize:(stepCount * stepSize); x]; end + +function derivativesTable = buildDerivatiesTable(x, equations) + derivativesTable = zeros(size(x)); + for eqnum = 1:size(equations, 1) + derivativesTable(eqnum, 1) = equations{eqnum}(x(:, 1)); + end +end + +function [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, equations, derivativesTable) + for step = 1 : stepCount + [x, derivativesTable] = rk4stepLoop(x, step, equations, stepSize, derivativesTable); + end +end + +function [x, derivativesTable] = rk4stepLoop(x, step, equations, stepSize, derivativesTable) + stepValue = x(:, step); + [x, derivativesTable] = equationsLoop(x, equations, stepValue, stepSize, step, derivativesTable); + +end + +function [x, derivativesTable] = equationsLoop(x, equations, stepValue, stepSize, step, derivativesTable) + for equationNumber = 1 : 2 + Phi = RK4Phi(equations{equationNumber}, stepValue, stepSize); + x(equationNumber, step + 1) = x(equationNumber, step) + stepSize * Phi; + + derivativesTable(equationNumber, step + 1) = equations{equationNumber}(x(:, step + 1)); + end +end \ No newline at end of file diff --git a/ENUME/projectC/rk4auto.m b/ENUME/projectC/rk4auto.m index 154283c7..df2b198e 100644 --- a/ENUME/projectC/rk4auto.m +++ b/ENUME/projectC/rk4auto.m @@ -19,7 +19,7 @@ function [x, sizes, errors] = rk4auto(functs, init, interval, initstep, eps_rel, % advance output function for eqnum = 1:size(functs, 1) % generic single-step iteration - phi = rk4phi(functs{eqnum}, stepval, stepsize); + phi = RK4Phi(functs{eqnum}, stepval, stepsize); x(eqnum, step + 1) = x(eqnum, step) + stepsize * phi; end @@ -30,7 +30,7 @@ function [x, sizes, errors] = rk4auto(functs, init, interval, initstep, eps_rel, % 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); + phi = RK4Phi(functs{eqnum}, stepval, stepsize / 2); stepval(eqnum) = stepval(eqnum) + (stepsize / 2) * phi; end end diff --git a/ENUME/projectC/rk4phi.m b/ENUME/projectC/rk4phi.m deleted file mode 100644 index 7abc7caa..00000000 --- a/ENUME/projectC/rk4phi.m +++ /dev/null @@ -1,8 +0,0 @@ -% 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 diff --git a/ENUME/projectC/task1lsapprox.m b/ENUME/projectC/task1lsapprox.m deleted file mode 100644 index bcc2ee3d..00000000 --- a/ENUME/projectC/task1lsapprox.m +++ /dev/null @@ -1,66 +0,0 @@ -function task1() - for polydeg = 0:3 - dataPoints = functionDataPoints(); - % obtain factors of approximating polynomial - [factors, error, gramcond] = approximate(dataPoints, 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(dataPoints(:, 1), dataPoints(:, 2)); - - % plot approximation - x = dataPoints(1):0.05:dataPoints(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 -end - - - -% find the approximating polynomial of the given degree -function [factors, error, gramcond] = approximate(dataPoints, polydeg) - % define the A matrix used for solving the error minimization problem - A = zeros(size(dataPoints, 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) = dataPoints(row, 1) ^ (col - 1); - end - end - - % solve least-square problem using QRdash decomposition - [Q, eqsys, invqtq] = QRDecomposition(A); - eqsys(:, end + 1) = invqtq * Q' * dataPoints(:, 2); - factors = backSubstitution(eqsys); - - % calculate error and condition number of Gram's matrix - error = norm(dataPoints(:, 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 diff --git a/ENUME/projectC/task2ode.m b/ENUME/projectC/task2.m similarity index 80% rename from ENUME/projectC/task2ode.m rename to ENUME/projectC/task2.m index 6994c99d..80a99b3a 100644 --- a/ENUME/projectC/task2ode.m +++ b/ENUME/projectC/task2.m @@ -13,7 +13,7 @@ algorithms = { }; % solve ODE using different algorithms and step sizes -for alg = 1:size(algorithms, 1) +for alg = 1 : 2 [algname, algfunc, stepsizes] = algorithms{alg, :}; % solve using the given algorithm for each step size @@ -32,8 +32,6 @@ for alg = 1:size(algorithms, 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' @@ -43,14 +41,11 @@ for alg = 1:size(algorithms, 1) % 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' @@ -72,9 +67,6 @@ 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 = { @@ -86,9 +78,6 @@ for stat = stats' 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 @@ -99,6 +88,3 @@ 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');