Consider the ODE problem
y′y(0)=−100y,t≥0=1 whose exact solution is
y(t)=exp(−100t) import numpy as np
from matplotlib import pyplot as plt
Right hand side function and exact solution
def f(t,y):
return -100.0*y
def yexact(t):
return np.exp(-100.0*t)
t0, y0, T = 0.0, 1.0, 0.2
Forward Euler
This implements Euler method
yn=yn−1+hf(tn−1,yn−1)=(1−100h)yn−1 def feuler(t0,T,y0,h):
N = int((T-t0)/h)
y = np.zeros(N)
t = np.zeros(N)
y[0] = y0
t[0] = t0
for n in range(1,N):
y[n] = y[n-1] + h*f(t[n-1],y[n-1])
t[n] = t[n-1] + h
te = np.linspace(t0,T,100)
ye = yexact(te)
plt.plot(t,y,te,ye,'--')
plt.legend(('Numerical','Exact'))
plt.xlabel('t')
plt.ylabel('y')
plt.title('FE, step size = ' + str(h));
We first try a step size slightly larger than the stability limit.
h = 1.1*2.0/100
feuler(t0,T,y0,h)
Now we try the maximum allowed time step.
h = 2.0/100
feuler(t0,T,y0,h)
The solution is oscillating between +1 and -1. Next we reduce h by a small amount.
h = 0.95*2.0/100
feuler(t0,T,y0,h)
Now the solution is decaying but in an oscillatory manner. To get monotonic decay, we need to use h<1001.
h = 0.95*1.0/100
feuler(t0,T,y0,h)
Backward Euler
The scheme is given by
yn=yn−1+hfn,yn=1+100hyn−1 def beuler(t0,T,y0,h):
N = int((T-t0)/h)
y = np.zeros(N)
t = np.zeros(N)
y[0] = y0
t[0] = t0
for n in range(1,N):
y[n] = y[n-1]/(1.0+100*h)
t[n] = t[n-1] + h
te = np.linspace(t0,T,100)
ye = yexact(te)
plt.plot(t,y,te,ye,'--')
plt.legend(('Numerical','Exact'))
plt.xlabel('t')
plt.ylabel('y')
plt.title('BE, step size = ' + str(h));
This scheme should be stable for any step size.
h = 1.1*2.0/100
beuler(t0,T,y0,h)
h = 0.95*1.0/100
beuler(t0,T,y0,h)
Trapezoidal method
The method is given by
yn=yn−1+2h(fn−1+fn)=yn−1+2hλ(yn−1+yn) yn=1−2hλ1+2hλyn−1 def cn(t0,T,y0,h):
N = int((T-t0)/h)
y = np.zeros(N)
t = np.zeros(N)
lam = -100.0
coef = (1.0 + 0.5*h*lam)/(1.0 - 0.5*h*lam)
y[0] = y0
t[0] = t0
for n in range(1,N):
y[n] = coef*y[n-1]
t[n] = t[n-1] + h
te = np.linspace(t0,T,100)
ye = yexact(te)
plt.plot(t,y,te,ye,'--')
plt.legend(('Numerical','Exact'))
plt.xlabel('t')
plt.ylabel('y')
plt.title('TR, step size = ' + str(h));
This scheme should be stable for any time step.
h = 1.1*2.0/100
cn(t0,T,y0,h)
h = 0.95*2.0/100
cn(t0,T,y0,h)