Attempting to perform the Runge-Kutta-Fehlberg Algorithm

Question

Riley 2023-2-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1906105-attempting-to-perform-the-runge-kutta-fehlberg-algorithm

编辑： Torsten 2023-2-4

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

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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-2-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1906105-attempting-to-perform-the-runge-kutta-fehlberg-algorithm#answer_1163450

编辑：Torsten 2023-2-4

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Jim Riggs 2023-2-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1906105-attempting-to-perform-the-runge-kutta-fehlberg-algorithm#answer_1163400

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Attempting to perform the Runge-Kutta-Fehlberg Algorithm

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Attempting to perform the Runge-Kutta-Fehlberg Algorithm

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论