Topic 14.5: Runge Kutta Fehlberg (Matlab)

Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

The following is an implementation in Matlab where Heun's method is used to approximate the error for Euler's method. The scaling factor is restricted to halving or doubling the size of the interval.

function [ts, ys] = heuler( f, tint, y0, h, eps_step )
    % tint is a vector of two entries represting the interval [t_0, t_end]
    t0 = tint(1);
    tend = tint(2);

    % We will start with vectors with one element and
    % continue to add further elements
    ts = t0;
    ys = y0;

    k = 1;

    while ( ts(k) < tend )
        K0 = f(ts(k), ys(k));
        K1 = f(ts(k) + h, ys(k) + h*K0);

	% Euler's approximation
        y1 = ys(k) + h*K0;

	% Heun's approximation
        z1 = ys(k) + h*(K0 + K1)/2;

        s = sqrt( eps_step/abs( y1 - z1 ) );

        if ( s < 1 )
            % h is too large
            h = h/2;
        elseif ( s < 2 )
            % h is appropriate
            ts( k + 1 ) = ts(k) + h;
            ys( k + 1 ) = y1;
            k = k + 1;
        else % s >= 2
            % h is too small
            ts( k + 1 ) = ts(k) + h;
            ys( k + 1 ) = y1;
            k = k + 1;
            h = 2*h;
        end

	% If ts(k) + h > tend, set h = tend - ts(k)
        h = min( h, tend - ts(k) );
    end
end