mirror of
https://github.com/kuhyx/WUT_Computer_Science.git
synced 2026-07-04 18:23:15 +02:00
feat: beautyfying task2.m and finishing code part
This commit is contained in:
parent
5b71034546
commit
245acae5b6
@ -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
|
||||
@ -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
|
||||
@ -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
|
||||
116
ENUME/projectC/task2.asv
Normal file
116
ENUME/projectC/task2.asv
Normal file
@ -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
|
||||
@ -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
|
||||
Loading…
Reference in New Issue
Block a user