Topic 14.5: Runge Kutta Fehlberg (Maple)

Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

The Maple code for the modified Euler-Heun adaptive method is:

f := (t, y) -> -2*y + exp(-2*(t - 6)^2);
dsolve( {D(y)(t) = f( t, y(t) ), y(0) = 1} );

t[0] := 0.0;
y[0] := 1.0;
h := 0.1;
eps[abs] := 0.001;

for i from 1 do
	K[0] := f( t[i - 1], y[i - 1] );
	yt := y[i - 1] + h*K[0];

	K[1] := f( t[i - 1] + h, yt );
	zt := y[i - 1] + h*(K[0] + K[1])/2;

	s := sqrt( eps[abs] / abs( yt - zt ) );

	if s < 0.75 then
		h := max( h/2, 0.0125 );
	elif s > 1.5 then
		h := min( 2*h, 0.4 );
	end;

	t[i] := t[i - 1] + h;
	y[i] := y[i - 1] + h*K[0];

	if t[i] > 10 then
		break;
	end if;
end do:

plots[pointplot]( [seq( [t[k], y[k]], k = 1..i )] );

The Maple code for the modified Runge-Kutta-Fehlberg method is:

f := (t, y) -> -2*y + exp(-2*(t - 6)^2);
dsolve( {D(y)(t) = f( t, y(t) ), y(0) = 1} );

t[0] := 0.0;
y[0] := 1.0;
h := 0.1;
eps[abs] := 0.001;

for i from 1 do
	K[1] := f( t[i-1],             y[i-1] );
	K[2] := f( t[i-1] +   (1/4)*h, y[i-1] +   (1/4)*h*K[1] );
	K[3] := f( t[i-1] +   (3/8)*h, y[i-1] +   (3/8)*h*(        1/4*K[1] +         3/4*K[2]) );
	K[4] := f( t[i-1] + (12/13)*h, y[i-1] + (12/13)*h*(    161/169*K[1] -     600/169*K[2] +     608/169*K[3]) );
	K[5] := f( t[i-1] +         h, y[i-1] +         h*(  8341/4104*K[1] -  32832/4104*K[2] +  29440/4104*K[3] -   845/4104*K[4]) );
	K[6] := f( t[i-1] +   (1/2)*h, y[i-1] +   (1/2)*h*(-6080/10260*K[1] + 41040/10260*K[2] - 28352/10260*K[3] + 9295/10260*K[4] - 5643/10260*K[5]) );

	yt := y[i - 1] + h*(  2375*K[1] +  11264*K[3] +  10985*K[4] - 4104*K[5] )/20520;
	zt := y[i - 1] + h*( 33440*K[1] + 146432*K[3] + 142805*K[4] - 50787*K[5] + 10260*K[6] )/282150;

	s := root[4]( (eps[abs]*h)/(2*abs(yt - zt)) );

	if s < 0.75 then
		h := max( h/2, 0.025 );
		continue;
	end if;

	t[i] := t[i - 1] + h;
	y[i] := yt;

	if t[i] >= 10 then
		break;
	end if;

	if s > 1.5 then
		h := min( 2*h, 1.6 );
	end if;
end do:

plots[pointplot]( [seq( [t[k], y[k]], k = 1..i )] );