Problem Sheet 1 Question 2b

The general form of the population growth differential equation

(98)\[\begin{equation} y^{'}=y-t, \ \ (0 \leq t \leq 2) \end{equation}\]

with the initial condition

(99)\[\begin{equation}y(0)=2\end{equation}\]

For N=4 with the analytic (exact) solution

(100)\[\begin{equation} y= e^{t}+t+1.\end{equation}\]

Apply Euler’s Method to approximate the solution of the given initial value problems using the indicated number of time steps. Compare the approximate solution with the given exact solution, and compare the actual error with the theoretical error.

## Library
import numpy as np
import math 

%matplotlib inline
import matplotlib.pyplot as plt # side-stepping mpl backend
import matplotlib.gridspec as gridspec # subplots
import warnings
import pandas as pd

warnings.filterwarnings("ignore")

General Discrete Interval

The continuous time \(a\leq t \leq b \) is discretised into \(N\) points seperated by a constant stepsize

(101)\[\begin{equation} h=\frac{b-a}{N}.\end{equation}\]

Specific Discrete Interval

Here the interval is \(0\leq t \leq 2\) with \(N=4\)

(102)\[\begin{equation} h=\frac{2-0}{4}=0.5.\end{equation}\]

This gives the 5 discrete points with stepsize h=0.5:

(103)\[\begin{equation} t_0=0, \ t_1=0.5, \ ... t_{4}=2. \end{equation}\]

This is generalised to

(104)\[\begin{equation} t_i=0+i0.5, \ \ \ i=0,1,...,4.\end{equation}\]

The plot below illustrates the discrete time steps from 0 to 2.

### Setting up time
t_end=2
t_start=0
N=4
h=(t_end-t_start)/(N)
t=np.arange(t_start,t_end+0.01,h)
fig = plt.figure(figsize=(10,4))
plt.plot(t,0*t,'o:',color='red')
plt.plot(t[0],0*t[0],'o',color='green')


plt.title('Illustration of discrete time points for h=%s'%(h))
plt.plot();

Euler Solution

The Euler difference equation is

(105)\[\begin{equation} w_{i+1}=w_i+h(w_i-t_i), \end{equation}\]

for \(i=0,1,2,3\) with \(w_0=2.\)

IC=2

INTITIAL_CONDITION=IC
w=np.zeros(N+1)
w[0]=INTITIAL_CONDITION
for i in range (0,N):
    w[i+1]=w[i]+h*(w[i]-t[i])

y=np.exp(t)+t+1 # Exact Solution

Global Error

The upper bound of the global error is:

(106)\[\begin{equation} |y(t_i)-w_i|\leq\frac{Mh}{2L} |e^{Lt_i}-1| \end{equation}\]

where \(|y''(t)|\leq M \) for \(t \in (0,2)\) and \(L\) is the Lipschitz constant. Generally we do not have access to \(M\) or \(L\), but as we have the exact solution we can find \(M\) and \(L\).

Local Error

The global error consistent of an exponential term and the local error. The local error for the Euler method is

(107)\[\begin{equation} |y(t_1)-w_1|<\tau=\frac{Mh}{2}, \end{equation}\]

by differentiation the exact solution \( y= y= e^{t}+t+1\) twice

(108)\[\begin{equation}y'= e^{t}+1,\end{equation}\]

(109)\[\begin{equation}y''=e^{t}.\end{equation}\]

Then we choose the largest possible value of \(y''\) on the interval of interest \((0,2)\), which gives,

(110)\[\begin{equation}|y''(t)|\leq e^{2}=7.38<8, \end{equation}\]

therefore

(111)\[\begin{equation}M=8.\end{equation}\]

This gives the local error

(112)\[\begin{equation}|y(t_1)-w_1|<\tau=\frac{Mh}{2}=\frac{8h}{2}=4h.\end{equation}\]

Lipschitz constant

The Lipschitz constant \(L\) is from the Lipschitz condition,

(113)\[\begin{equation}\left| \frac{\partial f(t,y)}{\partial y}\right|\leq L. \end{equation}\]

The constant can be found by taking partical derivative of \(f(t,y)=y-t\) with respect to \(y\)

(114)\[\begin{equation}\frac{\partial f(t,y)}{\partial y}=1\end{equation}\]

(115)\[\begin{equation}L=1.\end{equation}\]

Specific Global Upper Bound Formula

Puting these values into the upper bound gives:

(116)\[\begin{equation} |y(t_i)-w_i|\leq\frac{8h}{2} |e^{t_i}-1|. \end{equation}\]

Upper_bound=8*h/(2*1)*(np.exp(t)-1) # Upper Bound

Results

fig = plt.figure(figsize=(12,5))
# --- left hand plot
ax = fig.add_subplot(1,4,1)
plt.plot(t,w,'o:',color='k')
plt.plot(t[0],w[0],'o',color='green')

#ax.legend(loc='best')
plt.title('Numerical Solution (w) h=%s'%(h))

# --- right hand plot
ax = fig.add_subplot(1,4,2)
plt.plot(t,y,color='blue')
plt.title('Analytic Solution (y)')

#ax.legend(loc='best')
ax = fig.add_subplot(1,4,3)
plt.plot(t,np.abs(w-y),'o:',color='red')
plt.title('Exact Error |w-y|')

# --- title, explanatory text and save


ax = fig.add_subplot(1,4,4)
plt.plot(t,np.abs(w-y),'o:',color='red')
plt.plot(t,Upper_bound,color='purple')
plt.title('Upper Bound')

# --- title, explanatory text and save
fig.suptitle(r"$y'=y-t$", fontsize=20)
plt.tight_layout()
plt.subplots_adjust(top=0.85)    

d = {'time t_i': t,    'Euler (w_i) ':w,    'Exact (y_i) ':y,    'Exact Error( |y_i-w_i|) ':np.abs(y-w),'Upper Bound':Upper_bound}
df = pd.DataFrame(data=d)
df

	time t_i	Euler (w_i)	Exact (y_i)	Exact Error( |y_i-w_i|)	Upper Bound
0	0.0	2.0000	2.000000	0.000000	0.000000
1	0.5	3.0000	3.148721	0.148721	1.297443
2	1.0	4.2500	4.718282	0.468282	3.436564
3	1.5	5.8750	6.981689	1.106689	6.963378
4	2.0	8.0625	10.389056	2.326556	12.778112