Download
<< Back
OutputCode
%% Q3 Solution
%% a) Logistic Method - Adams-Moulton 2-step
%
% We solve the logistic problem
%
% $$x'=(1-x)x, x(0)=2$$
%
% using Adams-Moulton 2, Adams-Moulton 2 with AB2 initial guess,
% Adams-Bashfort 2, and Adams-Bashfort 3
function q3_solution
figure;
axis([0 3 1 2]);
hold on;
xlabel('t');
ylabel('x');
[t,x]=am2(@logistic,0,3, 2, .1);
plot(t,x(:,1),'ro','DisplayName','Adams-Moulton 2');
[t,x]=am2_mod(@logistic,0,3, 2, .1);
plot(t,x(:,1),'mx','DisplayName','Adams-Moulton 2 (AB2 Init.)');
[t,x]=ab2(@logistic,0,3, 2, .1);
plot(t,x(:,1),'b+','DisplayName','Adams-Bashfort 2');
[t,x]=ab3(@logistic,0,3, 2, .1);
plot(t,x(:,1),'ks','DisplayName','Adams-Bashfort 3');
legend('Location', 'NorthEast');
%% b) Linear Oscillator - Adams-Moulton 2-step
%
% We solve linear oscillator with $a = 0$, $b = 9$, $c = 10$, $\omega = 1$
% and $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]=am2(@oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'ro','DisplayName','Adams-Moulton 2');
[t,x]=am2_mod(@oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'mx','DisplayName','Adams-Moulton 2 (AB2 Init.)');
[t,x]=ab2(@oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'b+','DisplayName','Adams-Bashfort 2');
[t,x]=ab3(@oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'ks','DisplayName','Adams-Bashfort 3');
legend('Location', 'NorthEast');
%% c) 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]=am2(@linsystem,0,.1, [2;1], .001);
plot(t,x(:,1),'ro','DisplayName','Adams-Moulton 2');
[t,x]=am2_mod(@linsystem,0,.1, [2;1], .001);
plot(t,x(:,1),'mx','DisplayName','Adams-Moulton 2 (AB2 Init.)');
[t,x]=ab2(@linsystem,0,.1, [2;1], .001);
plot(t,x(:,1),'b+','DisplayName','Adams-Bashfort 2');
[t,x]=ab3(@linsystem,0,.1, [2;1], .001);
plot(t,x(:,1),'ks','DisplayName','Adams-Bashfort 3');
ylim([-5 15]);
xlim([0 .12]);
legend('Location', 'NorthEast');
%% Convergence Analysis for Adams-Moulton
% We perform convergence analysis for 2-step Adams-Moulton
fprintf('Adams-Moulton 2-step = ');
conv_analysis(@am2);
fprintf('Adams-Moulton 2-step (AB2 Init.) = ');
conv_analysis(@am2_mod);
%%
% We perform convergence analysis for 2-step Adams-Moulton with various
% tolerances ($10^{-2}$, ..., $10^{-8}$) for the fixed point iteration
%
% *Note:* the following uses a modified version of |conv_analysis| which
% suppresses the output. This is purely for simplicity of the published
% output, so can be replaced by the normal |conv_analysis|
conv = zeros(7,2);
conv(:,1) = 10.^(-2:-1:-8);
figure;
for i=1:size(conv,1)
conv(i,2) = conv_analysis(@(field, t0, T, x0, h)...
am2_tol(field, t0, T, x0, h, conv(i,1)), Output=0, NewFigure=0);
hold on;
end
legend('tol=10^{-2}','tol=10^{-3}','tol=10^{-4}','tol=10^{-5}',...
'tol=10^{-6}','tol=10^{-7}','tol=10^{-7}',Location='SouthEast');
disp(array2table(conv,VariableNames={'Tolerance','Order of Convergence'}));
%%
% We notice that if the tolerance is too large it limits the rate of
% convergence due to the fact that as the step size decreases the
% convergence is limited by the tolerance as opposed to step size
end
%% Modified Adams-Moulton 2 to support passing tolerance as an argument
% <include>am2_tol.m</include>
