Personal Website:
Scott Congreve

Download << Back OutputCode

%% Q1 Solution

%% a) Linear Oscillator
%
% We solve linear oscillator with $a = 0$, $b = 9$, $c = 10$, $\omega = 1$
% and $x_0=(2,1)^T$ using BDF2, BDF3, BDF4, Adams-Moulton 2, and ODE23
function q3_solution

figure
axis([0 10 -4 4]);
hold on;
xlabel('x_1');
ylabel('x_2');

[t,x]=bdf2(@oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'b*','DisplayName','BDF2');

[t,x]=bdf3(@oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'r^','DisplayName','BDF3');

[t,x]=bdf4(@oscillator,0,10, [2;1], .1, 0.001);
plot(t,x(:,1),'mo','DisplayName','BDF4');

[t,x]=am2(@oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'c+','DisplayName','Adams-Moulton 2');

[t,x]=ode23(@oscillator, [0,10], [2;1]);
plot(t,x(:,1),'ks','DisplayName','ode23');

legend('Location', 'NorthEast');

%% b) Linear System - Adams-Moulton 2-step
%
% We solve the stiff linear system with $x_0=(2,1)^T$
% using Adams-Moulton 2, Adams-Moulton 2 with AB2 initial guess,
% Adams-Bashfort 2, and Adams-Bashfort 3

figure
axis([0 10 -4 4]);
hold on;
xlabel('x_1');
ylabel('x_2');

[t,x]=bdf2(@linsystem,0,0.05, [2;1], .001);
plot(t,x(:,1),'b*','DisplayName','BDF2');

[t,x]=bdf3(@linsystem,0,0.05, [2;1], .001);
plot(t,x(:,1),'r^','DisplayName','BDF3');

[t,x]=bdf4(@linsystem,0,0.05, [2;1], .001, 0.001);
plot(t,x(:,1),'mo','DisplayName','BDF4');

[t,x]=am2(@linsystem,0,0.05, [2;1], .001);
plot(t,x(:,1),'g+','DisplayName','Adams-Moulton 2');

[t,x]=ode23(@linsystem, [0,0.05], [2;1]);
plot(t,x(:,1),'ks','DisplayName','ode23');

ylim([1 8]);
xlim([0 .06]);
legend('Location', 'NorthEast');

%% Convergence Analysis for BDF2, BDF3, BDF4
% We first perform convergence analysis with the default tolerance of $0.001$

tol = 0.001;
fprintf('BDF2 = ');
conv_analysis(@(field, t0, T, x0, h)bdf2_tol(field, t0, T, x0, h, tol));
fprintf('BDF3 = ');
conv_analysis(@(field, t0, T, x0, h)bdf3_tol(field, t0, T, x0, h, tol));
fprintf('BDF4 = ');
conv_analysis(@(field, t0, T, x0, h)bdf4(field, t0, T, x0, h, tol));

%%
% We now perform convergence analysis with a tolerance of $10^{-6}$
tol = 1e-6;
fprintf('BDF2 = ');
conv_analysis(@(field, t0, T, x0, h)bdf2_tol(field, t0, T, x0, h, tol));
fprintf('BDF3 = ');
conv_analysis(@(field, t0, T, x0, h)bdf3_tol(field, t0, T, x0, h, tol));
fprintf('BDF4 = ');
conv_analysis(@(field, t0, T, x0, h)bdf4(field, t0, T, x0, h, tol));

%%
% Finally we perform convergence analysis with a tolerance of $10^{-10}$
tol = 1e-10;
fprintf('BDF2 = ');
conv_analysis(@(field, t0, T, x0, h)bdf2_tol(field, t0, T, x0, h, tol));
fprintf('BDF3 = ');
conv_analysis(@(field, t0, T, x0, h)bdf3_tol(field, t0, T, x0, h, tol));
fprintf('BDF4 = ');
conv_analysis(@(field, t0, T, x0, h)bdf4(field, t0, T, x0, h, tol));

%%
% We note that the convergence rate seems to be limited by the tolerance
% for the nonlinear solve; most notable, look at the plots the method
% appears to stop converging once the error is down to roughly the
% tolerance level.

%% BDF4 Implementation

function [t,x] = bdf4(field, t0, T, x0, h, tol)

m = 4;
n = ceil((T-t0)/h);

t = t0+h*(0:n).';
x = ones(n+1,length(x0));
[to,xo] = rk_classical(field, t0, t0+(m-1)*h, x0, h);
t(1:m) = to(1:m);
x(1:m,:) = xo(1:m,:);

for i=m:n
    f_old = feval(field, t(i)+h, x(i,:).');
    x(i+1,:) = (48*x(i,:).' - 36*x(i-1,:).' + 16*x(i-2,:).' ...
        - 3*x(i-3,:).' + 12*h*f_old)/25;
    f_new = feval(field, t(i)+h, x(i+1,:).');
    while norm(f_old-f_new, 'inf') >= tol
        f_old = f_new;
        x(i+1,:) = (48*x(i,:).' - 36*x(i-1,:).' + 16*x(i-2,:).' ...
            - 3*x(i-3,:).' + 12*h*f_old)/25;
        f_new = feval(field, t(i)+h, x(i+1,:).');
    end
end

end

%% BDF3 Implementation (with Tolerance)

function [t,x] = bdf3_tol(field, t0, T, x0, h, tol)

m = 3;
n = ceil((T-t0)/h);

t = t0+h*(0:n).';
x = ones(n+1,length(x0));
[to,xo] = rk_classical(field, t0, t0+(m-1)*h, x0, h);
t(1:m) = to(1:m);
x(1:m,:) = xo(1:m,:);

for i=m:n
    f_old = feval(field, t(i)+h, x(i,:).');
    x(i+1,:) = (18*x(i,:).' - 9*x(i-1,:).' + 2*x(i-2,:).' + 6*h*f_old)/11;
    f_new = feval(field, t(i)+h, x(i+1,:).');
    while norm(f_old-f_new, 'inf') >= tol
        f_old = f_new;
        x(i+1,:) = (18*x(i,:).' - 9*x(i-1,:).' + 2*x(i-2,:).' + 6*h*f_old)/11;
        f_new = feval(field, t(i)+h, x(i+1,:).');
    end
end
end

%% BDF2 Implementation (with Tolerance)

function [t,x] = bdf2_tol(field, t0, T, x0, h, tol)

m = 2;
n = ceil((T-t0)/h);

t = t0+h*(0:n).';
x = ones(n+1,length(x0));
[to,xo] = rk_classical(field, t0, t0+(m-1)*h, x0, h);
t(1:m) = to(1:m);
x(1:m,:) = xo(1:m,:);

for i=m:n
    f_old = feval(field, t(i)+h, x(i,:).');
    x(i+1,:) = (4*x(i,:).' - x(i-1,:).' + 2*h*f_old)/3;
    f_new = feval(field, t(i)+h, x(i+1,:).');
    while norm(f_old-f_new, 'inf') >= tol
        f_old = f_new;
        x(i+1,:) = (4*x(i,:).' - x(i-1,:).' + 2*h*f_old)/3;
        f_new = feval(field, t(i)+h, x(i+1,:).');
    end
end
end

end