<a href="https://colab.research.google.com/github/john-s-butler-dit/Numerical-Analysis-Python/blob/master/Chapter%2002%20-%20Higher%20Order%20Methods/2_Problem%20Sheet/02_PS%20Question%202a.ipynb" target="_parent"></a>

# Problem Sheet 2 Question 2a - 3rd Order Taylor

The general form of the population growth differential equation
\begin{equation}y^{'}=t-y, \ \ (0 \leq t \leq 4),\end{equation}
with the initial condition
\begin{equation}y(0)=1.\end{equation}
For N=4
with the analytic (exact) solution
\begin{equation} y= 2e^{-t}+t+1.\end{equation}


Apply __3rd Order Taylor 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.




```python
## 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
\begin{equation} h=\frac{b-a}{N}.\end{equation}
## Specific Discrete Interval
Here the interval is $0\leq t \leq 4$ with $N=4$ 
\begin{equation} h=\frac{4-0}{4}=1.\end{equation}
This gives the 5 discrete points with stepsize h=1:
\begin{equation} t_0=0, \ t_1=1, \ ... t_{4}=4. \end{equation}
This is generalised to 
\begin{equation} t_i=0+i1, \ \ \ i=0,1,...,4.\end{equation}
The plot below illustrates the discrete time steps from 0 to 4.


```python
### Setting up time
t_end=4
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();
```

## 3rd Order Taylor Solution

The 3rd Order Taylor difference equation is given by
\begin{equation}w_{i+1}=w_i+h\left[f(t_i,w_i)+\frac{h}{2}f'(t_i,w_i)+\frac{h^2}{6}f''(t_i,w_i)\right],\end{equation}
where
\begin{equation} f(t_i,w_i)=t_i-w_i,\end{equation}
\begin{equation}f'(t_i,w_i)=1-t_i+w_i,\end{equation}
and
\begin{equation} f''(t_i,w_i)=-1+t_i-w_i,\end{equation}


which gives
\begin{equation} w_{i+1}=w_i+h(t_i-w_i)+\frac{h^2}{2}(1-t_i+w_i)+\frac{h^3}{6}(-1+t_i-w_i), \end{equation}
for $i=0,1,2,3$ with the initial condition,
\begin{equation}w_0=1.\end{equation}


```python
IC=1

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


```python
y=2*np.exp(-t)+t-1 # Exact Solution
```

## Global Error
The upper bound of the global error is:
\begin{equation} |y(t_i)-w_i|\leq\frac{Mh^3}{24}\frac{1}{L} |e^{Lt_i}-1| \end{equation}
where
$|y''''(t)|\leq M $ for $t \in (0,4)$ 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 3rd order Taylor method is
\begin{equation} |y(t_1)-w_1|<\tau=\frac{Mh^3}{24} \end{equation}

by differentiation the exact solution 
\begin{equation} y= 2e^{-t}+t+1\end{equation}
four times
\begin{equation}y'=-2e^{-t}+1,\end{equation}
\begin{equation}y''=2e^{-t},\end{equation}
\begin{equation}y'''=-2e^{-t},\end{equation}
and
\begin{equation}y''''=2e^{-t}.\end{equation}

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

\begin{equation}|y''(t)|\leq |2e^{0}|=2, \end{equation}
therefore
\begin{equation}M=2.\end{equation}
This gives the local error
\begin{equation} |y(t_1)-w_1|<\tau=\frac{Mh^3}{24}=\frac{2h^3}{24}=\frac{h^3}{12}.\end{equation}
### Lipschitz constant
The Lipschitz constant $L$ is from the Lipschitz condition,
\begin{equation}\left| \frac{\partial f(t,y)}{\partial t}\right|\leq L. \end{equation}
The constant can be found by taking partical derivative of $f(t,y)=t-y$ with respect to $y$
\begin{equation}\frac{\partial f(t,y)}{\partial y}=-1\end{equation}
\begin{equation}L=|-1|=1.$$

### Specific Global Upper Bound Formula
Puting these values into the upper bound gives:
\begin{equation} |y(t_i)-w_i|\leq\frac{h^3}{12} |e^{t_i}-1|. \end{equation}


```python
Upper_bound=h*h*h/(12)*(np.exp(t)-1) # Upper Bound
```

## Results


```python
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'=t-y$", fontsize=20)
plt.tight_layout()
plt.subplots_adjust(top=0.85)    
```


```python

d = {'time t_i': t,    ' 2nd Order Taylor (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
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>time t_i</th>
      <th>2nd Order Taylor (w_i)</th>
      <th>Exact (y_i)</th>
      <th>Exact Error( |y_i-w_i|)</th>
      <th>Upper Bound</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.0</td>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1.0</td>
      <td>0.666667</td>
      <td>0.735759</td>
      <td>0.069092</td>
      <td>0.143190</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2.0</td>
      <td>1.222222</td>
      <td>1.270671</td>
      <td>0.048448</td>
      <td>0.532421</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3.0</td>
      <td>2.074074</td>
      <td>2.099574</td>
      <td>0.025500</td>
      <td>1.590461</td>
    </tr>
    <tr>
      <th>4</th>
      <td>4.0</td>
      <td>3.024691</td>
      <td>3.036631</td>
      <td>0.011940</td>
      <td>4.466513</td>
    </tr>
  </tbody>
</table>
</div>




```python

```

088ca063c5656e2104496c06aa6b85a744259638
ipynb
Jupyter Notebook
Chapter 02 - Higher Order Methods/2_Problem Sheet/202a_Problem_Sheet.ipynb
jjcrofts77/Numerical-Analysis-Python
97e4b9274397f969810581ff95f4026f361a56a2
MIT
2019-09-05T21:39:12.000Z
2022-03-26T14:00:25.000Z
Chapter 02 - Higher Order Methods/2_Problem Sheet/202a_Problem_Sheet.ipynb
jjcrofts77/Numerical-Analysis-Python
97e4b9274397f969810581ff95f4026f361a56a2
MIT
Chapter 02 - Higher Order Methods/2_Problem Sheet/202a_Problem_Sheet.ipynb
jjcrofts77/Numerical-Analysis-Python
97e4b9274397f969810581ff95f4026f361a56a2
MIT
2021-06-17T15:34:04.000Z
2022-01-14T14:53:43.000Z
Qwen/Qwen-72B
1. YES
2. YES