Example 4.3

A second order equation

Contents

Setup

close all;
clear h w Tmax f T y_true Ye Te i A Yi Ti Yt Tt C D;

h = 0.15;
w = 4;
Tmax = 6;

Analyic solution

f = @(x) cos(w*x);
T = linspace(0,6,100);
y_true = f(T);

Explicit Euler

Ye = [1;0];
Te = 0;
i = 1;
A = [1 h; -w^2*h 1];
while (Te(i) <= Tmax )
    Ye(:,i+1) = A*Ye(:,i);
    Te(i+1) = Te(i)+h;
    i = i+1;
end

Implicit Euler

Yi = [1;0];
Ti = 0;
i = 1;
B = [1 -h; w^2*h 1];
while (Ti(i) <= Tmax )
    Yi(:,i+1) = B\Yi(:,i);
    Ti(i+1) = Ti(i)+h;
    i = i+1;
end

Trapezoidal

Yt = [1;0];
Tt = 0;
i = 1;
C = [1 -h/2; w^2*h/2 1];
D = [1 h/2; -w^2*h/2 1];
while (Tt(i) <= Tmax )
    Yt(:,i+1) = C\D*Yt(:,i);
    Tt(i+1) = Tt(i)+h;
    i = i+1;
end

Plot generation

figure(1);clf
plot(T,y_true,'k')
hold on
plot(Te,Ye(1,:),'r:')
plot(Ti,Yi(1,:),'g-.')
plot(Tt,Yt(1,:),'--')
axis([0 6 -2.6 2.6])
xlabel('t')
ylabel('y(t)')
legend('Exact','Explicit Euler','Implicit Euler','Trapezoidal')

Published with MATLAB® 7.5