ex08ch2.m

these codes are for solving OED with matlab

function ex08ch2
% Investigate zero stability for two third order formulas.
% Apply them to y' = -y, y(0) = 1 on [0,1] with step size
% h = 1/2^i for i = 2:10 and starting values taken from
% the analytical solution. The maximum error is monitored.
% Method m = 1 is AB3, method m = 2 is a formula that is not
% zero stable. The coefficients for method m are stored in
% c(m,:) and d(m,:). The explicit formulas then have the form
%
% y_{n+1} =   [c(m,1)*y_n + c(m,2)*y_{n-1} + c(m,k)*y_{n-2)] +
%           h*[d(m,1)*f_n + d(m,2)*f_{n-1} + d(m,3)*f_(n-2)]
c = [   1    0    0 
       -3/2  3  -1/2 ];
d = [  23/12  -16/12  5/12
        3        0     0   ];
maxerror = zeros(10,2);
for m = 1:2
    for i = 2:10
        N = 2^i;
        h = 1/N;
        tn = 2*h;
        % y = [y_n, y_{n-1}, y_{n-2}]
        y = [exp(-2*h), exp(-h), exp(-0)];
        f = -y;
        for n = 3:N
            tnp1 = tn + h;
            ynp1 = dot(c(m,:),y) + h * dot(d(m,:),f);
            maxerror(i,m) = max(abs(ynp1 - exp(-tnp1)),maxerror(i,m));
            y = [ynp1, y(1:2)];
            f = -y;
            tn = tnp1;
        end
    end
end
fprintf('\nFor h = 1/2^i, the maximum errors were\n\n')
fprintf('    i     AB3      Unstable\n')
fprintf('  ---------------------------\n')
for i = 2:10
    fprintf('%5i %10.2e %10.2e\n',i,maxerror(i,:))
end