6.3. First-order methods#

Let \(y_{j}\) denote approximation \(y\) at grid points for \(j = 1,2,3,\ldots,n\) (i.e. \(y_{j} = y(t_{j})\)) and some positive integer values of \(n\), we can use several methods. Consider the following IVP:

(6.11)#\[\begin{split}\begin{cases} \dfrac{dy}{dt} = f(t,y),& a \leq t \leq b, \\ y(a) = \alpha,& \end{cases}\end{split}\]

In this section, we focus on several first-order methods to approximate the solution of this IVP numerically. Here, it is assumed that the grid points are equally distributed throughout the interval \([a, b]\). Let \(n\) be a positive integer, we have, the following step size

\[\begin{equation*} h = \Delta t = \frac{b-a}{n} \end{equation*}\]

and each grid point can be identified as

\[\begin{equation*} t_{j} = a +jh,\qquad j = 0,1,2,\ldots,n. \end{equation*}\]

6.3.1. (Forward) Euler method#

Assume that \(y\in C^{2}[a,~b]\). Then, using Taylor polynomials and series from Section \ref{Taylor_Polynomials}, for \(j = 0, 1, 2, \ldots,n - 1\), we have,

(6.12)#\[\begin{equation} y(t_{j+1}) = y(t_{j}) + hy'(t_{j}) + \frac{h^2}{2} y''(\xi_{j}), \end{equation}\]

for some \(\xi_{j} \in(t_{j},~t_{j+1})\). Since \(y(t)\) satisfies the differential equation \eqref{IVP_01},

(6.13)#\[y(t_{j+1}) = y(t_{j}) + hf(t_{j}, y(t_{j})) + \underbrace{\frac{h^2}{2} y''(\xi_{j})}_{\text{Remainder Term}}.\]

Deleting the remainder term, Euler’s method generates \(w_{i} \approx y(t_i)\) for \(i = 1,~2,~3,\ldots,~n\). Therefore, Euler’s method approximates the solution of the IVP (6.11) at each grid point starting from \(t_0 = a\) where the solution is known as \(\alpha\). This process is carried forward iteratively until point \(t_{n-1}\).

Let \(w_{i} \approx y(t_i)\) for \(i = 1,~2,~3,\ldots,~n\), then the (Forward) Euler’s method is

(6.14)#\[\begin{split}\begin{cases} w_{j+1} = w_{j} + hf(t_{j},~w_{j}),\quad j = 0,~1,~2,\ldots,~n-1\\ w_{0} = \alpha. \end{cases}\end{split}\]
import pandas as pd 
import numpy as np

def ForwardEuler(f, y0, a, b, h= False, N=False):
    '''
    Parameters
    ----------
    f : function
        DESCRIPTION. the ODE y'=f(t,y(x))
    y0 : float
        the initial value.
    a : float
        DESCRIPTION. a is the left side of interval [a, b]
    b : float
        DESCRIPTION. b is the right side of interval [a, b]
    h : float, optional
        DESCRIPTION. The default is False. stepsize
    N : int, optional
        DESCRIPTION. The default is False. number of points.

    Returns
    -------
    Table : dataframe
        DESCRIPTION. a summary of the algorithm output

    '''
    if N:
        h = (b-a)/(N)
    if h:
        N = int((b-a)/h)
    t = np.linspace(a, b, N+1)
    y = np.zeros(t.shape, dtype=float)
    y[0] = y0
    # loop
    for j in range(N):
        y[j+1] = y[j] + h*f(t[j], y[j])
    Table = pd.DataFrame({'t':t, 'y':y})
    return Table

function [Table] = ForwardEuler(f, y0, a, b, N)
%{
Parameters
----------
f : function
    DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
    the initial value.
a : float
    DESCRIPTION. a is the left side of interval [a, b]
b : float
    DESCRIPTION. b is the right side of interval [a, b]
N : int
    DESCRIPTION. Number of points.

Returns
-------
Table : dataframe
    DESCRIPTION. a summary of the algorithm output

Example:
f = @(t, y) t.*exp(-t.^2)-2.*t.*y
a = 0
b = 1
N = 10
y0 = 1
%}
h = (b-a)/(N);
t = linspace(a, b, N+1)';
y = zeros(length(t),1);
y(1) = y0;
% loop
for j=1:N
    y(j+1) = y(j) + h*f(t(j), y(j));
end
Table = table(t,y);
end

Example: Consider the following IVP,

\[\begin{equation*} \begin{cases} y'+2ty=te^{-t^2},& 0 \leq t \leq 1,\\ y(0) = 1, \end{cases} \end{equation*}\]

with exact solution \(\displaystyle y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}\). Use the forward Euler method with \(n = 10\) for solving this IVP.

Solution:

For \(n = 10\): \([0,~1]\) can be discretized using

\[\begin{equation*} h = \frac{b - a}{n} = \frac{1 - 0}{10} = 0.1 \end{equation*}\]

as follows,

\[\begin{equation*} \left\{ 0,~0.1,~0.2,~0.3,~0.4,~0.5,~0.6,~0.7,~0.8,~0.9,~1.0 \right\} \end{equation*}\]

Observe that, here,

\[\begin{align*} y'= \underbrace{te^{-t^2}-2ty}_{f(t,~y)},\quad 0 \leq t \leq 1. \end{align*}\]

Now using \eqref{eq.06.02}, we get,

\[\begin{align*} w_{0} &= 1,\\ w_{1} &= w_{0} + hf(t_{0}, w_{0}) = 1 + (0.1)\,f(0,~1) = 1,\\ w_{2} &= w_{1} + hf(t_{1}, w_{1}) = 1 + (0.1)\,f(0.1,~1) = 0.9899004983,\\ \vdots &\\ w_{10} &= w_{9} + hf(t_{9}, w_{9}) = 0.6468407511 + (0.1)\,f(0.9,~0.6468407511) = 0.5704466419. \end{align*}\]

Furthermore, using the exact solution, we can calculate the exact values at these discretized points. Let \(y_{j} = y(t_{j})\) for \(0 \leq j \leq 10\) denote the exact values at the discretized points. Using

\[\begin{equation*} y\left(t_{j} \right) = \left(1 + \dfrac{t_{j}^{2}}{2}\right) e^{-t_{j}^{2}}, \quad j = 0,1,2,\ldots,10 \end{equation*}\]

we have,

\[\begin{equation*} y_{0} = 1.000000,~y_{1} = 0.995000,~y_{2} = 0.980005,~\dots,~y_{10} = 0.551819. \end{equation*}\]

Finally, the absolute error can be calculated using \(| w_{i}- y_{i}|\) for \(j = 0,1,2,\ldots,10\). The results can be summarized as follows,

import sys
sys.path.insert(0,'..')
import hd_tools as hd

import numpy as np
import pandas as pd
from hd_IVP_Algorithms import ForwardEuler  

# f(t, y(t)):
f = lambda t, y: t*np.exp(-t**(2))-2*t*y
(a, b) = (0, 1)
# the exact solution y(t)
y_exact = lambda t: (1+(t**2)/2)*np.exp(-t**2)
y0 = 1
# Table
Table = ForwardEuler(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] =  np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
                                                'border-color': 'DarkGreen'}).format({'Error': "{:.6e}"}))
  t y Exact Error
1 0.100000 1.000000 0.995000 4.999917e-03
2 0.200000 0.989900 0.980005 9.895270e-03
3 0.300000 0.969520 0.955058 1.446218e-02
4 0.400000 0.938767 0.920315 1.845169e-02
5 0.500000 0.897751 0.876151 2.160050e-02
6 0.600000 0.846916 0.823258 2.365822e-02
7 0.700000 0.787147 0.762720 2.442705e-02
8 0.800000 0.719830 0.696026 2.380419e-02
9 0.900000 0.646841 0.625026 2.181517e-02
10 1.000000 0.570447 0.551819 1.862748e-02

6.3.1.1. Error Bounds for (Forward) Euler’s Method#

The following lemma is Lemma 5.7 from [Burden and Faires, 2005].

Lemma

For all \(x \geq -1\) and any positive \(m\), we have

(6.15)#\[\begin{equation} 0 \leq (1 + x)^m \leq \exp(mx),\quad x \geq - 1,~m>0. \end{equation}\]

Proof: Using Taylor’s Theorem from \eqref{Taylor_App} with \(f(x) = e^{x}\), \(x_0 = 0\), and \(n = 1\), we have,

(6.16)#\[\begin{equation} e^x = 1 + x + \frac{1}{2}x^2 e^{\xi}, \end{equation}\]

where \(\xi\) is a number between \(x\) and zero. Therefore,

(6.17)#\[\begin{equation} 0 \leq 1 + x \leq 1 + x + \frac{1}{2}x^2e^{\xi} = e^{x}. \end{equation}\]

Since \(1 + x \geq 0\), we have

(6.18)#\[\begin{equation} 0\leq (1 + x)^m \leq (e^{x})^{m} = e^{mx}. \end{equation}\]

The following lemma is Lemma 5.8 from [Burden and Faires, 2005].

Lemma

Assume that \(\{ a_{i} \}_{i =0}^{k}\) is a sequence satisfying \( a_{0} \geq -t/s\) for some positive values \(s\) and \(t\), and

(6.19)#\[\begin{equation}\label{Lemma5.9_eq01} a_{i+1} \leq (1 + s)a_{i} + t,\quad i = 0,1,2,\ldots, k-1, \end{equation}\]

then,

(6.20)#\[\begin{equation} a_{i+1} \leq \exp((i+1)s) \left( a_{0} + \frac{t}{s} \right) - \frac{t}{s}. \end{equation}\]

Proof: Let’s fixed index \(i\) and then if follows from the above lemma that,

(6.21)#\[\begin{align} a_{i+1} &\leq (1 + s)a_{i} + t \notag \\ &\leq (1 + s)[(1 + s)a_{i-1} + t] + t = (1 + s)^{2}a_{i-1} + [1 + (1 + s)]t \notag \\ &\leq (1 + s)^3a_{i-2} + \left[ 1 + (1 + s) + (1 + s)^{2}\right]t \notag \\ &\vdots \\ &\leq (1 + s)^{i+1}a_{0} + \left[1 + (1 + s) + (1 + s)^2 +\ldots + (1 + s)^{i}\right]t. \end{align}\]

On the other hand,

(6.22)#\[\begin{equation} 1 + (1 + s) + (1 + s)^2 +\ldots + (1 + s)^{i} = \sum_{j = 1}^{i}(1 + s)^{j} \end{equation}\]

represent a geometric series with ratio \((1 + s)\). Thus,

(6.23)#\[\begin{equation} \frac{1 - (1+s)^{i+1}}{1 - (1+s)} = \frac{1}{s}\left[(1 + s)^{i+1} -1 \right]. \end{equation}\]

Therefore,

(6.24)#\[\begin{equation} a_{i+1} \leq (1+ s)^{i+1} a_{0} + \frac{(1+s)^{i+1} - 1}{s}t = (1+s)^{i+1}\left(a_{0} + \frac{t}{s} \right) - \frac{t}{s}. \end{equation}\]

Now, it follows from using Lemma \ref{lemma_5.7} with \(x = 1 + s\) that,

(6.25)#\[\begin{equation} a_{i+1} \leq \exp((i+1)s) \left( a_{0} + \frac{t}{s} \right) - \frac{t}{s}. \end{equation}\]

\end{proof}

The following theorem is Theorem 5.9 from [Burden and Faires, 2005].

Theorem

Assume that domain \(D = \{ (t, y):~a\leq t \leq b,~-\infty<y <\infty\}\) and \(f\in C(D)\). If \(f\) satisfies a Lipschitz condition on \(D\) in the variable \(y\) with constant \(L\) and assume that constant \(M\) exists such that

(6.26)#\[\begin{equation} | y''(t)| \leq M, \quad t\in [a, b], \end{equation}\]

where \(y(t)\) denotes the unique solution to the IVP (6.11).

Assume that \(w_0\), \(w_1\), \ldots , \(w_n\) are the approximations generated by Euler’s method for some positive integer \(n\). Therefore,

(6.27)#\[| y(t_{i}) - w_{i}| \leq \frac{hM}{2L}\left[ \exp(L(t_{i} - a)) -1 \right],\quad i = 0,~1,~2, \ldots ,~n.\]

Proof:

In case \hlg{\(i = 0\)}, the result is true since \(y(t_{0}) = w_{0} = \alpha\). From (6.13) that,

(6.28)#\[\begin{equation} y(t_{j+1}) = y(t_{j}) + hf(t_{j}, y(t_{j})) + \frac{h^2}{2} y''(\xi_{i}), \end{equation}\]

for \(i = 0, 1, . . . , n - 1\), and from the equations in (6.14),

(6.29)#\[\begin{equation} w_{i+1} = w_{i} + hf(t_{i},~w_{i})\quad i = 0,~1,~2,\ldots,~n-1. \end{equation}\]

Now, let \(y_{i} = y(t_{i})\) and \(y_{i+1} = y(t_{i+1})\), then

(6.30)#\[\begin{equation} y_{i+1} - w_{i+1} = y_{i} - w_{i} + hf(t_{i},~y_{i}) - hf(t_{i},~w_{i}) + \frac{h^2}{2} y''(\xi_{i}). \end{equation}\]

Therefore,

(6.31)#\[\begin{equation} |y_{i+1} - w_{i+1}| \leq |y_{i} - w_{i}| + h|f(t_{i},~y_{i}) - f(t_{i},~w_{i})| + \frac{h^2}{2} y''(\xi_{i}). \end{equation}\]

As \(f\) satisfies a Lipschits condition in the second variable, \(y\), with constant \(L\), and \(|y''(t)|\leq M\). Therefore,

(6.32)#\[\begin{equation} |y_{i+1} - w_{i+1}| \leq ( 1 + hL)|y_{i} - w_{i}| + \frac{h^2\,M}{2}. \end{equation}\]

At this point, we use Lemma \ref{lemma_5.8} and let \(s = hL\), \(t = (h^2\, M)/2\), and \(a_{j} = |y_{j} - w_{j}|\), for reach \(j = 0,1,\ldots, n\), it follows that

(6.33)#\[\begin{equation} |y_{i+1} - w_{i+1}| \leq e^{(i+1)hL} \left( |y_{0} - w_{0} | + \frac{h^2\,M}{2hL} \right) - \frac{h^2\,M}{2hL}. \end{equation}\]

As \(|y_{0} - w_{0}| = 0\) and \((i+1)h = t_{i+1} - t_{0} = t_{i+1} - a\). Therefore,

(6.34)#\[\begin{equation} | y(t_{i+1}) - w_{i+1}| \leq \frac{hM}{2L}\left[ \exp(L(t_{i+1} - a)) -1 \right],\quad i = 1,~ 2, \ldots , n-1. \end{equation}\]

and

(6.35)#\[\begin{equation} | y(t_{i}) - w_{i}| \leq \frac{hM}{2L}\left[ \exp(L(t_{i} - a)) -1 \right],\quad i = 0,~1,~2, \ldots ,~n. \end{equation}\]

Example: Consider the following IVP,

\[\begin{equation*} \begin{cases} y' = y - t^2 + 1,& 0 \leq t \leq 2,\\ y(0) = 0.5, \end{cases} \end{equation*}\]

with exact solution

\[\begin{equation*} y(t) = (t + 1)^2 - 0.5\,\exp(t). \end{equation*}\]

  • a. Use the forward Euler method with \(h = 0.2\) and solve this IVP.

  • b. Find bounds for the approximation errors and compare these bounds to the absolute errors.

Solution:

a. Discretizing \([0,~2]\) using \(h = 0.2\):

\[\begin{equation*} \left\{0,~0.2,~0.4,~0.6,~0.8,~1.0,~1.2,~1.4,~1.6,~1.8,~2.0\right\} \end{equation*}\]

Note that here \(n = 10\) as \(\displaystyle{ n = \frac{b - a}{h} = \frac{2 - 0}{0.2} = 10}\). On the other hand,

\[\begin{equation*} f(t, y) = y - t^2 + 1,\quad 0 \leq t \leq 2,\\ \end{equation*}\]

Therefore,

\[\begin{align*} w_{0} &= 0.5,\\ w_{1} &= w_{0} + hf(t_{0}, w_{0}) = 0.8000000000 ,\\ w_{2} &= w_{1} + hf(t_{1}, w_{1}) = 1.1520000000,\\ \vdots &\\ w_{10} &= w_{9} + hf(t_{9}, w_{9}) = 4.8657845043. \end{align*}\]

Similarly, the exact values at the discretized points can be calculated using the exact solution at \(t_0 = 0\), \(t_1 = 0.2\), \(t_2 = 0.4\),\ldots, \(t_8 = 1.6\), \(t_9 = 1.8\), \(t_{10} = 2.0\). The results can be summarized as follows,

import numpy as np
import pandas as pd
# f(t, y(t)):
cot = lambda t: 1/np.tan(t)
f = lambda t, y: y - t**2 + 1
(a, b) = (0, 2)
# the exact solution y(t)
y_exact = lambda t: (t + 1)**2 - 0.5* np.exp(t)
y0 = 0.5
# Table
Table = ForwardEuler(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] =  abs(Table['Exact'] - Table['y'])
display(Table.style.set_properties(subset=['t'],
                                   **{'background-color': 'PaleGreen', 'color': 'Black',
                                      'border-color': 'DarkGreen'}).format({'Error':'{:.6e}'}))
  t y Exact Error
0 0.000000 0.500000 0.500000 0.000000e+00
1 0.200000 0.800000 0.829299 2.929862e-02
2 0.400000 1.152000 1.214088 6.208765e-02
3 0.600000 1.550400 1.648941 9.854060e-02
4 0.800000 1.988480 2.127230 1.387495e-01
5 1.000000 2.458176 2.640859 1.826831e-01
6 1.200000 2.949811 3.179942 2.301303e-01
7 1.400000 3.451773 3.732400 2.806266e-01
8 1.600000 3.950128 4.283484 3.333557e-01
9 1.800000 4.428154 4.815176 3.870225e-01
10 2.000000 4.865785 5.305472 4.396874e-01

b. Note that,

\[\begin{align*} f(t, y) = y - t^2 + 1 \quad \Rightarrow \quad \frac{\partial}{\partial y} f(t, y) = 1, \end{align*}\]

for all values of \(y\). On the other hand,

\[\begin{align*} y(t) = (t + 1)^2 - 0.5\,\exp(t) \quad \Rightarrow \quad y'(t) = 2\,t-\frac{{\mathrm{e}}^t}{2}+2 \quad \Rightarrow \quad y''(t) = 2-\frac{{\mathrm{e}}^t}{2}. \end{align*}\]

The maximum \(\left|2-\dfrac{e^t}{2}\right|\) on \([0,~2]\) takes place at \(t = 2\) (why!?), and

\[\begin{align*} \max_{0\leq t \leq 2}|y''(t)| = \left|2-\frac{{\mathrm{e}}^2}{2}\right| = \dfrac{e^2}{2} - 2. \end{align*}\]

Using the inequality in the error bound (6.27) for Euler’s method with \(h = 0.2\), \(L = 1\), and \(M = \dfrac{e^2}{2} - 2\) gives

\[\begin{align*} | y_{i} - w_{i}| \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_i} - 1\right)\quad i = 0,1,2,\ldots,10. \end{align*}\]

Therefore,

\[\begin{align*} | y(0) - w_{0}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_0} - 1\right) = 0,\\ | y(0.2) - w_{1}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_1} - 1\right) = 0.0375173000,\\ | y(0.4) - w_{2}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_2} - 1\right) = 0.0833411000,\\ | y(0.6) - w_{3}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_3} - 1\right) = 0.1393100000,\\ | y(0.8) - w_{4}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_4} - 1\right) = 0.2076710000,\\ | y(1.0) - w_{5}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_5} - 1\right) = 0.2911680000,\\ | y(1.2) - w_{6}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_6} - 1\right) = 0.3931500000, \end{align*}\]
\[\begin{align*} | y(1.4) - w_{7}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_7} - 1\right) = 0.5177120000,\\ | y(1.6) - w_{8}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_8} - 1\right) = 0.6698520000,\\ | y(1.8) - w_{9}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_9} - 1\right) = 0.8556770000,\\ | y(2.0) - w_{10}| & \leq 0.1 \left( \frac{e^2}{2} - 2\right) \left(e^{t_{10}} - 1\right) = 1.0826400000. \end{align*}\]

Now, we can compare the absolute error from above table with these error bounds.

Example: In the previous example, investigate the order of convergence numerically.

Solution:

For measuring error, we use the following measurement tool:

\[\begin{align*} E_{h} = \max_{1 \leq j \leq \frac{b-a}{h}}\left| y(t_{j}) - y_{j}\right| \end{align*}\]
h = [2**(-i) for i in range(3, 12)]
Cols = ['h', 'N', 'Eh']
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)

for n in range(Table.shape[0]):
    TB = ForwardEuler(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
    Table.loc[n, 'Eh'] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
    
display(Table.style.set_properties(subset=['h', 'N'],
                                   **{'background-color': 'PaleGreen', 'color': 'Black',
                                      'border-color': 'DarkGreen'}).format({'h':'{:.4e}','Eh':'{:.4e}'}))
  h N Eh
0 1.2500e-01 16 2.9500e-01
1 6.2500e-02 32 1.5722e-01
2 3.1250e-02 64 8.1306e-02
3 1.5625e-02 128 4.1364e-02
4 7.8125e-03 256 2.0865e-02
5 3.9062e-03 512 1.0479e-02
6 1.9531e-03 1024 5.2510e-03
7 9.7656e-04 2048 2.6284e-03
8 4.8828e-04 4096 1.3150e-03 
from bokeh.plotting import show
hd.derivative_ConvergenceOrder(vecs = [Table['Eh'].values],
                            labels = ['Forward Euler method'],
                            xlabel = r"$$i$$",
                            ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}}  \right)$$",
                            title = 'Order of convergence: %s' % 'Forward Euler method',
                            legend_orientation = 'horizontal', ylim = [0.9, 1.02])

6.3.2. Backward (Implicit) Euler method#

The backward Euler’s method can be used to approximate the solution of the IVP (6.11). This method utilizes \(f(t_{j+1}, w_{j+1})\) instead of the current point \(f(t_{j}, w_{j})\):

(6.36)#\[w_{j+1} = w_{j} + hf(t_{j+1},~w_{j+1}).\]

where \(w_{0} = \alpha\) and \(w_{j+1}\) for \(j = 0, 1, 2, \ldots ,n - 1\) are approximated by solving the above equation in each step.

Observe that \(w_{j+1}\) appears on both sides of the equation (6.36) which necessitates the solution of an algebraic equation for the unknown [Atkinson, Han, and Stewart, 2011].

import pandas as pd 
import numpy as np
from sympy import symbols, solve

def BackwardEuler(f, y0, a, b, h= False, N=False):
    '''
    Parameters
    ----------
    f : function
        DESCRIPTION. the ODE y'=f(t,y(x))
    y0 : float
        the initial value.
    a : float
        DESCRIPTION. a is the left side of interval [a, b]
    b : float
        DESCRIPTION. b is the right side of interval [a, b]
    h : float, optional
        DESCRIPTION. The default is False. stepsize
    N : int, optional
        DESCRIPTION. The default is False. number of points.

    Returns
    -------
    Table : dataframe
        DESCRIPTION. a summary of the algorithm output

    '''
    if N:
        h = (b-a)/(N)
    if h:
        N = int((b-a)/h)
    t = np.linspace(a, b, N+1)
    y = np.zeros(t.shape, dtype=float)
    y[0] = y0
    b = symbols('b')
    for j in range(N):
        y[j+1] = float(solve(y[j] + h*f(t[j+1], b) - b)[0])
    Table = pd.DataFrame({'t':t, 'y':y})
    return Table

function [Table] = BackwardEuler(f, y0, a, b, N)
%{
Parameters
----------
f : function
    DESCRIPTION. the ODE y'=f(t,y(x))
y0 : float
    the initial value.
a : float
    DESCRIPTION. a is the left side of interval [a, b]
b : float
    DESCRIPTION. b is the right side of interval [a, b]
N : int
    DESCRIPTION. Number of points.

Returns
-------
Table : dataframe
    DESCRIPTION. a summary of the algorithm output

Example:
f = @(t, y) t.*exp(-t.^2)-2.*t.*y
a = 0
b = 1
N = 10
y0 = 1
%}
h = (b-a)/(N);
t = linspace(a, b, N+1)';
y = zeros(length(t),1);
y(1) = y0;
syms b
% loop
for j=1:N
    y(j+1) = vpasolve(y(j) + h*f(t(j+1), b) - b);
end
Table = table(t,y);
end

Example: Consider the following IVP,

\[\begin{align*} \begin{cases} y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\ y(0) = 1, \end{cases} \end{align*}\]

with exact solution

\[\begin{align*} y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}. \end{align*}\]

Use the backward Euler method for solving this IVP. Also, investigate the order of convergence numerically.

Solution: From one of our previous examples, we have,

\[\begin{align*} \begin{cases} y'+2ty=te^{-t^2},\quad 0 \leq t \leq 1,\\ y(0) = 1, \end{cases} \end{align*}\]

with exact solution \(\displaystyle{y\left(t \right) = \left(1 + \dfrac{t^{2}}{2}\right) e^{- t^{2}}}\). For \hlg{\(n = 10\)}: \([0,~1]\) can be discretized using

\[\begin{align*} h = \frac{b - a}{n} = \frac{1 - 0}{10} = 0.1 \end{align*}\]

as follows,

\[\begin{align*} \left\{ 0,~0.1,~0.2,~0.3,~0.4,~0.5,~0.6,~0.7,~0.8,~0.9,~1.0 \right\} \end{align*}\]

and using (6.36) and \(y_{0} = 1\), have

\[\begin{align*} w_{1} &= w_{0} + hf(t_{1}, w_{1}) \quad \Rightarrow \quad 1.0099004983374917 = 1.02\,w_{1} \quad \Rightarrow \quad w_{1} = 0.990099, \\ w_{2} &= w_{1} + hf(t_{2}, w_{2}) \quad \Rightarrow \quad 1.009314316564901 =1.04\,w_{2} \quad \Rightarrow \quad w_{2} = 0.980005,\\ \vdots \\ w_{10} &= w_{9} + hf(t_{10}, w_{10}) \quad \Rightarrow \quad 0.6430689394047624 =1.2\,w_{10} \quad \Rightarrow \quad w_{10} = 0.551819. \end{align*}\]

The results can be summarized as follows,

from hd_IVP_Algorithms import BackwardEuler  

# f(t, y(t)):
f = lambda t, y: t*np.exp(-t**(2))-2*t*y
(a, b) = (0, 1)
# the eact solution y(t)
y_exact = lambda t: (1+(t**2)/2)*np.exp(-t**2)
#
y0 = 1
# Table
Table = BackwardEuler(f = f, y0 = y0, a = a, b = b, N = 10)
Table['Exact'] = y_exact(Table['t'])
Table['Error'] =  np.abs(Table['Exact'] - Table['y'])
display(Table[1:].style.set_properties(subset=['Error'], **{'background-color': 'Lavender', 'color': 'Navy',
                                                'border-color': 'DarkGreen'}).format({'Error': "{:.4e}"}))
  t y Exact Error
1 0.100000 0.990099 0.995000 4.9016e-03
2 0.200000 0.970495 0.980005 9.5107e-03
3 0.300000 0.941427 0.955058 1.3631e-02
4 0.400000 0.903252 0.920315 1.7063e-02
5 0.500000 0.856539 0.876151 1.9612e-02
6 0.600000 0.802142 0.823258 2.1116e-02
7 0.700000 0.741251 0.762720 2.1469e-02
8 0.800000 0.675374 0.696026 2.0652e-02
9 0.900000 0.606281 0.625026 1.8745e-02
10 1.000000 0.535891 0.551819 1.5928e-02

Moreover, as for the error analysis, define,

\[\begin{align*} E_{h} = \max_{1 \leq j \leq \frac{b-a}{h}}\left| y(t_{j}) - y_{j}\right|, \end{align*}\]

then,

h = [2**(-i) for i in range(3, 8)]
Cols = ['h', 'N', 'Eh']
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table['h'] = h
Table['N'] = ((b-a)/Table['h']).astype(int)

for n in range(Table.shape[0]):
    TB = BackwardEuler(f = f, y0 = y0, a = a, b = b, h = Table['h'][n])
    Table.loc[n, 'Eh'] = np.max(np.abs(y_exact(TB['t'])[1:] - TB['y'][1:]))
    
display(Table.style.set_properties(subset=['h', 'N'],
                                   **{'background-color': 'PaleGreen', 'color': 'Black',
                                      'border-color': 'DarkGreen'}).format({'h':'{:.4e}','Eh':'{:.4e}'}))

hd.derivative_ConvergenceOrder(vecs = [Table['Eh'].values],
                            labels = ['Backward Euler method'],
                            xlabel = r"$$i$$",
                            ylabel = r"$$\ln \left( E_{h_{i}} / E_{h_{i-1}}  \right)$$",
                            title = 'Order of convergence: %s' % ' Backward Euler method',
                            legend_orientation = 'horizontal', ylim = [0.9, 1.1])
  h N Eh
0 1.2500e-01 8 2.6255e-02
1 6.2500e-02 16 1.3750e-02
2 3.1250e-02 32 7.0121e-03
3 1.5625e-02 64 3.5410e-03
4 7.8125e-03 128 1.7793e-03