From 245acae5b6f61df8fecd6a3f9acff25b918a1ae0 Mon Sep 17 00:00:00 2001 From: PolishPigeon Date: Sat, 8 Jan 2022 05:36:00 +0100 Subject: [PATCH] feat: beautyfying task2.m and finishing code part --- ENUME/projectC/{adamspc.m => AdamsPCMethod.m} | 2 +- ENUME/projectC/adamspc.asv | 70 ------- ENUME/projectC/auxplot.m | 37 ---- ENUME/projectC/task2.asv | 116 ++++++++++++ ENUME/projectC/task2.m | 176 ++++++++++-------- 5 files changed, 218 insertions(+), 183 deletions(-) rename ENUME/projectC/{adamspc.m => AdamsPCMethod.m} (98%) delete mode 100644 ENUME/projectC/adamspc.asv delete mode 100644 ENUME/projectC/auxplot.m create mode 100644 ENUME/projectC/task2.asv diff --git a/ENUME/projectC/adamspc.m b/ENUME/projectC/AdamsPCMethod.m similarity index 98% rename from ENUME/projectC/adamspc.m rename to ENUME/projectC/AdamsPCMethod.m index 2ba4838d..f9ff6e4b 100644 --- a/ENUME/projectC/adamspc.m +++ b/ENUME/projectC/AdamsPCMethod.m @@ -1,4 +1,4 @@ -function x = adamspc(functions, initialValues, interval, stepSize) +function x = AdamsPCMethod(functions, initialValues, interval, stepSize) [x, derrivatives, explicitCoefficients, implicitCoefficients, stepCount] = initialize(functions, initialValues, interval, stepSize); x = adamsPcLoop(x, derrivatives, explicitCoefficients, implicitCoefficients, stepCount, functions, stepSize); % append arguments to output diff --git a/ENUME/projectC/adamspc.asv b/ENUME/projectC/adamspc.asv deleted file mode 100644 index a10670a6..00000000 --- a/ENUME/projectC/adamspc.asv +++ /dev/null @@ -1,70 +0,0 @@ -function x = adamspc(functions, initialValues, interval, stepSize) - [x, derrivatives, explicitCoefficients, implicitCoefficients, stepCount] = initialize(functions, initialValues, interval, stepSize); - [stepSize, stepCount, x] = adamsPcLoop(x, derrivatives, explicitCoefficients, implicitCoefficients, stepCount, functions, stepSize); - % append arguments to output - x = [interval(1):stepSize:(stepCount * stepSize); x]; -end - -function [x, derrivatives, explicitCoefficients, implicitCoefficients, stepCount] = initialize(functions, initialValues, interval, stepSize) - % obtain first five steps from RK4 - [x, derrivatives] = RK4(functions, initialValues, interval, stepSize, 5); - x = x(2:end, :); - - % define coefficient - explicitCoefficients = [1901, -2774, 1616, -1274, 251] / 720; - % Constants that can be found on the Internet or in Mister Tatjewski - % book on page 177, notice how beta(3) = 1616 instead of 2616 - % I have found on the internet different value for this parameter and - % got better results with it so I decided to stick with it - implicitCoefficients = [475, 1427, -798, 482, -173, 27] / 1440; - % Constants that can be found on the Internet or in Mister Tatjewski - % book on page 178 - % build output based on preceding values - stepCount = ceil((interval(2) - interval(1)) / stepSize); -end - -function [stepSize, stepCount, x] = adamsPcLoop(x, derrivatives, explicitCoefficients, implicitCoefficients, stepCount, functions, stepSize) - for step = 6: (stepCount + 1) - % P - predictionOfX = adamsPredict(x, step, stepSize, explicitCoefficients, derrivatives); - - - - % correct - x(:, step) = x(:, step - 1); - for equationNumber = 1:size(functions, 1) - for previous = 1:5 - x(equationNumber, step) = x(equationNumber, step) ... - + stepSize * implicitCoefficients(previous + 1) * derrivatives(equationNumber, step - previous); - end - x(equationNumber, step) = x(equationNumber, step) ... - + stepSize * implicitCoefficients(1) * derrivativePrediction(equationNumber); - end - - % evaluate - for equationNumber = 1:size(functions, 1) - derrivatives(equationNumber, step) = functions{equationNumber}(x(:, step)); - end - end -end - -function predictionOfX = adamsPredict(x, step, stepSize, explicitCoefficients, derrivatives) - predictionOfX = x(:, step - 1); - for equationNumber = 1 : 2 - for previous = 1:5 - predictionOfX(equationNumber) = predictionOfX(equationNumber) ... - + stepSize * explicitCoefficients(previous) * derrivatives(equationNumber, step - previous); - end - end -end - -function evaluateAdamsValue = adamsEvaluate - evaluateAdamsValue = zeros(size(functions, 1), 1); - for equationNumber = 1:size(functions, 1) - evaluateAdamsValue(equationNumber) = functions{equationNumber}(predictionOfX); - end -end - -function adamsCorrect - -end \ No newline at end of file diff --git a/ENUME/projectC/auxplot.m b/ENUME/projectC/auxplot.m deleted file mode 100644 index 69b1eb45..00000000 --- a/ENUME/projectC/auxplot.m +++ /dev/null @@ -1,37 +0,0 @@ -% 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/task2.asv b/ENUME/projectC/task2.asv new file mode 100644 index 00000000..b102f533 --- /dev/null +++ b/ENUME/projectC/task2.asv @@ -0,0 +1,116 @@ + +function task2() + [ordinaryDifferentialEquations, initialValues, interval, algorithms] = initialize(); + solveAndPrintODEs(algorithms, ordinaryDifferentialEquations, initialValues, interval); + [result, sizes, errors] = solveAndPlotODEAutomatic(ordinaryDifferentialEquations, initialValues, interval); + plotStatistics(result, sizes, errors); + plotComparisonToMatlabFunction(ordinaryDifferentialEquations, interval, initialValues); +end + +function [ordinaryDifferentialEquations, initialValues, interval, algorithms] = initialize() + ordinaryDifferentialEquations = { + @(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) + }; + initialValues = [8; 9]; + interval = [0; 15]; + + algorithms = { + 'RK4 algorithm', @RK4, [0.01, 0.011]; + 'Adams PC algorithm', @AdamsPCMethod, [0.002, 0.01305] + }; +end + +function solveAndPrintODEs(algorithms, ordinaryDifferentialEquations, initialValues, interval) + for algorithm = 1 : 2 + [algorithmName, algorithmFunction, stepSizes] = algorithms{algorithm, :}; + stepResults = solveForEachStep(stepSizes, ordinaryDifferentialEquations, initialValues, interval, algorithmFunction); + plotAgainstTime(ordinaryDifferentialEquations, algorithmName, stepResults); + plotAgainst(algorithmName, stepResults); + + end +end + +function stepResults = solveForEachStep(stepSizes, ordinaryDifferentialEquations, initialValues, interval, algorithmFunction) + stepResults = cell(size(stepSizes, 2), 3); + stepNames = {'Optimal Step Size', 'Slightly Larger Step Size'}; + for stepNumber = 1:size(stepSizes, 2) + result = algorithmFunction(ordinaryDifferentialEquations, initialValues, interval, stepSizes(stepNumber)); + stepResults(stepNumber, :) = {stepSizes(stepNumber), stepNames{stepNumber}, result}; + end +end + +function beginPlot() + figure; + grid on; + hold on; +end + +function plotAgainstTime(ordinaryDifferentialEquations, algorithmName, stepResults) + for equationNumber = 1:size(ordinaryDifferentialEquations, 1) + beginPlot(); + + title([algorithmName, ', x_', num2str(equationNumber), ' ODE']); + + for stepresult = stepResults' + plot(stepresult{3}(1, :), stepresult{3}(equationNumber + 1, :)); + end + + hold off; + legend(stepResults{:, 2}); + end +end + +function plotAgainst(algorithmName, stepResults) + beginPlot(); + title([algorithmName, ' trajectory plot (x_2 compared to x_1)']); + for stepresult = stepResults' + plot(stepresult{3}(2, :), stepresult{3}(3, :)); + end + hold off; + legend(stepResults{:, 2}); +% %print(['report/', func2str(algfunc), 'traj'], '-dpdf'); +end + +function [result, sizes, errors] = solveAndPlotODEAutomatic(ordinaryDifferentialEquations, initialValues, interval) + [result, sizes, errors] = solveRKAutomatic(ordinaryDifferentialEquations, initialValues, interval); + plotTrajectory(result); +end + +function [result, sizes, errors] = solveRKAutomatic(ordinaryDifferentialEquations, initialValues, interval) + initialStepSize = 1e-5; + relativeEpsilon = 1e-9; + absoluteEpsilon = 1e-9; + [result, sizes, errors] = RK4Automatic(ordinaryDifferentialEquations, initialValues, interval, initialStepSize, relativeEpsilon, absoluteEpsilon); +end + +function plotTrajectory(result) + figure; + plot(result(2, :), result(3, :)); + grid on; + title('RK4 auto-step trajectory plot (x_2 against x_1)'); +end + +function plotStatistics(result, sizes, errors) + 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}); + end +end + +function plotComparisonToMatlabFunction(ordinaryDifferentialEquations, interval, initialValues) + % compare results with ODE45 + odefun = @(t, x) [ ordinaryDifferentialEquations{1}(x); ordinaryDifferentialEquations{2}(x) ]; + odeoptions = odeset('RelTol', 10e-10, 'AbsTol', 10e-10); + [~, x] = ode45(odefun, interval, initialValues, odeoptions); + figure; + plot(x(:, 1), x(:, 2)); + grid on; + title('ODE45 trajectory plot (x_2 against x_1)'); +end \ No newline at end of file diff --git a/ENUME/projectC/task2.m b/ENUME/projectC/task2.m index ce908f1c..ccc58b6c 100644 --- a/ENUME/projectC/task2.m +++ b/ENUME/projectC/task2.m @@ -1,90 +1,116 @@ -% 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] -}; +function task2() + [ordinaryDifferentialEquations, initialValues, interval, algorithms] = initialize(); + solveAndPrintODEs(algorithms, ordinaryDifferentialEquations, initialValues, interval); + [result, sizes, errors] = solveAndPlotODEAutomatic(ordinaryDifferentialEquations, initialValues, interval); + plotStatistics(result, sizes, errors); + plotComparisonToMatlabFunction(ordinaryDifferentialEquations, interval, initialValues); +end -% solve ODE using different algorithms and step sizes -for alg = 1 : 2 - [algname, algfunc, stepsizes] = algorithms{alg, :}; +function [ordinaryDifferentialEquations, initialValues, interval, algorithms] = initialize() + ordinaryDifferentialEquations = { + @(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) + }; + initialValues = [8; 9]; + interval = [0; 15]; - % 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}; + algorithms = { + 'RK4 algorithm', @RK4, [0.01, 0.011]; + 'Adams PC algorithm', @AdamsPCMethod, [0.002, 0.01305] + }; +end + +function solveAndPrintODEs(algorithms, ordinaryDifferentialEquations, initialValues, interval) + for algorithm = 1 : 2 + [algorithmName, algorithmFunction, stepSizes] = algorithms{algorithm, :}; + stepResults = solveForEachStep(stepSizes, ordinaryDifferentialEquations, initialValues, interval, algorithmFunction); + plotAgainstTime(ordinaryDifferentialEquations, algorithmName, stepResults); + plotAgainst(algorithmName, stepResults); + 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']); - - % plot component for each step size - for stepresult = stepresults' - plot(stepresult{3}(1, :), stepresult{3}(eqnum + 1, :)); +end + +function stepResults = solveForEachStep(stepSizes, ordinaryDifferentialEquations, initialValues, interval, algorithmFunction) + stepResults = cell(size(stepSizes, 2), 3); + stepNames = {'Optimal Step Size', 'Slightly Larger Step Size'}; + for stepNumber = 1:size(stepSizes, 2) + result = algorithmFunction(ordinaryDifferentialEquations, initialValues, interval, stepSizes(stepNumber)); + stepResults(stepNumber, :) = {stepSizes(stepNumber), stepNames{stepNumber}, result}; + end +end + +function beginPlot() + figure; + grid on; + hold on; +end + +function plotAgainstTime(ordinaryDifferentialEquations, algorithmName, stepResults) + for equationNumber = 1:size(ordinaryDifferentialEquations, 1) + beginPlot(); + + title([algorithmName, ', x_', num2str(equationNumber), ' ODE']); + + for stepresult = stepResults' + plot(stepresult{3}(1, :), stepresult{3}(equationNumber + 1, :)); end - - % finish plot + hold off; - legend(stepresults{:, 2}); + legend(stepResults{:, 2}); end - - % plot first two components against each other - figure; grid on; hold on; - title([algname, ' trajectory plot (x_2 against x_1)']); - - % plot for each step size - for stepresult = stepresults' +end + +function plotAgainst(algorithmName, stepResults) + beginPlot(); + title([algorithmName, ' trajectory plot (x_2 compared to x_1)']); + for stepresult = stepResults' plot(stepresult{3}(2, :), stepresult{3}(3, :)); end - - % finish plot hold off; - legend(stepresults{:, 2}); + legend(stepResults{:, 2}); % %print(['report/', func2str(algfunc), 'traj'], '-dpdf'); end -% solve ODE using RK4 with automatic step size -[result, sizes, errors] = RK4Automatic(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)'); - -% 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}); +function [result, sizes, errors] = solveAndPlotODEAutomatic(ordinaryDifferentialEquations, initialValues, interval) + [result, sizes, errors] = solveRKAutomatic(ordinaryDifferentialEquations, initialValues, interval); + plotTrajectory(result); 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)'); +function [result, sizes, errors] = solveRKAutomatic(ordinaryDifferentialEquations, initialValues, interval) + initialStepSize = 1e-5; + relativeEpsilon = 1e-9; + absoluteEpsilon = 1e-9; + [result, sizes, errors] = RK4Automatic(ordinaryDifferentialEquations, initialValues, interval, initialStepSize, relativeEpsilon, absoluteEpsilon); +end + +function plotTrajectory(result) + figure; + plot(result(2, :), result(3, :)); + grid on; + title('RK4 with auto step trajectory plot (x_2 against x_1)'); +end + +function plotStatistics(result, sizes, errors) + stats = { + "RK4 with auto step step size", "rk4sizes", sizes; + "RK4 with auto step approximation error", "rk4errors", errors + }; + for stat = stats' + figure; + plot(result(1, 2:(end - 1)), stat{3}); + grid on; + title(stat{1}); + end +end + +function plotComparisonToMatlabFunction(ordinaryDifferentialEquations, interval, initialValues) + % compare results with ODE45 + odefun = @(t, x) [ ordinaryDifferentialEquations{1}(x); ordinaryDifferentialEquations{2}(x) ]; + odeoptions = odeset('RelTol', 10e-10, 'AbsTol', 10e-10); + [~, x] = ode45(odefun, interval, initialValues, odeoptions); + figure; + plot(x(:, 1), x(:, 2)); + grid on; + title('ODE45 trajectory plot (x_2 against x_1)'); +end \ No newline at end of file