![Contents](../../images/content.png)
![Previous Chapter](../../images/backback.png)
![Start of Chapter](../../images/backstop.png)
![Previous Topic](../../images/back.png)
![Next Topic](../../images/fwd.png)
![Next Chapter](../../images/fwdfwd.png)
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
Copyright ©2006 by Douglas Wilhelm Harder. All rights reserved.