Effacer les filtres
Effacer les filtres

Attempting to perform the Runge-Kutta-Fehlberg Algorithm

17 vues (au cours des 30 derniers jours)
Riley
Riley le 3 Fév 2023
Modifié(e) : Torsten le 4 Fév 2023
I am attempting to code the Runge-Kutta-Fehlberg Algorithm for the solution of differential equations. I feel like my code is close, however on the example problem my code runs only 4 iterations and my updated h-step size exceeds the minimum value without having computed a low enough error estimation. Does anyone see the error in my code that would lead to perhaps R, h, or del being computed incorrectly? Perhaps my K's have something incorrect? I have stared at it for hours and can't locate it:
(I know my codes probably not following all the correct coding rules, but I'm really just looking to get the right solution and don't care about optimization apart from that)
syms t y(t)
eqn = diff(y,t) == (y^2+y)/t;
cond = y(1)==-2;
sol(t) = dsolve(eqn,cond)
sol(t) = 
%%%%%%%%%%%%%%%%%%%% 2c %%%%%%%%%%%%%%%%%%%%%
a = 1;
b = 3;
j = 1;
TOL = 10.^(-4);
hmax = .5;
hmin = .02;
h = hmax;
FLAG = 1;
t = a:hmin:b;
y = zeros(size(t));
y_4 = zeros(size(t));
tempy_4 = zeros(size(t));
y_5 = zeros(size(t));
tempy_4(1) = -2;
y_4(1) = -2;
y_5(1) = -2;
n = numel(y);
i = 1;
a1 = 1 /4;
b1 = 3 /8;
b2 = 3 /32;
b3 = 9 /32;
c1 = 12 /13;
c2 = 1932 /2197;
c3 = 7200 /2197;
c4 = 7296 /2197;
d1 = 439 /216;
d2 = 8;
d3 = 3680 /513;
d4 = 845 /4104;
e1 = 1 /2;
e2 = 8 /27;
e3 = 2;
e4 = 3544 /2565;
e5 = 1859 /4104;
e6 = 11 /40;
while FLAG == 1
fprintf('while statement begun \n');
f = @(t,y) (y.^2 + y)./t;
K1 = h *f(t(i), tempy_4(i));
K2 = h *f((t(i)+a1 *h), (tempy_4(i) + a1 *K1));
K3 = h *f((t(i)+b1 *h), (tempy_4(i) + b2 *K1 + b3 *K2));
K4 = h *f((t(i)+c1 *h), (tempy_4(i) + c2 *K1 - c3 *K2 + c4 *K3));
K5 = h *f((t(i)+h), (tempy_4(i) + d1 *K1 - d2 *K2 + d3 *K3 - d4 *K4));
K6 = h *f((t(i)+e1 *h), (tempy_4(i) - e2 *K1 + e3 *K2 - e4 *K3 + e5 *K4 - e6 *K5));
tempy_4(i+1) = y_4(i) + (25 /216) *K1 + (1408 /2565) *K3+(2197 /4104) *K4-(1 /5) *K5;
y_5(i+1) = y_5(i)+(16 /135) *K1 + (6656 /12835) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;
R = abs(y_5(i+1)-tempy_4(i+1)) /h;
if R <= TOL
fprintf('Step 2 \n');
t(i+1) = t(i) + h;
y_4(i+1) = y_4(i) + 25 /216 *K1 + 1408 /2565 *K3+2197 /4104 *K4-1 /5 *K5;
i = i+1;
end
del = .84.*(TOL./R).^(.25);
if del <= .1
fprintf('Step 3.1 \n');
h = .1*h;
elseif del >= 4
fprintf('Step 3.2 \n');
h = 4.*h;
else
fprintf('Step 3.3 \n');
h = del.*h;
end
if h > hmax
fprintf('Step 4 \n');
h = hmax;
end
if t(i) >= b
fprintf('Step 5.1 \n');
FLAG = 0;
elseif t + h > b
fprintf('Step 5.2 \n');
h = b - t;
elseif h < hmin
fprintf('minimum h exceeded \n');
FLAG = 0;
end
j = j + 1;
if j > 10
FLAG = 0;
end
end
while statement begun
Step 3.3
while statement begun
Step 3.3
while statement begun
Step 3.3
while statement begun
Step 3.3
minimum h exceeded
plot(t,y_4)

Réponses (2)

Torsten
Torsten le 3 Fév 2023
Modifié(e) : Torsten le 4 Fév 2023
A partial answer is given by the code below - your test ODE is solved correctly with fixed time stepsize.
Errors in your original code that are corrected are:
h = hmin;
instead of
h = hmax;
y_5(i+1) = y_4(i)+(16 /135) *K1 + (6656 /12835) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;
instead of
y_5(i+1) = y_5(i)+(16 /135) *K1 + (6656 /12835) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;
12825
instead of
12835
in the denominator (see answer by Jim Riggs)
The adaptive version gives wrong results. The main reason is that the time vector t is no longer valid because you change h from step to step. Thus evaluating your function with the t(i) arguments is false.
I suggest you start from the code below and try to incorporate the adaptive stepsize control step by step.
syms t y(t)
eqn = diff(y,t) == (y^2+y)/t;
cond = y(1)==-2;
sol(t) = dsolve(eqn,cond)
sol(t) = 
f_ana = matlabFunction(sol);
%%%%%%%%%%%%%%%%%%%% 2c %%%%%%%%%%%%%%%%%%%%%
a = 1;
b = 3;
j = 1;
TOL = 10.^(-4);
hmax = .5;
hmin = .02;
h = hmin;
FLAG = 1;
t = a:hmin:b;
y = zeros(size(t));
y_4 = zeros(size(t));
tempy_4 = zeros(size(t));
y_5 = zeros(size(t));
tempy_4(1) = -2;
y_4(1) = -2;
y_5(1) = -2;
n = numel(y);
i = 1;
a1 = 1 /4;
b1 = 3 /8;
b2 = 3 /32;
b3 = 9 /32;
c1 = 12 /13;
c2 = 1932 /2197;
c3 = 7200 /2197;
c4 = 7296 /2197;
d1 = 439 /216;
d2 = 8;
d3 = 3680 /513;
d4 = 845 /4104;
e1 = 1 /2;
e2 = 8 /27;
e3 = 2;
e4 = 3544 /2565;
e5 = 1859 /4104;
e6 = 11 /40;
%while FLAG == 1
fprintf('while statement begun \n');
while statement begun
f = @(t,y) (y.^2 + y)./t;
for i = 1:numel(t)-1
K1 = h *f(t(i), y_4(i));
K2 = h *f((t(i)+a1 *h), (y_4(i) + a1 *K1));
K3 = h *f((t(i)+b1 *h), (y_4(i) + b2 *K1 + b3 *K2));
K4 = h *f((t(i)+c1 *h), (y_4(i) + c2 *K1 - c3 *K2 + c4 *K3));
K5 = h *f((t(i)+h), (y_4(i) + d1 *K1 - d2 *K2 + d3 *K3 - d4 *K4));
K6 = h *f((t(i)+e1 *h), (y_4(i) - e2 *K1 + e3 *K2 - e4 *K3 + e5 *K4 - e6 *K5));
y_4(i+1) = y_4(i) + (25 /216) *K1 + (1408 /2565) *K3+(2197 /4104) *K4-(1 /5) *K5;
y_5(i+1) = y_4(i)+(16 /135) *K1 + (6656 /12825) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;
%R = abs(y_5(i+1)-tempy_4(i+1)) /h;
%if R <= TOL
% fprintf('Step 2 \n');
% t(i+1) = t(i) + h;
% y_4(i+1) = y_4(i) + 25 /216 *K1 + 1408 /2565 *K3+2197 /4104 *K4-1 /5 *K5;
% i = i+1;
%end
%del = .84.*(TOL./R).^(.25);
%if del <= .1
% fprintf('Step 3.1 \n');
% h = .1*h;
%elseif del >= 4
% fprintf('Step 3.2 \n');
% h = 4.*h;
%else
% fprintf('Step 3.3 \n');
% h = del.*h;
%end
%if h > hmax
% fprintf('Step 4 \n');
% h = hmax;
%end
%if t(i) >= b
% fprintf('Step 5.1 \n');
% FLAG = 0;
%elseif t + h > b
% fprintf('Step 5.2 \n');
% h = b - t;
%elseif h < hmin
% fprintf('minimum h exceeded \n');
% FLAG = 0;
%end
%j = j + 1;
%if j > 10
% FLAG = 0;
%end
end
hold on
plot(t,y_4)
plot(t,y_5)
plot(t,f_ana(t))
hold off
grid on

Jim Riggs
Jim Riggs le 3 Fév 2023
I found the error indicated below.
Note that your notation makes your code much harder to verify because you have separated the coefficient magnitudes from their signs. For example, C3 has a magnitude of 7200/2197, and a negative sign. So I have to look at two different places to verify that C3 is implemented correctly. If all of the K equations were implemented with "+" for each term, then you put the sign on the coefficient, it is much easier to verify; C3 = -7200/2197.

Catégories

En savoir plus sur Logical dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by