diff --git a/ENUME/projectC/adamspc.m b/ENUME/projectC/adamspc.m index 83018aa1..1e4b777d 100644 --- a/ENUME/projectC/adamspc.m +++ b/ENUME/projectC/adamspc.m @@ -1,6 +1,6 @@ function x = adamspc(functs, init, interval, stepsize) % obtain first five steps from RK4 - [x, derivs] = rk4(functs, init, interval, stepsize, 5); + [x, derivs] = RK4(functs, init, interval, stepsize, 5); x = x(2:end, :); % define coefficient tables diff --git a/ENUME/projectC/rk4.m b/ENUME/projectC/rk4.m deleted file mode 100644 index 7b6bb04e..00000000 --- a/ENUME/projectC/rk4.m +++ /dev/null @@ -1,49 +0,0 @@ -% solve ODE system using RK4 with constant step size -function [x, derivativesTable] = RK4(equations, initialValues, interval, stepSize, maxSteps) - x = initialValues; - - derivativesTable = buildDerivatiesTable(x, equations); - - % 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]; -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 deleted file mode 100644 index ba4bd559..00000000 --- a/ENUME/projectC/rk4auto.m +++ /dev/null @@ -1,61 +0,0 @@ -% automatic step size variant of RK4 -function [x, sizes, errors] = RK4Automatic(equations, initialValues, interval, initialStepSize, relativeEpsilon, absoluteEpsilon) - % set start points of output - args = interval(1); - x = initialValues; - - % initialize output plots - sizes = double.empty(); - errors = double.empty(); - - % integrate function until end of interval reached - stepsize = initialStepSize; - step = 0; - while 1 - % obtain the preceding function values - step = step + 1; - stepval = x(:, step); - - % advance output function - for eqnum = 1:size(equations, 1) - % generic single-step iteration - phi = RK4Phi(equations{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(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(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)) * relativeEpsilon + absoluteEpsilon; - 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