Illustration of implementing Forward Euler, computing errors, and exploring absolute stability.
The next cell sets things up so that NumPy
and matplotlib
are imported and plots appear inline.
%pylab inline
Populating the interactive namespace from numpy and matplotlib
def euler(f, t0, tfinal, eta, nsteps, verbose=True):
un = zeros(nsteps+1)
un[0] = eta
dt = float(tfinal) / nsteps
if verbose:
print "dt = %10.6f" % dt
tn = linspace(t0,tfinal,nsteps+1)
for n in range(nsteps):
un[n+1] = un[n] + dt*f(un[n],tn[n])
return tn,un
Solve $u'(t) = -\lambda (u(t) - \cos(t)) - \sin(t)$ with $u(0) = 1$.
The solution of $u(t) = \cos(t)$ regardless of the value of $\lambda$, but the behavior of numerical methods can change dramatically.
lam = -1.
def f(u,t):
return lam * (u - cos(t)) - sin(t)
This defines the function $f(u,t)$. Note we use lam
for $\lambda$ because lambda
is a reserved keyword in Python, used to define "lambda functions", a handy way to describe a one-line function on the fly. We'll use this to define the true solution function:
utrue = lambda t: cos(t)
Initial data and final time:
t0 = 0; tfinal = 10.; eta = 1.
Solve with Forward Euler and plot the solution:
nsteps = 50
tn,un = euler(f, t0, tfinal, eta, nsteps)
plot(tn,un,'bo-')
plot(tn,utrue(tn),'k')
dt = 0.200000
[<matplotlib.lines.Line2D at 0x105ea5990>]
Try the following:
lam = -1
and nsteps = 50, 100
lam = -100
and nsteps = 100, 400, 500
lam = -100
and nsteps = 490
and zoom in with xlim(6,7)
nvals = array([25,50,100,200,400,800,1600,3200])
dtvals = tfinal / nvals
errors_euler = []
for nsteps in nvals:
tn,un = euler(f, t0, tfinal, eta, nsteps, False)
en = abs(un[-1] - utrue(tfinal))
print "nsteps = %6i, error = %10.8f" % (nsteps, en)
errors_euler.append(en)
loglog(dtvals, errors_euler, 'o-')
xlabel('Time step')
ylabel('Error at t = %s' % tfinal)
ylim(1e-4, 1)
nsteps = 25, error = 0.15545439 nsteps = 50, error = 0.07313109 nsteps = 100, error = 0.03553682 nsteps = 200, error = 0.01752463 nsteps = 400, error = 0.00870295 nsteps = 800, error = 0.00433683 nsteps = 1600, error = 0.00216478 nsteps = 3200, error = 0.00108148
(0.0001, 1)
f = lambda u,t: u**2 - sin(t) - cos(t)**2
tfinal = 20
nsteps = 160
tn,un = euler(f, t0, tfinal, eta, nsteps)
plot(tn,un,'bo-')
plot(tn,utrue(tn),'k')
dt = 0.125000
[<matplotlib.lines.Line2D at 0x108341c50>]
Try other values of nsteps = 150, 140, 130, 120, 110
.