Download<< Back
%% BDF 2/3-Step - Linear Oscillator
%
% We solve linear oscillator with $a = 0$, $b = 9$, $c = 10$, $\omega = 1$
% and $x_0=(2,1)^T$ using BDF2, BDF3, Adams-Moulton 2 and ode23.
a = 0;
b = 9;
c = 10;
omega = 1;
oscillator = @(t,x) [0 1; -b -a]*x + [0; c*cos(omega*t)];
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),'ro','DisplayName','BDF2');
[t,x]=bdf3(oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'b+','DisplayName','BDF3');
[t,x]=am2(oscillator,0,10, [2;1], .1);
plot(t,x(:,1),'ms','DisplayName','Adams-Moulton 2');
% compare with ode23
[t,x]=ode23(oscillator,[0,10], [2;1]);
plot(t,x(:,1),'k-','DisplayName','ode23');
legend('Location', 'NorthEast');
%% BDF 2/3 - Stiff Linear Problem
%
% We solve a stiff linear problem using BDF2, BDF3, Adams-Moulton 2 and ode23.
linsystem = @(t,x) [998 1998; -999 -1999]*x;
figure
axis([0 .1 -1 7]);
hold on;
xlabel('t');
ylabel('x_1');
[t,x]=bdf2(linsystem,0,.05, [2;1], 1e-3);
plot(t,x(:,1),'go','DisplayName','BDF2');
[t,x]=bdf3(linsystem,0,.05, [2;1], 1e-3);
plot(t,x(:,1),'b+','DisplayName','BDF3');
[t,x]=am2(linsystem,0,.05, [2;1], 1e-3);
plot(t,x(:,1),'ms','DisplayName','Adams-Moulton 2');
% comparing with ode23
options=odeset('AbsTol',1e-3,'RelTol',1e-3);
[t,x]=ode23(linsystem,[0,.05], [2;1], options);
plot(t,x(:,1),'k-','DisplayName','ode23');
legend('Location', 'NorthEast');
%% Satellite ODE - BDF2/3
%
% We solve satellite ODE using BDF2, BDF3, Adams-Moulton 2 and ode23.
figure
hold on;
xlabel('t');
ylabel('x_1');
zlabel('x_2');
xlim([0 7]);
ylim([-5 5]);
zlim([-5 5]);
T = 6.19216933131963970674;
x0 = [1.2; 0; 0; -1.04935750983031990726];
[t,x]=bdf2(@sat_ode,0,T, x0, 1e-3);
plot3(t, x(:,1), x(:,2),'ro','DisplayName','BDF2');
[t,x]=bdf3(@sat_ode,0,T, x0, 1e-3);
plot3(t, x(:,1), x(:,2),'b+','DisplayName','BDF3');
[t,x]=am2(@sat_ode,0,T, x0, 1e-3);
plot3(t, x(:,1), x(:,2),'ms','DisplayName','Adams-Moulton 2');
% compare with ode23
[t,x]=ode23(@sat_ode,[0,T], x0);
plot3(t, x(:,1), x(:,2),'k-','DisplayName','ode23');
view(3);
legend('Location', 'NorthEast');
