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 x0, x1, x2, … , but also find polynomial approximations in the intervals between discrete mesh points. By joining grid points (x0, y0), (x1, y1), … 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 first course APMA0330
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
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
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).
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:
\[
y' = f(x,y ) , \quad y(x_0 ) = y_0 .
\]
A numerical method for solving the above initial value problem starts with choosing a set of mesh points
where \( x_n = x_0 + n\,h , \) n = 0, 1, …. For simplicity, we consider the set of uniformly distributed grid points starting from x0. We shall always denote by yn a number intended to approximate y(xn), the value of the exact solution of the initial value problem y' = f(x,y), y(x0) = y0 at x = xn.
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 xn ∈ [𝑎,b] :
Then the approximate numerical solution at mesh points \( \left( x_n , y_n \right) , \ n=0,1,2,\ldots , N, \) to the initial value problem \( y' = f(x,y), \ y(x_0 ) = y_0 \) can be evaluated according to the formula:
where 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
This relation shows that at a fixed point xn, the error in the approximation of the exact solution y by the Taylor algorithm of order m tends to zero like hm 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
where the increment function \( \displaystyle
\Phi_2 (h) = \frac{h}{2} \left. \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) \right\vert_{x=x_n , y= y_n} \)
is just adding the second order differential deviation to the next term in the same equation used for the first order Taylor expansion. We show how it works in the following examples.
Example:
We start with the initial value problem for the Riccati equation
\[
y' = y^2 - x^2 , \qquad y(0) = 1/2.
\]
Upon taking the mesh of points with uniform step size h = 0.1,
\( x_n = n\,h , \quad n=0,1,2,\ldots , \)
we denote by yn 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:
With this formula, we are able to make approximating increments:
\[
y_{n+1} = y_n + h \left( y^2_n - x^2_n \right) + h^2 \left( y_n^3 - x_n^2 y_n - x_n \right) .
\]
Now we ask Mathematica to perform calculations:
which gives the Taylor polynomial approximation at point x = 0.5 to be y5 = 0.619956, which is closed to the true value of 0.6177691815444183.
The results of applying Taylor's approximation of order 2 with h = 0.1
n
xn
yn
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 \) ):
Let us try to solve this problem numerically
using Taylor's algorithm of order 2. However, we check with a standard (very reliable) mathematica procedure, called NDSolve that provide a numerical solution to the given problem:
sol = NDSolve[{Derivative[1][y][x] == 1/(2*x - 3*y[x] + 5),
y[0] == 1}, y, {x, 0, 1.0}]
y[0.5] /. sol
{1.23391}
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.
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 y10 = 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:
The results of applying Taylor's approximation of order 2 with h = 0.1
n
xn
yn
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,
sol = NDSolve[{Derivative[1][u][x] == (x^2 - u[x]^2)*Sin[u[x]],
u[0] == 1}, u, {x, 0, 1.0}]
u[0.9] /. sol
Upon taking the mesh of points with uniform step size h = 0.1,
\( x_n = n\,h , \quad n=0,1,2,\ldots , \)
we denote by yn 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:
Naturally, we expect the values defined by the above Taylor's algorithm to approximate the corresponding values y(xn) 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
\[
y' = y^2 - x^2 , \qquad y(0) = 1/2.
\]
Upon taking the mesh of points with step size h = 0.1,
\( x_n = n\,h , \quad n=0,1,2,\ldots , \)
we denote by yn 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:
This gives the Taylor approximation of order three to be 0.617834, although the correct value is 0.6177691815444183. However, if you try step size 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.
The results of applying Taylor's approximation of order 3 with h = 0.1
n
xn
yn
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:
The approximation for this value is y10=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.
We can also apply the same approach from the previous example.
The results of applying Taylor's approximation of order 3 with h = 0.1
n
xn
yn
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.
The results of applying Taylor's approximation of order 3 with h = 0.1
n
xn
yn
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
To solve numerically the initial value problem
\[
y' = f(x,y) , \qquad y(x_0 ) = y_0 ,
\]
we introduce the grid points starting with the initial point: x0, x1, … , xN, where, for simplicity, the step size is assumed to be uniform. Then we express these points according to the formula xn = x0 + n h, n = 0,1, … , N, with fixed step size h. We shall always denote by yn a number intended to approximate y(xn), the value of the true solution of the given initial value problem at x = xn. If the method of Taylor expansion of order 4 is used, these values are calculated according to the following scheme:
This scheme is known as the 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
\[
y' = y^2 - x^2 , \qquad y(0) = 1/2.
\]
Upon taking the mesh of points with step size h = 0.1,
\( x_n = n\,h , \quad n=0,1,2,\ldots , \)
we denote by yn 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:
Entering 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:
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 y10 = 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
sol = NDSolve[{Derivative[1][u][x] == (x^2 - u[x]^2)*Sin[u[x]],
u[0] == 1}, u, {x, 0, 1.0}]
u[1.0] /. sol
We compare the fourth order Taylor series approximation y10 = 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
