Attempting to perform the Runge-Kutta-Fehlberg Algorithm

10 views (last 30 days)
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)

Answers (2)

Torsten
Torsten on 3 Feb 2023
Edited: Torsten on 4 Feb 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 on 3 Feb 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.

Community Treasure Hunt

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

Start Hunting!