6.4. Higher-Order Taylor Methods#
Remark
(Forward) Euler’s method is Taylor’s method of order one.
Assume that \(y(t)\) is the solution of the IVP (6.11) such that it has \((n+1)\) continuous derivatives. If expand \(y(t)\) in terms of its \(n\)th Taylor polynomial about \(t_{i}\) and then evaluate \(t_{i+1}\), we have [Burden and Faires, 2005],
where \(\xi_{i}\) is a number between \(t_{i}\) and \(t_{i+1}\). Now since,
It follows from (6.37) and (6.38) that,
The difference-equation method can be achieved by deleting the remainder term involving \(\xi_{i}\) in the above equation.
The following definition is Definition 5.11 from [Burden and Faires, 2005].
The following theorem isTheorem 5.12 from [Burden and Faires, 2005].
Proof: Observe that the equation \eqref{eq.5.17} can be expressed as,
where \(\xi_{i}\) is a number between \(t_{i}\) and \(t_{i+1}\). Thus, the local truncation error is
As \(y\in C^{n+1}[a,~b]\), we have \(y^{(n+1)}(t) = f^{(n)} (t,~y(t))\) bounded on \([a,b]\) and \(\tau_{i}(h) = \mathcal{O}(h^{n})\) for each \(i = 0, 1,\ldots, n\).
Example: Consider the following IVP, Apply Taylor’s method of orders two with \(n = 10\) to the initial-value problem problem
with exact solution
Solution:
For the method of order two we need the 1st, the 2nd and 3rd derivatives of
with respect to the variable \(t\). Thus,
Since \(y' = t^3+3\,t^2-y+1\),
Therefore,
and
import numpy as np
import pandas as pd
a = 0; b = 1
def Example1(N, a = a, b = b):
y_exact = lambda t: t**3 + 1;
alpha = 1
t = np.linspace(a, b, N+1)
h = (b-a)/(N)
f = lambda t, y: t**3 + 3*(t**2) - y + 1;
f1 = lambda t, y: -t**3 + 6*t + y -1;
fun = lambda t, y, h: y + h*f(t, y) + (h**2)*f1(t, y)/2
w = np.zeros(t.shape, dtype=float)
w[0] = alpha
for j in range(N):
w[j+1] = fun(t[j], w[j], h)
Table = pd.DataFrame({'t':t, 'y': y_exact(t),
'w':w, 'Absolute Error': abs(y_exact(t) - w)})
return Table
N = [2**(i) for i in range(3, 11)]
Eh_list = []
for n in N:
TB = Example1(N = n)
Eh_list.append(np.max(np.abs(TB['Absolute Error'])[1:]))
Table = pd.DataFrame(data = {'h':[(b-a)/n for n in N], 'n':N, 'Eh': Eh_list})
display(Table.style.set_properties(subset=['h'], **{'background-color': 'Lavender', 'color': 'Navy',
'border-color': 'DarkGreen'}).format({'Eh': "{:.6e}"}))
h | n | Eh | |
---|---|---|---|
0 | 0.125000 | 8 | 1.051778e-02 |
1 | 0.062500 | 16 | 2.547861e-03 |
2 | 0.031250 | 32 | 6.270429e-04 |
3 | 0.015625 | 64 | 1.555377e-04 |
4 | 0.007812 | 128 | 3.873265e-05 |
5 | 0.003906 | 256 | 9.664255e-06 |
6 | 0.001953 | 512 | 2.413705e-06 |
7 | 0.000977 | 1024 | 6.031316e-07 |