# Preface

This section gives an introduction to polynomial approximations based on Taylor polynomials. In principle, it is possible to not only determine the numerical approximations for the values of the exact solution at
the mesh points *x*_{0}, *x*_{1}, *x*_{2}, … , but also find polynomial approximations in the intervals between discrete mesh points. By joining grid points (*x*_{0}, *y*_{0}), (*x*_{1}, *y*_{1}), … with straight-line segments, we are able to
obtain an overall approximation to the solution curve corresponding to the solution of the given initial-value problem.

Return to computing page for the second course APMA0340

Return to Mathematica tutorial for the second course APMA0340

Return to Mathematica tutorial for the first course APMA0330

Return to the main page for the second course APMA0340

Return to the main page for the first course APMA0330

Return to Part III of the course APMA0330

## Glossary

The Taylor series method is of general applicability, and it is the standard to which the accuracy of various other numerical methods are compared. It can be devised to have any specified degree of accuracy under some general assumptions. The most important single result in numerical computations is Taylor's theorem, which we now state below.

** Theorem (Taylor): :** Let

*f*(

*x*) have

*n*+1 continuous derivatives on the closed interval [𝑎.

*b*] for some positive integer

*n*, and let \( x, x_0 \in [a,b] . \) Then

*x*and

*x*

_{0}. ⧫

Brook Taylor (1685--1731) was educated at St. John's College of the Cambridge University (England), from which he graduated in 1709. He wrote two very significant books, ‘Methodus incrementorum directa et inversa’ and ‘Linear
Perspective,’ which were published in 1715. He was the individual to invent ‘Integration of Parts’ and also a series called the **Taylor’s Expansion**. The Taylor’s Theorem is based on the letter written by him to Machin in which he told
about the origination of this idea. However, it appears that he did not appreciate its importance and he certainly did
not bother with a proof of Taylor's theorem, which was discovered
and published in 1694 by Johann Bernoulli (1667--1748). Taylor acknowledged that his work was based on that of Newton
and Kepler. However, Johann Bernoulli claimed his priority in discovering both of his main results (Taylor's expansion and integration by parts).

Johann Bernoulli | Brook Taylor |

# Polynomial Approximations

The Taylor series method is of general applicability, and it is a standard to which we compare the accuracy of the various other numerical methods for solving an initial value problem for ordinary differential equations:

*n*= 0, 1, …. For simplicity, we consider the set of uniformly distributed grid points starting from

*x*

_{0}. We shall always denote by

*y*

_{n}a number intended to approximate

*y*(

*x*

_{n}), the value of the exact solution of the initial value problem

*y*' =

*f*(

*x,y*),

*y*(

*x*

_{0}) =

*y*

_{0}at

*x*=

*x*

_{n}. In what follows, it will be assumed that the given slope function

*f*(

*x,y*) has as many derivatives in a neighborhood of the initial point as needed.

** Theorem:** Suppose that a function

*y*(

*x*) is a solution of the initial value problem \( y' = f(x,y), \ y(x_0 ) = y_0 \) and that it has

*m+1*continuous derivatives on the interval [a,b] containing grid points \( \{ x_n \}_{n\ge 0} . \) Assume that

*y(x)*has a Taylor series expansion of order

*m*about fixed value

*x*

_{n}∈ [𝑎,

*b*] :

*h*is a step size. This scheme is known as the

**Taylor algorithm**of order

*m*. ⧫

The approximate numerical solution to the initial value problem \( y' = f(x,y ) , \quad y(x_0 ) = y_0 \) over the interval [𝑎,*b*] is derived by using the Taylor series expansion on each subinterval \( [ x_n , x_{n+1}] . \) The general step for Taylor's method of order *m* is

*m*has the property that the final global error is of the order \( O\left( h^{m+1} \right) , \) that is,

*x*

_{n}, the error in the approximation of the exact solution

*y*by the Taylor algorithm of order

*m*tends to zero like

*h*

^{m}as the stepsize

*h*tends to zero. It is possible theoretically to choose

*m*as large as necessary to make this error as small as desired. However, in practice, we usually compute two sets of approximations using step sizes

*h*and

*h/2*and compare the results.

**Theorem (Precision of Taylor's method of order m): **
Suppose that the solution to the initial value problem \( y' = f(x,y ) , \quad y(x_0 ) = y_0 \) possesses

*m+1*continuous derivatives. Let \( \left( x_n , f_n \right) , \ n=0,1,2,\ldots , N, \) be the sequence of approximations generated by Taylor's method of order

*m*. Then

Second Order Polynomial Approximations

Since the first order Taylor series approximation is identical with Euler’s method, we start with the second order one:

**Example:**We start with the initial value problem for the Riccati equation

*h*= 0.1, \( x_n = n\,h , \quad n=0,1,2,\ldots , \) we denote by

*y*

_{n}approximate values of the unknown solution at these grid points. The initial point is (0, 0.5) is obtained from the given initial condition. For a second order Taylor's approximation, we only need to determine the first derivative of the slope function:

*Mathematica*to perform calculations:

y[0] = 1/2; h = 0.1; f[x_, y_] = y^2 - x^2;

Do[y[n + 1] = y[n] + h f[h n, y[n]] + h^2 *(y[n]^3 - n*h - n^2 *h^2 *y[n]), {n, 0, 10}]

y[5]

*x*= 0.5 to be

*y*

_{5}= 0.619956, which is closed to the true value of 0.6177691815444183.

n | x_{n} |
y_{n} |
Exact solution | Absolute error |
---|---|---|---|---|

1 | 0.1 | 0.52625 | 0.525974 | -0.000276406 |

2 | 0.2 | 0.553349 | 0.552738 | -0.000610995 |

3 | 0.3 | 0.579441 | 0.578417 | -0.00102377 |

4 | 0.4 | 0.60244 | 0.600901 | -0.0015396 |

5 | 0.5 | 0.619956 | 0.617769 | -0.00218714 |

6 | 0.6 | 0.629224 | 0.626228 | -0.00299599 |

7 | 0.7 | 0.627042 | 0.623052 | -0.00399052 |

8 | 0.8 | 0.609753 | 0.604575 | -0.00517862 |

9 | 0.9 | 0.573298 | 0.566763 | -0.00653424 |

10 | 1.0 | 0.513405 | 0.505431 | -0.00797472 |

**Example:**Consider the initial value problem \( y’=1/(2x-3y+5), \quad y(0)=1 \) and use the standard

*mathematica*command

`DSolve`

. However, *mathematica*fails to provide the explicit solution because it is expressed through a special function, called the

**Lambert function**(which we denote by

*W*and which is defined from the equation \( W\,e^W = x \) ):

*mathematica*procedure, called

`NDSolve`

that provide a numerical solution to the given problem:
y[0.5] /. sol

To apply Taylor's algorithm, we pick up a fixed step length *h* = 0.1, which means that we need to do 10 iterations to reach the end point *x* = 1.

h:=0.1;

f1[x_, y_] := 1 / (2 x - 3 y + 5)

fx[x_, y_] := D[f1[a, b], a] /. {a -> x, b -> y}

fy[x_, y_] := D[f1[a, b], b] /. {a -> x, b -> y}

f2[x_, y_] := fx[x, y] + fy[x, y] * f1[x, y]

y[0] = 1;

Do[y[n + 1] = y[n] + h * f1[h n, y[n]] + h^2 * f2[h n, y[n]]/2, {n, 0, 9}]

y[10]

At each step in the sequence, the expansion requires 3 addition, 1 subtraction, 3 multiplication, and 1 division operations. Since there are 8 steps to this sequence, in order to obtain the desired approximation, the code performs 80 operations, which would be hectic work to do by hand.
The approximation is y_{10} = 1.435290196. This algorithm requires at least 22 operations per step, which means that the entire sequence requires 220 steps. However, this is not significant since *Mathematica* does all the work.

To check the answer, we write a subroutine to evaluate numerically the initial value problem for arbitrary slope function, initial conditions, and step size:

Block[{xold = x0, yold = y0, sollist = {{x0, y0}}, h},

h = N[stepsize];

Do[xnew = xold + h;

fold = f /. {x -> xold, y -> yold};

\[Phi]2 = fold + (h/2)*((D[f, x] /. {x -> xold, y -> yold}) + (D[f, y] /. {x -> xold,

y -> yold})*(f /. {x -> xold, y -> yold}));

ynew = yold + h*\[Phi]2;

sollist = Append[sollist, {xnew, ynew}];

xold = xnew;

yold = ynew, {(xn - x0)/stepsize}];

Return[sollist]]

1.23393, 0.6, 1.27678, 0.7, 1.31833, 0.8, 1.35858, 0.9, 1.39756, 1., 1.43529}

Using* Mathematica* command **NDSolve**, we check the obtained value:

n | x_{n} |
y_{n} |
Exact solution | Absolute error |
---|---|---|---|---|

1 | 0.1 | 1.04938 | 1.04937 | -4.904×10^{-6} |

2 | 0.2 | 1.09747 | 1.09746 | -8.9437×10^{-6} |

3 | 0.3 | 1.14427 | 1.14426 | -0.000012028 |

4 | 0.4 | 1.18976 | 1.18974 | -0.000014039 |

5 | 0.5 | 1.23393 | 1.23391 | -0.0000149395 |

6 | 0.6 | 1.27678 | 1.27677 | -0.0000147313 |

7 | 0.7 | 1.31833 | 1.31832 | -0.0000134095 |

8 | 0.8 | 1.35858 | 1.35857 | -0.0000110424 |

9 | 0.9 | 1.39756 | 1.39755 | -7.70324×10^{-6} |

10 | 1.0 | 1.43529 | 1.43529 | -3.49865×10^{-6} |

**Example:**Let us take a look at the initial value problem \( y’=(x^2 -y^2) \sin (y), \quad y(0)=1 , \) for which it is impossible to find the actual solution. First, we apply the numerical solver to idenify the numerical values at grid points,

u[0.9] /. sol

*h*= 0.1, \( x_n = n\,h , \quad n=0,1,2,\ldots , \) we denote by

*y*

_{n}approximate values of unknown solution at these grid points. The initial point (0, 1) is obtained from the given initial condition. For second order Taylor's approximation, we need only to determine the first derivative of the slope function:

*Mathematica*:

f2[x_,y_] = D[f[x,y],x] + D[f[x,y],y] f[x,y]

*Mathematica*to do all calculations

Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2*f2[h*n, y[n]]/2, {n, 0, 10}]

y[10]

f1[x_,y_]=D[f[x,y],x]+f[x,y]*D[f[x,y],y]

y[1] = 1; h =0.1;

a = 0; b = 1; n = Floor[(b-a)/h];

Do[x[i] = 0.0 + (i-1)*h, {i,1,n+1}]

Do[{k1=f1[x[i],y[i]],

y[i+1] = y[i] + h*f[x[i],y[i]] +h*h/2*k1}, {i,1,n}]

Do[Print[x[i], " ", y[i]],{i,1,n+1}]

n | x_{n} |
y_{n} |
Exact solution | Absolute error |
---|---|---|---|---|

1 | 0.1 | 0.925207 | 0.924517 | -0.000689672 |

2 | 0.2 | 0.865145 | 0.864125 | -0.00101984 |

3 | 0.3 | 0.817529 | 0.816377 | -0.0011522 |

4 | 0.4 | 0.779706 | -0.0011734 | |

5 | 0.5 | 0.754281 | 0.753152 | -0.0011288 |

6 | 0.6 | 0.737247 | 0.736205 | -0.00104159 |

7 | 0.7 | 0.729638 | 0.728715 | -0.000922592 |

8 | 0.8 | 0.731621 | 0.730845 | -0.000776388 |

9 | 0.9 | 0.743648 | 0.743043 | -0.000605431 |

10 | 1.0 | 0.766427 | 0.76601293 | -0.000414068 |

Third Order Polynomial Approximations

The third order Taylor approximation is obtained by adding a third order differential deviation to the equation for the second order expansion.

_{3}adds just one term to Φ

_{2}:

*y*(

*x*

_{n}) of the exact solution, particularly when step size

*h*is small.

We will show how it works in two examples that were considered previously. As usual, we use a uniform step size *h*.

**Example:**We reconsider the initial value problem for the Riccati equation

*h*= 0.1, \( x_n = n\,h , \quad n=0,1,2,\ldots , \) we denote by

*y*

_{n}approximate values of unknown solution at these grid points. The initial point (0, 0.5) is obtained from the given initial condition. In order to apply the Taylor approximation, we need to determine the derivatives of the slope function:

*Mathematica*:

op[x_,y_] = Expand[D[#,x] + f[x,y]*D[#,y]]&

phi2[x_,y_] = op[x,y][op[x,y][f[x,y]]]

Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2/2*phi1[n*h, y[n]] + 1/6*h^3*phi2[h*n, y[n]] , {n, 0, 5}] y[5]

*h*= 0.5 and make one step, the corresponding approximation becomes 0.622396. With smaller step size

*h*= 0.01, we get 0.617769, but it requires 50 steps.

Now we calculate how many arithmetic operations are need to perform the third order Taylor approximation in this particular example. To evaluate *f*(*x,y*) requirs 3 arithmetic operations (2 multiplications an d1 subruction). To evaluate `phi1`

, 7 arithmetic operations are needed, and for `phi2`

we need 13 arithmetic operations. Adding all these numbers and using 11 arithmetic operations for the algorithm, we come with 34 arithmetic operations at every step.

n | x_{n} |
y_{n} |
Exact solution | Absolute error |
---|---|---|---|---|

1 | 0.1 | 0.525979 | 0.525974 | -5.57311×10^{-6} |

2 | 0.2 | 0.552752 | 0.552738 | -0.0000139251 |

3 | 0.3 | 0.578443 | 0.578417 | -0.0000259783 |

4 | 0.4 | 0.600943 | 0.600901 | -0.0000426144 |

5 | 0.5 | 0.617834 | 0.617769 | -0.0000644693 |

6 | 0.6 | 0.626319 | 0.626228 | -0.0000915256 |

7 | 0.7 | 0.623174 | 0.623052 | -0.000122428 |

8 | 0.8 | 0.604728 | 0.604575 | -0.000153514 |

9 | 0.9 | 0.566941 | 0.566763 | -0.000177769 |

10 | 1.0 | 0.505615 | 0.505431 | -0.000184346 |

**Example:**We reconsider the IVP: \( y’=1/(2x-3y+5), \quad y(0)=1 . \) To use the third order polynomial approximation, we choose a uniform grid with constant step length

*h*= 0.1. In order to apply the Taylor approximation, we need to determine the derivatives of the slope function:

Here is the syntax for the third order Taylor approximation:

op[x_, y_] = Expand[D[#, x] + f[x, y]*D[#, y]] &

phi1[x_, y_] = op[x, y][f[x, y]]

phi2[x_, y_] = op[x, y][op[x, y][f[x, y]]]

y[0] = 1 ;h = 0.1; x[0] = 0

Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2*phi1[h*n, y[n]]/2 + h^3*phi2[h*n, y[n]]/6, {n, 0, 10}]

y[10]

_{10}=1.435283295 while the true value is 1.4352866048707.

The total number of numerical operations here is 420. It is obvious at this point why using mathematical programs such as *Mathematica* is a desired approach for such problems.

n | x_{n} |
y_{n} |
Exact solution | Absolute error |
---|---|---|---|---|

1 | 0.1 | 1.04938 | 1.04937 | -4.904×10^{-6} |

2 | 0.2 | 1.09747 | 1.09746 | -8.9437×10^{-6} |

3 | 0.3 | 1.14427 | 1.14426 | -0.000012028 |

4 | 0.4 | 1.18976 | 1.18974 | -0.0000140392 |

5 | 0.5 | 1.23393 | 1.23391 | -0.0000149395 |

6 | 0.6 | 1.27678 | 1.27677 | -0.0000147313 |

7 | 0.7 | 1.31833 | 1.31832 | -0.0000134095 |

8 | 0.8 | 1.35858 | 1.35857 | -0.0000110424 |

9 | 0.9 | 1.39756 | 1.39755 | -7.70324×10^{-6} |

10 | 1.0 | 1.43529 | 1.43529 | -3.49865×10^{-6} |

**Example:** We reconsider the IVP: \( y’=(x^2 -y^2) \sin (y), \quad y(0)=1 . \) To use the third order polynomial approximation, we choose a uniform grid with constant step length *h = 0.1*.

taylor3[f_, {x_, x0_, xn_}, {y_, y0_}, stepsize_] :=

Block[{xold = x0, yold = y0, sollist = {{x0, y0}}, h},

h = N[stepsize];

Do[xnew = xold + h;

fold = f /. {x -> xold, y -> yold};

\[Phi]2 =

fold + (h/2)*((D[f, x] /. {x -> xold,

y -> yold}) + (D[f, y] /. {x -> xold,

y -> yold})*(f /. {x -> xold, y -> yold}));

\[Alpha] = (D[f, {x, 2}] /. {x -> xold, y -> yold}) +

2*(D[f, x, y] /. {x -> xold, y -> yold})*(f /. {x -> xold,

y -> yold}) + (D[f, {y, 2}] /. {x -> xold,

y -> yold})*(f /. {x -> xold, y -> yold})^2 + (D[f,

x] /. {x -> xold, y -> yold})*(D[f, y] /. {x -> xold,

y -> yold}) + (D[f, y] /. {x -> xold,

y -> yold})^2*(f /. {x -> xold, y -> yold});

ynew = yold + h*\[Phi]2 + h^3/3!*\[Alpha];

sollist = Append[sollist, {xnew, ynew}];

xold = xnew;

yold = ynew, {(xn - x0)/stepsize}];

Return[sollist]]

{0.6, 2.36595}, {0.7, 2.73957}, {0.8, 2.9658}, {0.9, 3.07578}, {1., 3.12003}}

Another example of code (with the same output):

f1[x_,y_]=D[f[x,y],x]+f[x,y]*D[f[x,y],y]

f2[x_,y_]=D[f[x,y],x,x]+2*f[x,y]*D[f[x,y],x,y] + D[f[x,y],y,y]*f[x,y]*f[x,y] + D[f[x,y],x]*D[f[x,y],y] +f[x,y]*(D[f[x,y],y])^2

y[1] = 1; h =0.1;

a = 0; b = 1; n = Floor[(b-a)/h];

Do[x[i] = 0.0 + (i-1)*h, {i,1,n+1}]

Do[{k1=f1[x[i],y[i]],

k2=f2[x[i],y[i]],

y[i+1] = y[i] + h*f[x[i],y[i]] +h*h/2*k1 + k2/6*h^3}, {i,1,n}]

Do[Print[x[i], " ", y[i]],{i,1,n+1}]

n | x_{n} |
y_{n} |
Exact solution | Absolute error |
---|---|---|---|---|

1 | 0.1 | 0.92444 | 0.924517 | 0.0000773961 |

2 | 0.2 | 0.864011 | 0.864125 | 0.00011362 |

3 | 0.3 | 0.816248 | 0.816377 | 0.000128852 |

4 | 0.4 | 0.77957 | 0.779706 | 0.000133904 |

5 | 0.5 | 0.753018 | 0.753152 | 0.000134541 |

6 | 0.6 | 0.736072 | 0.736205 | 0.000133674 |

7 | 0.7 | 0.728582 | 0.728715 | 0.000132868 |

8 | 0.8 | 0.730712 | 0.730845 | 0.000132538 |

9 | 0.9 | 0.74291 | 0.743043 | 0.000132268 |

10 | 1.0 | 0.765882 | 0.76601293 | 0.000130676 |

Fourth Order Polynomial Approximations

*x*

_{0},

*x*

_{1}, … ,

*x*

_{N}, where, for simplicity, the step size is assumed to be uniform. Then we express these points according to the formula

*x*

_{n}=

*x*

_{0}+

*n h*,

*n*= 0,1, … ,

*N*, with fixed step size

*h*. We shall always denote by

*y*

_{n}a number intended to approximate

*y*(

*x*

_{n}), the value of the true solution of the given initial value problem at

*x*=

*x*

_{n}. If the method of Taylor expansion of order 4 is used, these values are calculated according to the following scheme:

**Taylor algorithm of order 4**. It approximates the solution of the initial value problem \( y' = f(x,y), \ y(a ) = y_0 \) over the interval [𝑎,

*b*] by evaluating derivatives

*y''*,

*y'''*, and

*y''''*and using the Taylor polynomial at each step of size

*h*.

**Example:**We reconsider the initial value problem for the Riccati equation

*h*= 0.1, \( x_n = n\,h , \quad n=0,1,2,\ldots , \) we denote by

*y*

_{n}approximate values of unknown solution at these grid points. The initial point is (0, 0.5) is obtained from the given initial condition. In order to apply the Taylor approximation, we need to determine the derivatives of the slope function

*f*(

*x,y*) =

*y*² -

*x*² up to the order 3:

*Mathematica*:

op[x_,y_] = Expand[D[#,x] + f[x,y]*D[#,y]]&

phi2[x_,y_] = op[x,y][op[x,y][f[x,y]]]]

phi3[x_, y_] = op[x, y][op[x, y][op[x, y][f[x, y]]]]

*h*.

Do[y[n + 1] = y[n] + h f[n*h, y[n]] + h^2/2*phi1[n*h, y[n]] + 1/6*h^3*phi2[h*n, y[n]]+ 1/24*h^4*phi3[h*n, y[n]], {n, 0, 5}]

`y[5]`

, we get the fourth order Taylor approximation to be
0.617772; compare with the true value: 0.6177691815444183.
■ **Example:**Consider the initial value problem \( y’=1/(2x-3y+5), \quad y(0)=1 \) and let us try to find out

*y*(1). In order to apply the Taylor approximation of order four, we need to determine the derivatives of the slope function:

We set for simplicity
a fixed step length *h* = 0.1, which means that we need to do 10 iterations to reach the end point *x* = 1.

op[x_, y_] = Expand[D[#, x] + f[x, y]*D[#, y]] &

phi1[x_, y_] = op[x, y][f[x, y]]

phi2[x_, y_] = op[x, y][op[x, y][f[x, y]]]]

phi3[x_, y_] = op[x, y][op[x, y][op[x, y][f[x, y]]]]]

]

y[0] = 1.0; h = 0.1; x[0] = 0;]

Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2*phi1[h*n, y[n]]/2 + h^3*phi2[h*n, y[n]]/6 + h^4*phi3[h*n, y[n]]/24 , {n, 0, 10}]]

y[10]

At each step in the sequence, the expansion requires 3 addition, 1 subtraction, 3 multiplication, and 1 division operations. Since there are 8 steps to this sequence, in order to obtain the desired approximation, the code performs 80 operations, which would be hectic work to do by hand.
The approximation is y_{10} = 1.435290196. This algorithm requires at least 22 operations per step, which means that the entire sequence requires 220 steps. However, this is not significant since *Mathematica* does all the work.

**Example:**Let us take a look at the initial value problem \( y’=(x^2 -y^2) \sin (y), \quad y(0)=1 , \) for which it is impossible to find the actual solution. Hence, we use a standard

*Mathematica*numerical solver

*Mathematica*:

op[x_, y_] = Expand[D[#, x] + f[x, y]*D[#, y]] & phi1[x_, y_] = op[x, y][f[x, y]] phi2[x_, y_] = op[x, y][op[x, y][f[x, y]]] phi3[x_, y_] = op[x, y][op[x, y][op[x, y][f[x, y]]]]

*Mathematica*to do all calculations

Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2/2*phi1[n*h, y[n]] + 1/6*h^3*phi2[h*n, y[n]] + 1/24*h^4*phi3[h*n, y[n]], {n, 0, 10}]

Do[Print[y[n]], {n, 10}]

*y*

_{10}= 0.766021 with the true value that gives an error 8.08804×10

^{-6}. So we see that Taylor's approximation of order 4 gives 5 correct significant digits. ■

Return to Mathematica page

Return to the main page (APMA0330)

Return to the Part 1 (Plotting)

Return to the Part 2 (First Order ODEs)

Return to the Part 3 (Numerical Methods)

Return to the Part 4 (Second and Higher Order ODEs)

Return to the Part 5 (Series and Recurrences)

Return to the Part 6 (Laplace Transform)

Return to the Part 7 (Boundary Value Pr