From 5b71034546a86a82df4305799db1324152a759f5 Mon Sep 17 00:00:00 2001 From: PolishPigeon Date: Sat, 8 Jan 2022 04:25:00 +0100 Subject: [PATCH] beautyfying adamspc.m --- ENUME/projectC/RK4Automatic.asv | 92 -------------------------- ENUME/projectC/adamspc.asv | 70 ++++++++++++++++++++ ENUME/projectC/adamspc.m | 110 +++++++++++++++++++++----------- ENUME/projectC/equationsLoop.m | 0 4 files changed, 141 insertions(+), 131 deletions(-) delete mode 100644 ENUME/projectC/RK4Automatic.asv create mode 100644 ENUME/projectC/adamspc.asv delete mode 100644 ENUME/projectC/equationsLoop.m diff --git a/ENUME/projectC/RK4Automatic.asv b/ENUME/projectC/RK4Automatic.asv deleted file mode 100644 index e74ee104..00000000 --- a/ENUME/projectC/RK4Automatic.asv +++ /dev/null @@ -1,92 +0,0 @@ -% 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/adamspc.asv b/ENUME/projectC/adamspc.asv new file mode 100644 index 00000000..a10670a6 --- /dev/null +++ b/ENUME/projectC/adamspc.asv @@ -0,0 +1,70 @@ +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/adamspc.m b/ENUME/projectC/adamspc.m index 1e4b777d..2ba4838d 100644 --- a/ENUME/projectC/adamspc.m +++ b/ENUME/projectC/adamspc.m @@ -1,47 +1,79 @@ -function x = adamspc(functs, init, interval, stepsize) +function x = adamspc(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 + 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, derivs] = RK4(functs, init, interval, stepsize, 5); + [x, derrivatives] = RK4(functions, initialValues, 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; - + % 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); - 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 + stepCount = ceil((interval(2) - interval(1)) / stepSize); +end + +function x = adamsPcLoop(x, derrivatives, explicitCoefficients, implicitCoefficients, stepCount, functions, stepSize) + for step = 6: (stepCount + 1) + [x, derrivatives] = PECE(x, derrivatives, explicitCoefficients, implicitCoefficients, step, functions, stepSize); + end +end + +function [x, derrivatives] = PECE(x, derrivatives, explicitCoefficients, implicitCoefficients, step, functions, stepSize) + % P + predictionOfX = adamsPredict(x, step, stepSize, explicitCoefficients, derrivatives); - % evaluate - for eqnum = 1:size(functs, 1) - derivs(eqnum, step) = functs{eqnum}(x(:, step)); + % E + derrivativePrediction = zeros(size(functions, 1), 1); + derrivativePrediction = adamsEvaluate(functions, predictionOfX, derrivativePrediction); + + % C + x = adamsCorrect(step, functions, stepSize, implicitCoefficients, derrivatives, derrivativePrediction, x); + + % E + derrivatives = adamsEvaluateTwo(functions, step, x, derrivatives); +end + + +function predictionOfX = adamsPredict(x, step, stepSize, explicitCoefficients, derrivatives) + % predict + 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 - - % append arguments to output - x = [interval(1):stepsize:(stepcount * stepsize); x]; end + +function derrivativePrediction = adamsEvaluate(functions, predictionOfX, derrivativePrediction) + for equationNumber = 1:size(functions, 1) + derrivativePrediction(equationNumber) = functions{equationNumber}(predictionOfX); + end +end + +function x = adamsCorrect(step, functions, stepSize, implicitCoefficients, derrivatives, derrivativePrediction, x) + 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 +end + +function derrivatives = adamsEvaluateTwo(functions, step, x, derrivatives) + for equationNumber = 1:size(functions, 1) + derrivatives(equationNumber, step) = functions{equationNumber}(x(:, step)); + end +end \ No newline at end of file diff --git a/ENUME/projectC/equationsLoop.m b/ENUME/projectC/equationsLoop.m deleted file mode 100644 index e69de29b..00000000