I am implementing Runge-Kutta method with adaptive step in matlab. I get different results as compared to matlab's own ode45 and my own implementation of Runge-Kutta method with fixed step. What am I doing wrong in my code? Is it possible?
function [ result ] = rk4_modh( f, int, init, h, h_min )
%
% f - function handle
% int - interval - pair (x_min, x_max)
% init - initial conditions - pair (y1(0),y2(0))
% h_min - lower limit for h (step length)
% h - initial step length
% x - independent variable ( for example time )
% y - dependent variable - vertical vector - in our case ( y1, y2 )
function [ k1, k2, k3, k4, ka, y ] = iteration( f, h, x, y )
% core functionality performed within loop
k1 = h * f(x,y);
k2 = h * f(x+h/2, y+k1/2);
k3 = h * f(x+h/2, y+k2/2);
k4 = h * f(x+h, y+k3);
ka = (k1 + 2*k2 + 2*k3 + k4)/6;
y = y + ka;
end
% constants
% relative error
eW = 1e-10;
% absolute error
eB = 1e-10;
s = 0.9;
b = 5;
% initialization
i = 1;
x = int(1);
y = init;
while true
hy = y; hx = x;
%algorithm
[ k1, k2, k3, k4, ka, y ] = iteration( f, h, x, y );
% error estimation
for j=1:2
[ hk1, hk2, hk3, hk4, hka, hy ] = iteration( f, h/2, hx, hy );
hx = hx + h/2;
end
err(:,i) = abs(hy - y);
% step adjustment
e = abs( hy ) * eW + eB;
a = min( e ./ err(:,i) )^(0.2);
mul = a * s;
if mul >= 1
% step length admitted
keepH(i) = h;
k(:,:,i) = [ k1, k2, k3, k4, ka ];
previous(i,:) = [ x+h, y' ]; %'
i = i + 1;
if floor( x + h + eB ) == int(2)
break;
else
h = min( [mul*h, b*h, int(2)-x] );
x = x + keepH(i-1);
end
else
% step length requires further adjustments
h = mul * h;
if ( h < h_min )
error('Computation with given precision impossible');
end
end
end
result = struct( 'val', previous, 'k', k, 'err', err, 'h', keepH );
end
The function in question is:
function [ res ] = fun( x, y )
%
res(1) = y(2) + y(1) * ( 0.9 - y(1)^2 - y(2)^2 );
res(2) = -y(1) + y(2) * ( 0.9 - y(1)^2 - y(2)^2 );
res = res'; %'
end
The call is:
res = rk4( @fun, [0,20], [0.001; 0.001], 0.008 );
The resulting plot for x1 :
The result of ode45( @fun, [0, 20], [0.001, 0.001] ) is:
