From 3ae8bd014c11129d0e0d5d4bee274bf9a38c8a53 Mon Sep 17 00:00:00 2001 From: PolishPigeon Date: Sat, 8 Jan 2022 03:24:09 +0100 Subject: [PATCH] Beautifying RK4Automatic.m --- ENUME/projectC/{rk4.asv => RK4.m} | 32 ++++++----- ENUME/projectC/RK4Automatic.asv | 92 +++++++++++++++++++++++++++++ ENUME/projectC/RK4Automatic.m | 96 +++++++++++++++++++++++++++++++ ENUME/projectC/auxplot.m | 4 -- ENUME/projectC/equationsLoop.m | 0 ENUME/projectC/rk4.m | 2 +- ENUME/projectC/rk4auto.m | 18 +++--- ENUME/projectC/task2.m | 4 +- 8 files changed, 217 insertions(+), 31 deletions(-) rename ENUME/projectC/{rk4.asv => RK4.m} (53%) create mode 100644 ENUME/projectC/RK4Automatic.asv create mode 100644 ENUME/projectC/RK4Automatic.m create mode 100644 ENUME/projectC/equationsLoop.m diff --git a/ENUME/projectC/rk4.asv b/ENUME/projectC/RK4.m similarity index 53% rename from ENUME/projectC/rk4.asv rename to ENUME/projectC/RK4.m index ba1e8c93..7b6bb04e 100644 --- a/ENUME/projectC/rk4.asv +++ b/ENUME/projectC/RK4.m @@ -1,8 +1,8 @@ % solve ODE system using RK4 with constant step size -function [x, derivativesTable] = rk4(functions, initialValues, interval, stepSize, maxSteps) +function [x, derivativesTable] = RK4(equations, initialValues, interval, stepSize, maxSteps) x = initialValues; - derivativesTable = buildDerivatiesTable(x, functions); + derivativesTable = buildDerivatiesTable(x, equations); % Calculate stepCount stepCount = ceil((interval(2) - interval(1)) / stepSize); @@ -13,35 +13,37 @@ function [x, derivativesTable] = rk4(functions, initialValues, interval, stepSiz % then choose smaller number between maxSteps and stepCount and choose % it for stepCount - [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, functions, derivativesTable); + [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, equations, derivativesTable); % append arguments to output x = [interval(1):stepSize:(stepCount * stepSize); x]; end -function derivativesTable = buildDerivatiesTable(x, functions) +function derivativesTable = buildDerivatiesTable(x, equations) derivativesTable = zeros(size(x)); - for eqnum = 1:size(functions, 1) - derivativesTable(eqnum, 1) = functions{eqnum}(x(:, 1)); + for eqnum = 1:size(equations, 1) + derivativesTable(eqnum, 1) = equations{eqnum}(x(:, 1)); end end -function [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, functions, derivativesTable) +function [x, derivativesTable] = rk4Loop(x, stepCount, stepSize, equations, derivativesTable) for step = 1 : stepCount - [x, derivativesTable] = rk4stepLoop(x, step, functions, stepSize, derivativesTable); + [x, derivativesTable] = rk4stepLoop(x, step, equations, stepSize, derivativesTable); end end -function [x, derivativesTable] = rk4stepLoop(x, step, functions, stepSize, derivativesTable) +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 - % generic single-step iteration - phi = RK4Phi(functions{equationNumber}, stepValue, stepSize); - x(equationNumber, step + 1) = x(equationNumber, step) + stepSize * phi; + Phi = RK4Phi(equations{equationNumber}, stepValue, stepSize); + x(equationNumber, step + 1) = x(equationNumber, step) + stepSize * Phi; - % update derivatives table - derivativesTable(equationNumber, step + 1) = functions{equationNumber}(x(:, step + 1)); + derivativesTable(equationNumber, step + 1) = equations{equationNumber}(x(:, step + 1)); end end \ No newline at end of file diff --git a/ENUME/projectC/RK4Automatic.asv b/ENUME/projectC/RK4Automatic.asv new file mode 100644 index 00000000..e74ee104 --- /dev/null +++ b/ENUME/projectC/RK4Automatic.asv @@ -0,0 +1,92 @@ +% automatic step size variant of RK4 +function [x, sizes, errors] = RK4Automatic(equations, initialValues, interval, initialStepSize, relativeEpsilon, absoluteEpsilon) + + [functionArguments, x, sizes, errors, stepSize, step, flag] = Initialize(interval, initialValues, initialStepSize); + + [functionArguments, x, sizes, errors] = RK4AutomaticLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag); + % append arguments to output + x = [functionArguments; x]; +end + +function [functionArguments, x, sizes, errors, stepSize, step, flag] = Initialize(interval, initialValues, initialStepSize) + functionArguments = interval(1); + x = initialValues; + + sizes = double.empty(); + errors = double.empty(); + + stepSize = initialStepSize; + step = 0; + flag = 1; +end + +function [functionArguments, x, sizes, errors] = RK4AutomaticLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag) + while flag + [functionArguments, x, sizes, errors, stepSize, step, interval, flag] = insideWhileLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag); + end +end + +function [functionArguments, x, sizes, errors, stepSize, step, interval, flag] = insideWhileLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag) + [step, stepValue] = initializeSteps(step, x); + x = equationsLoop(x, equations, stepValue, stepSize, step); + [flag, functionArguments] = stopAlgorithm(functionArguments, stepSize, step, interval); + if ~flag + return; + end + stepValue = calculateNextStep(equations, stepValue, stepSize); + [stepCorrectionFactor, errors] = calculateStepCorrection(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors); + stepSize = 0.9 * stepCorrectionFactor * stepSize; + sizes(step) = stepSize; +end + +function v + +function [step, stepValue] = initializeSteps(step, x) + step = step + 1; + stepValue = x(:, step); +end + +function x = equationsLoop(x, equations, stepValue, stepSize, step) + for equationNumber = 1 : 2 + Phi = RK4Phi(equations{equationNumber}, stepValue, stepSize); + x(equationNumber, step + 1) = x(equationNumber, step) + stepSize * Phi; + end +end + +function [flag, functionArguments] = stopAlgorithm(functionArguments, stepSize, step, interval) + functionArguments(step + 1) = functionArguments(step) + stepSize; + if functionArguments(end) >= interval(2) + flag = 0; + else + flag = 1; + end +end + +function stepValue = calculateNextStep(equations, stepValue, stepSize) + for halfStep = 1 : 2 + for equationsNumber = 1:size(equations, 1) + Phi = RK4Phi(equations{equationsNumber}, stepValue, stepSize / 2); + stepValue(equationsNumber) = stepValue(equationsNumber) + (stepSize / 2) * Phi; + end + end +end + +function [stepCorrectionFactor, errors] = calculateStepCorrection(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors) + stepCorrectionFactor = Inf; % Initialize stepCorrectionFactor + [stepCorrectionFactor, errors] = calculateStepCorrectionLoop(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors, stepCorrectionFactor); + stepCorrectionFactor = stepCorrectionFactor ^ (1/5); +end + +function [stepCorrectionFactor, errors] = calculateStepCorrectionLoop(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors, stepCorrectionFactor) + for equationsNumber = 1 : 2 + approximationError = abs(stepValue(equationsNumber) - x(equationsNumber, step + 1)) / 15; + errors(step) = approximationError; + + epsilon = abs(stepValue(equationsNumber)) * relativeEpsilon + absoluteEpsilon; + equationAlpha = epsilon / approximationError; + + if equationAlpha < stepCorrectionFactor + stepCorrectionFactor = equationAlpha; + end + end +end diff --git a/ENUME/projectC/RK4Automatic.m b/ENUME/projectC/RK4Automatic.m new file mode 100644 index 00000000..a45495e3 --- /dev/null +++ b/ENUME/projectC/RK4Automatic.m @@ -0,0 +1,96 @@ +function [x, sizes, errors] = RK4Automatic(equations, initialValues, interval, initialStepSize, relativeEpsilon, absoluteEpsilon) + [functionArguments, x, sizes, errors, stepSize, step, flag] = Initialize(interval, initialValues, initialStepSize); + [functionArguments, x, sizes, errors] = RK4AutomaticLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag); + x = [functionArguments; x]; +end + +function [functionArguments, x, sizes, errors, stepSize, step, flag] = Initialize(interval, initialValues, initialStepSize) + functionArguments = interval(1); + x = initialValues; + + sizes = double.empty(); + errors = double.empty(); + + stepSize = initialStepSize; + step = 0; + flag = 1; +end + +function [functionArguments, x, sizes, errors] = RK4AutomaticLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag) + while flag + [functionArguments, x, sizes, errors, stepSize, step, interval, flag] = insideWhileLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag); + end +end + +function [functionArguments, x, sizes, errors, stepSize, step, interval, flag] = insideWhileLoop(functionArguments, x, sizes, errors, stepSize, step, equations, relativeEpsilon, absoluteEpsilon, interval, flag) + [step, stepValue, x, functionArguments, flag] = calculateXandFunctionArguments(step, x, equations, stepSize, functionArguments, interval); + if ~flag + return; + end + [stepSize, errors] = calculateStepAndErrors(equations, stepValue, stepSize, x, step, relativeEpsilon, absoluteEpsilon, errors); + sizes(step) = stepSize; + +end + +function [step, stepValue, x, functionArguments, flag] = calculateXandFunctionArguments(step, x, equations, stepSize, functionArguments, interval) + [step, stepValue] = initializeSteps(step, x); + x = equationsLoop(x, equations, stepValue, stepSize, step); + [flag, functionArguments] = stopAlgorithm(functionArguments, stepSize, step, interval); +end + +function [stepSize, errors] = calculateStepAndErrors(equations, stepValue, stepSize, x, step, relativeEpsilon, absoluteEpsilon, errors) + stepValue = calculateNextStep(equations, stepValue, stepSize); + [stepCorrectionFactor, errors] = calculateStepCorrection(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors); + stepSize = 0.9 * stepCorrectionFactor * stepSize; + +end + +function [step, stepValue] = initializeSteps(step, x) + step = step + 1; + stepValue = x(:, step); +end + +function x = equationsLoop(x, equations, stepValue, stepSize, step) + for equationNumber = 1 : 2 + Phi = RK4Phi(equations{equationNumber}, stepValue, stepSize); + x(equationNumber, step + 1) = x(equationNumber, step) + stepSize * Phi; + end +end + +function [flag, functionArguments] = stopAlgorithm(functionArguments, stepSize, step, interval) + functionArguments(step + 1) = functionArguments(step) + stepSize; + if functionArguments(end) >= interval(2) + flag = 0; + else + flag = 1; + end +end + +function stepValue = calculateNextStep(equations, stepValue, stepSize) + for halfStep = 1 : 2 + for equationsNumber = 1:size(equations, 1) + Phi = RK4Phi(equations{equationsNumber}, stepValue, stepSize / 2); + stepValue(equationsNumber) = stepValue(equationsNumber) + (stepSize / 2) * Phi; + end + end +end + +function [stepCorrectionFactor, errors] = calculateStepCorrection(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors) + stepCorrectionFactor = Inf; % Initialize stepCorrectionFactor + [stepCorrectionFactor, errors] = calculateStepCorrectionLoop(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors, stepCorrectionFactor); + stepCorrectionFactor = stepCorrectionFactor ^ (1/5); +end + +function [stepCorrectionFactor, errors] = calculateStepCorrectionLoop(stepValue, x, step, relativeEpsilon, absoluteEpsilon, errors, stepCorrectionFactor) + for equationsNumber = 1 : 2 + approximationError = abs(stepValue(equationsNumber) - x(equationsNumber, step + 1)) / 15; + errors(step) = approximationError; + + epsilon = abs(stepValue(equationsNumber)) * relativeEpsilon + absoluteEpsilon; + equationAlpha = epsilon / approximationError; + + if equationAlpha < stepCorrectionFactor + stepCorrectionFactor = equationAlpha; + end + end +end diff --git a/ENUME/projectC/auxplot.m b/ENUME/projectC/auxplot.m index b8012870..69b1eb45 100644 --- a/ENUME/projectC/auxplot.m +++ b/ENUME/projectC/auxplot.m @@ -1,7 +1,3 @@ -% 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; diff --git a/ENUME/projectC/equationsLoop.m b/ENUME/projectC/equationsLoop.m new file mode 100644 index 00000000..e69de29b diff --git a/ENUME/projectC/rk4.m b/ENUME/projectC/rk4.m index 0f880480..7b6bb04e 100644 --- a/ENUME/projectC/rk4.m +++ b/ENUME/projectC/rk4.m @@ -1,5 +1,5 @@ % solve ODE system using RK4 with constant step size -function [x, derivativesTable] = rk4(equations, initialValues, interval, stepSize, maxSteps) +function [x, derivativesTable] = RK4(equations, initialValues, interval, stepSize, maxSteps) x = initialValues; derivativesTable = buildDerivatiesTable(x, equations); diff --git a/ENUME/projectC/rk4auto.m b/ENUME/projectC/rk4auto.m index df2b198e..ba4bd559 100644 --- a/ENUME/projectC/rk4auto.m +++ b/ENUME/projectC/rk4auto.m @@ -1,15 +1,15 @@ % automatic step size variant of RK4 -function [x, sizes, errors] = rk4auto(functs, init, interval, initstep, eps_rel, eps_abs) +function [x, sizes, errors] = RK4Automatic(equations, initialValues, interval, initialStepSize, relativeEpsilon, absoluteEpsilon) % set start points of output args = interval(1); - x = init; + x = initialValues; % initialize output plots sizes = double.empty(); errors = double.empty(); % integrate function until end of interval reached - stepsize = initstep; + stepsize = initialStepSize; step = 0; while 1 % obtain the preceding function values @@ -17,9 +17,9 @@ function [x, sizes, errors] = rk4auto(functs, init, interval, initstep, eps_rel, stepval = x(:, step); % advance output function - for eqnum = 1:size(functs, 1) + for eqnum = 1:size(equations, 1) % generic single-step iteration - phi = RK4Phi(functs{eqnum}, stepval, stepsize); + phi = RK4Phi(equations{eqnum}, stepval, stepsize); x(eqnum, step + 1) = x(eqnum, step) + stepsize * phi; end @@ -29,21 +29,21 @@ 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); + for eqnum = 1:size(equations, 1) + phi = RK4Phi(equations{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) + for eqnum = 1:size(equations, 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; + epsilon = abs(stepval(eqnum)) * relativeEpsilon + absoluteEpsilon; eqalpha = epsilon / delta; % minimum alpha wins diff --git a/ENUME/projectC/task2.m b/ENUME/projectC/task2.m index 80a99b3a..ce908f1c 100644 --- a/ENUME/projectC/task2.m +++ b/ENUME/projectC/task2.m @@ -8,7 +8,7 @@ interval = [0; 15]; % define available algorithms algorithms = { - 'RK4', @rk4, [0.01, 0.013408]; + 'RK4', @RK4, [0.01, 0.013408]; 'Adams PC', @adamspc, [0.002, 0.01305] }; @@ -59,7 +59,7 @@ for alg = 1 : 2 end % solve ODE using RK4 with automatic step size -[result, sizes, errors] = rk4auto(sysfuncts, initvalues, interval, ... +[result, sizes, errors] = RK4Automatic(sysfuncts, initvalues, interval, ... 1e-5, 10e-10, 10e-10); % plot trajectory