Figure 2.25:

Approximate solutions to \dot {x}=-x using explicit and implicit Euler methods with \Delta t=2.1, along with the exact solution x(t)=e^{-t}.

Code for Figure 2.25

Text of the GNU GPL.

main.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
% Forward Euler
lambda=complex(-1,0);
% stability limit for (-1,0): dt<2
% stability limit for (-1,5): dt<1/13~0.077
dt=2.1;
G=1+lambda*dt;
xe(1)=1;
xi(1)=1;
t(1)=0;
tend=20;
Nsteps=round(tend/dt);
%Forward Euler
for i=2:Nsteps
    xe(i)=G*xe(i-1);
%     xe(i)=exp(lambda*dt)*xe(i-1);
    t(i)=t(i-1)+dt;
end
% plot(t,real(xe),'b',t,real(x),'m')
hold on;
% pause

% Backward Euler
G=1/(1-lambda*dt);
for i=2:Nsteps;
    xi(i)=G*xi(i-1);
%    xe(i)=exp(lambda*dt)*xe(i-1);
%    t(i)=t(i-1)+dt;
end
% plot(t,real(x),'r')
%
% generate exact solution with small time steps
tlist=0:tend/1000:tend;
xex=real(exp(lambda*tlist));
plot(t,real(xe),'-o',t,real(xi),'-x',tlist,xex,'-');
xlabel('t');
ylabel('x(t)');
legend('EE','IE','Exact');

data1 = [t', xe', xi'];
data2 = [tlist', xex'];

save "TimeIntegrationExample.dat" data1 data2