Personal Website:
Scott Congreve

Download<< Back

function run_gauss2
% Implicit RK: gauss2 
%
% the Butcher tableau
%   a=sqrt(3)
%         1/2-a/6   1/4      1/4-a/6       
%         1/2+a/6   1/4+a/6  1/4           
%                   1/2      1/2                    

% Logistic equation  
% x'=(a-b*x)*x ,  x(0)=2;   a=1, b=1

figure
axis([0 3 0 2]); hold on
[t,x]=ode23(@logistic,[0,3], 2);

plot(t,x(:,1));
title('Trajectory')
xlabel('t');
ylabel('x');

[t,x]=gauss2(@logistic,0,3, 2, 1/2); plot(t,x(:,1),'ro')
[t,x]=gauss2(@logistic,0,3, 2, 1/4); plot(t,x(:,1),'yx')
[t,x]=gauss2(@logistic,0,3, 2, 1/8); plot(t,x(:,1),'g+')

% Convergence analysis
fprintf('\nGauss 2 Convergence Analysis\n');
fprintf('============================\n\n');
conv_analysis(@gauss2)

% Comparing  gauss2 and rk_classical
%     both methods are of order p = 4
%
% test problem: lin1p.m

x0 = 2;
[to,xo]=ode23(@lin1p,[0,2], x0);

fprintf('\nGauss2 and RK Classical Comparison\n');
fprintf('==================================\n\n');
fprintf('  h   | Classical RK |   Gauss 2   \n');
fprintf('------+--------------+-------------\n');
for h = [2 2.7 2.75 2.8 3]/50
    figure
    axis([-.1 2 -3 3]);
    hold on;
    title(sprintf('trajectory (h=%.3f)', h));
    xlabel('time');
    ylabel('x');
    plot(to,xo,'k','DisplayName','ode23');

    tic;
    [t,x]=rk_classical(@lin1p,0,2, x0, h); 
    time_rk = toc;
    plot(t,x,'b','DisplayName','Classical RK');

    tic;
    [t,x]=gauss2(@lin1p,0,2, x0, h);
    time_gauss2 = toc;
    plot(t,x,'r','DisplayName','Gauss 2');
    legend('Location', 'NorthEast');
    fprintf('%5.3f |  %10.8f  |  %10.8f  \n', h, time_rk, time_gauss2);
end