Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to computing page for the fourth course APMA0360
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to Mathematica tutorial for the fourth course APMA0360
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
Return to the main page for the course APMA0360
Return to Part III of the course APMA0330


This chapter gives an introduction to numerical methods used in solving ordinary differential equations of first order. Its purpose is not to give a comprehensive treatment of this huge topic, but to acquaint the reader with the basic ideas and techniques used in numerical solutions of differential equations. It also tries to dissolve a myth that Mathematica is not suitable for numerical calculations with numerous examples where iteration is an essential part of algorithms.

Introduction to Numerical Methods

Ordinary differential equations are at the heart of our perception of the physical universe. Up to this point, we have discussed analytical methods for solving differential equations that involve finding an exact expression for the solution. It turns out that the majority of engineering and applied science problems cannot be solved by methods discussed so far or when they can be solved, the formulas become so complicated and cumbersome that it is very difficult to interpret them. For this reason development of numerical methods for their solutions is one of the oldest and most successful areas of numerical computations. A universal numerical method for solving differential equations does not exist. Every numerical solver can be flanked by some problem involving differential equation. Hence, it is important to understand the limitations of every numerical algorithm and apply it wisely. This chapter gives a friendly introduction to this topic and addresses basic ideas and techniques for their realizations. It cannot be considered as a part of the comprehensive numerical analysis course and the reader is advised to study other sources to pursue a deeper knowledge.

Actually, this chapter consists of two parts. The first six sections are devoted to solving nonlinear algebraic and transcendental equations that have many applications in all fields of engineering and it is one of the oldest and most basic problems in mathematics. They also serve as an introduction to the iterative methods that are widely used in all areas of numerical analysis and applied mathematics. Mathematica has built-in functions Find Root, NSolve, Solve, Root that can solve the root-finding problems. However, programming the algorithms outlined below has great educational value. Besides, if you know what is under the hood, you can use it more efficiently or avoid misuse. We present several algorithms that can be and some will be extended for solving differential equations numerically that constitute the second part of this chapter. Since solving nonlinear equations is one of the most important problems that has many applications in all fields of engineering and computer science, we start this chapter with an introduction to this topic. Note that we do not discuss specialized algorithms for finding roots of polynomials or solving systems of algebraic equations.

Example 1: Consider a sphere of solid material floating in liquid (or water). Archimedes’ principle states that the buoyancy force is equal to the weight of the replaced liquid. Let V s = (4/3)πr 3 be the volume of the sphere of radius r, and let V w be the volume of water it displaces when it is partially submerged. In static equilibrium, the weight of the sphere is balanced by the buoyancy force
\[ \rho_s g\,V_s = \rho_w g\,V_w , \]
where ρs is the density of the sphere material, g is the acceleration due to gravity, ρw is the density of water. The volume Vh of water displaced when a sphere is submerged to a depth h is
\[ V_h = \frac{\pi}{3}\,h^2 \left( 3r-h \right) . \]

circle = Graphics[{Black, Thick, Circle[{0, 0}, 1.0]}];
line = Graphics[{Black, Thick, Line[{{-1.1, -0.3}, {1.2, -0.3}}]}];
line2 = Graphics[{Black, Line[{{0, -1}, {1.2, -1}}]}];
w = Plot[{-0.3, -1.3}, {x, -1.13, 1.13}, Filling -> {1 -> {2}}, FillingStyle -> Cyan];
ar = Graphics[{Black, Arrowheads[{{-0.04}, {0.04}}], Arrow[{{1.25, -1}, {1.25, -0.3}}]}];
th = Graphics[{Black, Text[Style["h", 18], {1.35, -0.7}]}];
fill = Graphics[{LightBlue, Polygon[Table[{0.99 Cos[theta], 0.99 Sin[theta]}, {theta, Pi + 0.3, 2 Pi - 0.29, 0.05}]]}];
Show[line, line2, w, ar, th, circle, fill]

Floating sphere

Applying Archimedes’ principle produces the following equation in term of h

\[ h^3 - 3r\,h^2 + 4\rho r^3 = 0. \tag{1.1} \]
Given values of r and the specific gravity of the sphere material ρ = ρsw, the solutions h of Eq.(1.1) can be obtained by a numerical method because the analytical formula is very complicated.
End of Example 1

Example 2: Let us consider the mechanical system represented by the four rigid rods a₁, a₂, a₃, and a₄, presented in the figure below:

rods = Graphics[{Thick, Black, Line[{{0, 0}, {0.7, 0.7}, {1.1, 0.9}, {1.5, 0}}]}]; xa = Graphics[{Thick, Black, Arrowheads[0.08], Arrow[{{-0.2, 0}, {1.95, 0}}]}]; ya = Graphics[{Thick, Black, Arrowheads[0.07], Arrow[{{0.0, 0}, {0, 1.25}}]}]; ca = Graphics[{Thick, Red, Circle[{0, 0}, 0.4, {0, Pi/4}]}]; caa = Graphics[{Thick, Red, Arrowheads[0.04], Arrow[{{0.32, 0.23}, {0.28, 0.28}}]}]; cb = Graphics[{Red, Thick, Circle[{1.5, 0}, 0.28, {0, 1.98}]}]; cbb = Graphics[{Red, Thick, Arrowheads[0.04], Arrow[{{1.45, 0.273}, {1.4, 0.26}}]}]; ta = Graphics[{Red, Text[Style["\[Alpha]", 16], {0.45, 0.17}]}]; tb = Graphics[{Red, Text[Style["\[Beta]", 16], {1.6, 0.13}]}]; c1 = Graphics[{Black, Circle[{0, 0}, 0.05]}]; c2 = Graphics[{Black, Circle[{0.7, 0.7}, 0.05]}]; c3 = Graphics[{Black, Circle[{1.1, 0.9}, 0.05]}]; c4 = Graphics[{Black, Circle[{1.5, 0}, 0.05]}]; a1 = Graphics[{Black, Line[{{-0.03, 0.04}, {0.66, 0.73}}]}]; a11 = Graphics[{Black, Line[{{0.03, -0.035}, {0.72, 0.65}}]}]; a2 = Graphics[{Black, Line[{{0.67, 0.74}, {1.08, 0.945}}]}]; a22 = Graphics[{Black, Line[{{0.73, 0.66}, {1.095, 0.84}}]}]; a3 = Graphics[{Black, Line[{{1.155, 0.91}, {1.55, 0.03}}]}]; a33 = Graphics[{Black, Line[{{1.07, 0.85}, {1.45, 0.0}}]}]; a4 = Graphics[{Black, Line[{{-0.03, 0.046}, {1.5, 0.046}}]}]; a44 = Graphics[{Black, Line[{{0.0, -0.046}, {1.5, -0.046}}]}]; ta1 = Graphics[{Black, Text[Style[Subscript[a, 1], 16], {0.3, 0.5}]}]; ta2 = Graphics[{Black, Text[Style[Subscript[a, 2], 16], {0.8, 0.95}]}]; ta3 = Graphics[{Black, Text[Style[Subscript[a, 3], 16], {1.5, 0.45}]}]; ta4 = Graphics[{Black, Text[Style[Subscript[a, 4], 16], {0.8, -0.13}]}]; tx = Graphics[{Black, Text[Style["x", 16], {1.85, -0.08}]}]; ty = Graphics[{Black, Text[Style["y", 16], {0.1, 1.08}]}]; Show[rods, xa, ya, ca, caa, cb, cbb, ta, tb, c1, c2, c3, c4, a1, a11, \ a2, a22, a3, a33, a4, a44, ta1, ta2, ta3, ta4, tx, ty]

System with four rods

For any admissible value of the angle β, let us determine the value of the corresponding angle α between the rods a₄ and a₁. The vector identity holds:

\[ {\bf a}_4 = {\bf a}_1 + {\bf a}_2 + {\bf a}_3 . \]

Since the rod a₄ is always aligned horizontally, we can deduce the following relationship between β and α:

\[ \frac{a_4}{a_1}\,\cos\beta - \frac{a_4}{a_1}\,\cos \alpha - \cos \left( \beta - \alpha \right) + \frac{a_4^1 + a_4^2 - a_2^2 - a_3^2}{2 a_1 a_3} = 0 , \tag{2.1} \]
where 𝑎i, i = 1, 2, 3, 4, is the known length of the i-th rod. This is called the Freudenstein equation (named after an American physicist and engineer Ferdinand Freudenstein (1926 – 2006) who is known as the "Father of Modern Kinematics."), and we can rewrite it as f(α) = 0, wjere
\[ f(\alpha ) = \frac{a_4}{a_1}\,\cos\beta - \frac{a_4}{a_1}\,\cos \alpha - \cos \left( \beta - \alpha \right) + \frac{a_4^1 + a_4^2 - a_2^2 - a_3^2}{2 a_1 a_3} \tag{2.2} \]
A solution in explicit form is available only for special values of β. We would also like to mention that a solution does not exist for all values of β, and may not even be unique. To solve the equation for any given β lying between 0 and π we should invoke numerical methods.
End of Example 2

Let f be a real single-valued continuously differentiable function of a real variable. If f(α) = 0, then α is said to be a zero of f or, equivalently, a root of the equation f(x) = 0. In the first part of this chapter, we present several algorithms for finding roots of such equations iteratively:

\[ x_{k+1} = \phi \left( x_k \right) , \qquad k=0,1,2,\ldots , \]
where xk is an approximation to the zero and ϕ is an iteration function. The iterative method requires the knowledge of an initial guess x0 (there may be more initial guesses for some algorithms). An initial guess x0 can usually be found by using the context in which the problem first arose or plotting a simple graph of y = f(x) will often suffice for estimating x0. The function ϕ may depend on derivatives of f up to order r in order to advance the iterative method to order r + 1.

Let { xn } be a sequence generated by an iteration function for solving the equation f(x) = 0. Let α be its root, f(α) = 0, and let εn = xn - α be the error at the n-th stage. If there exists a real number p and a positive constant Cp such that
\[ \left\vert \varepsilon_{n+1} \right\vert \le C \left\vert \varepsilon_{n} \right\vert^p , \]
then p is called the order at which the sequence { xn } converges to the root α of the equation f(x) = 0. The least possible value of C is called the asymptotic error constant.
Since the root α is unknown, the definition of the order above has only theoretical value. In practical implementation, we will mostly use the computational order of convergence (COC for short)
\[ r_c = \frac{\ln \left\vert f\left( x_k \right) / f\left( x_{k-1} \right) \right\vert}{\ln \left\vert f\left( x_{k-1} \right) / f\left( x_{k-2} \right) \right\vert} . \]
To illustrate the concept of an iterative method, you may consider the problem of finding the reciprocal of a number 𝑎 ≠ 0 without a machine division. As you will see shortly, finding a root of the function \( \displaystyle f(x) = a - 1/x = 0 \) using Newton's method leads to the recurrence \( \displaystyle x_{k+1} = x_k \left( 2 - a\,x_k \right) , \quad k=0,1,2,\ldots . \) When you make some numerical experiments, it will take only a few iterations to find the reciprocal of 𝑎 with a reasonable number of correct decimal places. In each stage, this algorithm makes only two multiplications and one subtraction.

Now we give a basic classification of iterative methods following J.F. Traub (1932--2015).

The following existence and uniqueness statements are based on the famous intermediate value theorem.

Intermediate Value Theorem: : Let f(x) be continuous on the finite closed interval [𝑎,b] and define

\[ m = \mbox{Infimum}_{a\le x \le b} \,f(x) , \qquad M = \mbox{Supremum}_{a\le x \le b} \,f(x) . \]
Then for any number γ on the interval [m,M], there is at least one point in the closed interval [𝑎,b] for which f(ζ) = γ. In particular, there are points \( \displaystyle \underline{x} \) and \( \displaystyle \overline{x} \) in [𝑎,b] for which
\[ f \left(\underline{x} \right) = m \qquad \mbox{and} \qquad f \left( \overline{x} \right) = M . \]

Theorem: : Suppose that f(x) ∈ C¹[𝑎,b] and \( \displaystyle f' (x) - p\,f(x) \ne 0 , \) where p is a real number, then the equation f(x) = 0 has at most one root in [𝑎,b].

Theorem: : If f(x) ∈ C¹[𝑎,b], f(𝑎) f(b) < 0 and \( \displaystyle f' (x) - p\,f(x) \ne 0 , \) where p is a real number, then the equation f(x) = 0 has a unique root in (𝑎,b).

It would be very nice if discrete models could provide calculated solutions to (ordinary and partial) differential equations exactly, but of course they do not. In fact in general they could not, even in principle, since the solution depends on an infinite amount of initial data. Instead, the best we can hope for is that the errors introduced by discretization will be small when those initial data are reasonably well-behaved.

In the second part of this chapter, we discuss some simple numerical methods applicable to first-order ordinary differential equations in normal form subject to the prescribed initial condition:

\[ y' = f(x,y), \qquad y(x_0 ) = y_0 . \]
We always assume that the slope function \( f(x,y) \) is smooth enough so that the given initial value problem has a unique solution and the slope function is differentiable so that all formulas make sense. Since it is introductory, we do not discuss discretization errors in great detail---it is a subject for another course. The user should be aware that the magnitude of global errors may depend on the behavior of local errors in ways that ordinary analysis of discretization and rounding errors cannot predict. The greatest part of this chapter is devoted to some finite difference methods that show their robustness and usefulness, and can guide a user as to how numerical methods work. Also, presented codes have limited applications because we mostly deal with uniform truncation meshes/grids that practically are not efficient. We do not discuss adaptive Runge--Kutta--Fehlberg algorithms of higher order that are used now in most applications, as well as spline methods. Instead, we present some series approximations, including the Adomian Decomposition Method (ADM), as an introduction to applications of recurrences in practical life.

The finite difference approximations have a more complicated and intricate ""physics"" than the equations they are designed to simulate. In general, recurrences usually have more properties than their continuous analogous counterparts. However, this irony is no paradox; the finite differences are used not because the numbers they generate have simple properties, but because those numbers are simple to compute.

Numerical methods for the solution of ordinary differential equations may be put in two categories---numerical integration (e.g., predictor-corrector) methods and Runge--Kutta methods. The advantages of the latter are that they are self-starting and easy to program for digital computers but neither of these reasons is very compelling when library subroutines can be written to handle systems of ordinary differential equations. Thus, the greater accuracy and the error-estimating ability of predictor-corrector methods make them desirable for systems of any complexity. However, when predictor-corrector methods are used, Runge--Kutta methods still find application in starting the computation and in changing the interval of integration.

In this chapter, we discuss only discretization methods that have the common feature that they do not attempt to approximate the actual solution y = ϕ(x) over a continuous range of the independent variable. Instead, approximate values are sought only on a set of discrete points \( x_0 , x_1 , x_2 , \ldots . \) This set is usually called the grid or mesh. For simplicity, we assume that these mesh points are equidistant: \( x_n = a+ n\,h , \quad n=0,1,2,\ldots . \) The quantity h is then called the stepsize, stepwidth, or simply the step of the method. Generally speaking, a discrete variable method for solving a differential equation consists of an algorithm which, corresponding to each lattice point xn, furnishes a number yn which is to be regarded as an approximation to the true value \( \varphi (x_n ) \) of the actual solution at the point xn.

A complete analysis of a differential equation is almost impossible without utilizing computers and corresponding graphical presentations. We are not going to pursue a complete analysis of numerical methods, which usually requires:

In this chapter, we concentrate our attention on presenting some simple numerical algorithms for finding solutions of the initial value problems for first order differential equations in normal form:

\[ y' = f(x,y), \qquad y(x_0 ) = y_0 , \]
assuming that the given problem has a unique solution.


Fundamental Inequality

We start with estimating elements of sequences.

Theorem: Suppose that elements of a sequence \( \left\{ \xi_n \right\} , \ n=0,1,2,\ldots , \) satisfy the inequalities of the form

\[ | \xi_{n+1} | \le A\, | \xi_n | + B , \]
where A and B are certain nonnegative constants independent of n. Then
\[ | \xi_{n} | \le A^n \, | \xi_0 | + B \times \begin{cases} \frac{A^n -1}{A-1} , & \quad \mbox{if } A \ne 1, \\ n , & \quad \mbox{if } A =1. \end{cases} \qquad\qquad⧫ \]

Theorem: If A is a quantity of the form 1 + δ , then \( A = 1 + \delta < e^{\delta} \) and we have

\[ | \xi_{n} | \le e^{n\delta} \, | \xi_0 | + B \, \frac{e^{n\delta} -1}{\delta} . \]
The last inequality is more readily amendable for simplification.    ⧫

We need to recall a very important definition.

A function f is said to satisfy the Lipschitz's condition (or be Lipschitz continuouss) on some interval if there exists a positive constant L such that
\[ | f(y_1 ) - f(y_2 ) | \le L\, |y_1 - y_2 | \]
for every pair of points y1 and y2 from the given interval.   ▣

Rudolf Otto Sigismund Lipschitz (1832--1903) was born in Königsberg, Germany (now Kaliningrad, Russia). He entered the University of Königsberg at the age of 15 and completed his studies at the University of Berlin, from which he was awarded a doctorate in 1853. By 1864 he was a full professor at the University of Bonn, where he remained the rest of his professional life.

Theorem (Fundamental Inequality): Let f(x,y) be a continuous function in the rectangle \( [a,b] \times [c,d] \) and satisfying the Lipschitz condition

\[ | f(x,y_1 ) - f(x,y_2 ) | \le L\, |y_1 - y_2 | \]
for some positive constant L and all pairs y1, y2 uniformly in x. Let \( y_1 (x) \quad\mbox{and}\quad y_2 (x) \) be two continuous piecewise differentiable functions satisfying the inequalities
\[ | y'_1 (x ) - f(x, y_1 (x)) | \le \epsilon_1 , \qquad | y'_2 (x) -f(x, y_2 (x))| \le \epsilon_2 \]
with some positive constants ε1,2. If in addition, these functions differ by the small amount δ > 0 at some point:
\[ | y_1 (x_0 ) - y_2 (x_0 ) | \le \delta , \]
\[ | y_1 (x ) - y_2 (x ) | \le \delta \,e^{L| x- x_0 |} + \frac{\epsilon_1 + \epsilon_2}{L} \left( e^{L|x- x_0 |} -1 \right) . \qquad\qquad⧫ \]


Horner's Rule

As an illustration of an iteration procedure, we consider a very important problem that has a long history. We will discuss the marvelous problem of polynomial evaluation at any given point x = x0 that embraces different topics---synthetic division and Taylor polynomials (that is our future topic in chapter 5). Let us consider a polynomial in the variable x with real coefficients pn(x) of degree n

\[ p_n (x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n \qquad \mbox{or} \qquad p_n (x) = a_n x^n +a_{n-1} x^{n-1} + \cdots a_1 x + a_0 . \]
This polynomial is written in standard representation (either in increasing or decreasing form---both are common in the literature) and contains n additions and also n multiplications of coefficients 𝑎k by monomials xk. Such monomials can be computed either with a naive algorithm by sequential multiplications or utilizing a standard command raising to a power that is available to you by a computer package (for instance, Mathematica). Since you don't have control over the latter and we do not discuss optimization properties of the standard command, we will stick with the naive algorithm. In this case, you have to evaluate each monomial separately by multiplying x by itself. To obtain xk, you have to multiply x by itself k-1 times, so the total number of operations for evaluating all monomials becomes
\[ 0 + 1 + 2 + \cdots + n-1 = \frac{n(n-1)}{2} . \]
Therefore, the total number of arithmetic operations (that are called flops in scientific computations) required for polynomial evaluation in standard form is
\[ 2n+ \frac{n(n-1)}{2} = \frac{n(n+3)}{2} . \]

To check our conclusion, let us write a third-degree polynomial

\[ p_3 (x) = a + b*x + c*x*x + d*x*x*x . %\quad\mbox{without memory, and } \quad q_3 (x) = a + b*x + c*x^2 + d*x^3 %\quad\mbox{with memory}. \]
To evaluate p3(x), we need 3 additions and 6 multiplications of total
\[ \left( 2n + \frac{n(n-1)}{2} \right)_{n=3} = 6 + \frac{3\cdot 2}{2} = 9 . \]

Now we present a mesmerizing mathematical method called Horner's rule (or Horner's method, Horner's scheme, etc) for polynomial evaluation that requires fewer flops than any other algorithm. This method was described by the British mathematician William George Horner in 1819. However, it was published in 1803 by the Italian mathematician Paolo Ruffini (1765--1822); it can be traced back many hundreds of years to Chinese and Persian mathematicians (see historical notes in the article by A. Pathan and T. Collyer). Horner's rule is the most efficient method of evaluating a dense polynomial at a particular value, both in terms of the number of operations and even in terms of the number of registers. Thus, in any application where such evaluations are required, it is fast and efficient, and usually overlooked. Moreover, it also allows one to determine the values of function's derivatives at the same point recursively.

We can rewrite the same polynomial pn(x) in nested form by factoring out each power of x as far as it will go to obtain

\[ p_n (x) = a_0 + x \left( a_1 + x \left( a_2 + \cdots + x \left( a_{n-1} + a_n x \right) \cdots \right) \right) . \]
It requires n additions and n multiplications, so the total number of arithmetic operations is 2n, which is less than the standard evaluation, even with access to the memory.

Instead of utilizing the nested form for polynomial evaluation, we employ Ruffini's idea and exploit the polynomial remainder theorem that states that any polynomial f(x) of degree n can be represented as

\[ f(x) = \left( x - \alpha \right) q(x) + r , \]
where r is the remainder and q(x) is called the quotient polynomial, which has degree n-1. If the remainder is zero, the quotient polynomial is referred to as the deflated polynomial. Once you set x = α, you get
\[ f(\alpha ) = r . \]
Taking the derivative, we obtain
\[ f' (x ) = q(x) + \left( x - \alpha \right) q'(x) \qquad \Longrightarrow \qquad f' \left( \alpha \right) = q \left( \alpha \right) . \]
So we see that the derivative of the function f at the same point x = α is the value of the remainder polynomial evaluated at the same point.
Algorithm Ruffini's rule for the derivative of the polynomial
  1. Set p = a[n] and q = 0
  2. Do steps 3 and 4 for i from n-1 to 0, decreasing by 1
  3. set q = p + x0 * q
  4. set p = a[i] + x0 * p
  5. The value of f(x0) is p and the value of its derivative f'(x0) is q

The nested form of a polynomial contains the basic operation involving addition and multiplication: 𝑎 + x*( ). We saw the same operation previously when one polynomial is divided into another one. As we know from precalculus, the quotient polynomial is obtained by Synthetic Division. Although Mathematica has a dedicated command: HornerForm, we implement Synthetic Division directly.
synthetic[l_List, n_, x_] := Module[{b, k, rem}, Table[b[k], {k, 2, n}];
b[n] = l[[n]];
For[k = n - 1, 1 <= k, k--, b[k] = x*b[k + 1] + l[[k]]];
rem = b[1]; Print["remainder = ", rem] ;
B = Table[b[k], {k, 2, n}]]

In the code above, l is the list of coefficients written in increasing order:

\[ l = \left\{ a_0 , a_1 , \ldots , a_n \right\} \qquad\mbox{corresponding to polynomial} \quad p_n (x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n . \]
Parameter n corresponds to the length of the array l, so if we have a polynomial of degree n, then the length of the array l is set to be n+1. The last entry, x, is the point of evaluation of the given polynomial, pn(x). The output of this routine is the remainder, abbreviated as rem, and the array B = [ q1 q2qn = 𝑎n ] that contains the list of coefficients of the quotient polynomial. Since this polynomial is one degree less than the given one, pn(x), we will denote by \( p_{n-1}(x) = a_n x^{n-1} + q_{n-1} x^{n-2} + \cdots + q_2 x + q_1 \) the quotient polynomial. Synthetic division is usually organized in the table.

Coefficients:     𝑎n     𝑎n-1     𝑎n-2         𝑎1     𝑎0
x0         x0*𝑎n     x0*qn-1 = x0(x0*𝑎n + 𝑎n-1)         x0*q2     x0*q1
Quotient:     𝑎n     qn-1 = x0*𝑎n + 𝑎n-1     qn-2 = x0(x0*𝑎n + 𝑎n-1) + 𝑎n-2         q1 = x0*q2 + 𝑎1     Remainder: x0*q1 + 𝑎0


Let us summarize what we achieved with synthetic division. With 2n flops (n additions and n multiplications), we determine the value of a polynomial with real coefficients of degree n

\[ f (x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n , \]
at an arbitrary point x = x0 along with the quotient polynomial pn-1(x):
\[ f (x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n = \left( x- x_0 \right) p_{n-1} (x) + f \left( x_0 \right) . \]
Since the quotient polynomial has degree n-1, we can apply synthetic division again and obtain
\[ p_{n-1} (x) = \left( x- x_0 \right) p_{n-2} (x) + p_{n-1} \left( x_0 \right) = \left( x- x_0 \right) p_{n-2} (x) + f' \left( x_0 \right) , \]
where f' is the derivative of the function f. With the next synthetic division we obtain the value of the second derivative at the point x = x0:
\[ p_{n-2} (x) = \left( x- x_0 \right) p_{n-3} (x) + \frac{1}{2!}\,f'' \left( x_0 \right) . \]
As a result, with a sequential application of synthetic division to each quotient polynomial, we obtain the remainders equal to the corresponding values of derivatives of the given polynomial f(x) evaluated at the point x = x0. Symbolically, it can be written as
\[ p_n \left( x; r \right) = \left( x- r \right) p_{n-1} (x;r) + p_n \left( r; r \right) \qquad \mbox{and} \qquad p'_n \left( r; r \right) = p_{n-1} (r;r) . \]
This allows us to build the Taylor polynomial:
\[ f (x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n = f \left( x_0 \right) + f' \left( x_0 \right) \left( x - x_0 \right) + \frac{1}{2!}\,f'' \left( x_0 \right) \left( x - x_0 \right)^2 + \cdots + \frac{1}{n!}\,f^{(n)} \left( x_0 \right) \left( x - x_0 \right)^n , \]
where each coefficient of the monomial (x - x0)k is the remainder obtained by synthetic division upon application to the corresponding quotient polynomial pn-k(x).
HornerForm[8*x^4 - 3*x^2 + 5*x - 7]
-7 + x (5 + x (-3 + 8 x^2))

Example: We consider an example of polynomial evaluation at a point and its derivative using Horner's rule. Let

\[ f(x) = 2\,x^5 - 3\,x^4 + 5\,x^2 -7\,x + 9 . \]
We choose x0 = 3, and expand the given function f(x) into Taylor's series centered at x0:
\[ f(x) = 2\,x^5 - 3\,x^4 + 5\,x^2 -7\,x + 9 = f(3) + f' (3) \left( x- 3 \right) + \frac{f'' (3)}{2!} \left( x- 3 \right)^2 + \frac{f''' (3)}{3!} \left( x- 3 \right)^3 + \frac{f^{(4)} (3)}{4!} \left( x- 3 \right)^4 + \frac{f^{(5)} (3)}{5!} \left( x- 3 \right)^5 . \]
So we need to determine the values of its derivatives at one given point. Using synthetic division, we obtain
Coefficients:     2     -3     0     5     -7     9
x0 = 3         3*2 = 6     3*3 = 9     3*9 = 27     32*3 = 96     89*3 = 267
Quotient:     2     3*2 - 3 = 3     9+0 = 9     3*9 + 5 = 32     3*32 - 7 = 89     Remainder = 276

We check with Mathematica
synthetic[{9, -7, 5, 0, -3, 2}, 6, 3]
remainder = 276
{89, 32, 9, 3, 2}
Next application of synthetic division gives
Coefficients:     2     3     9     32     89
x0 = 3         3*2 = 6     3*9 = 27     3*36 = 108     3*140 = 420
Quotient:     2     3*2 + 3 = 9     27 + 9 = 36     3*36 + 32 = 140     Remainder = 509

synthetic[{89, 32, 9, 3, 2}, 5, 3]
remainder = 509
{140, 36, 9, 2}
So we know that
\[ f(3) = 276 \qquad\mbox{and} \qquad f' (3) = 509. \]
We proceed:
Coefficients:     2     9     36     140
x0 = 3         3*2 = 6     3*15 = 45     3*81 = 243
Quotient:     2     3*2 + 9 = 15     45 + 36 = 81     Remainder = 383

synthetic[{140, 36, 9, 2}, 4, 3]
remainder = 383
{81, 15, 2}
So     \( \displaystyle \frac{f'' (3)}{2} = 383 .\)
Coefficients:     2     15     81
x0 = 3         3*2 = 6     3*21 = 63
Quotient:     2     3*2 + 15 = 21     Remainder = 144

synthetic[{81, 15, 2}, 3, 3]
remainder = 144
{21, 2}
The last step is easy:
Coefficients:     2     21
x0 = 3         3*2 = 6
Quotient:     2     Remainder = 27

synthetic[{21, 2}, 2, 3]
remainder = 27
Therefore, we found all remainders and all derivatives evaluated at x = 3:
\[ f(x) = 2\,x^5 - 3\,x^4 + 5\,x^2 -7\,x + 9 = 2 \left( x-3 \right)^5 + 27 \left( x-3 \right)^4 + 144 \left( x-3 \right)^3 + 386 \left( x-3 \right)^2 + 509 \left( x-3 \right) + 276 . \]
Horner's rule for polynomial evaluation
% the polynomial coefficients are stored in the array a(j), j=0,1,..,n
% px = value of polynomial upon completion of the code

px = a(n)
for k = n-1 downto 0
px = a(k) + px*x

A similar algorithm for evaluating the derivative of a polynomial:

Horner's rule for polynomial derivative evaluation
% the polynomial coefficients are stored in the array a(j), j=0,1,..,n
% dp = value of derivative upon completion of the code

dp = n*a(n)
for k = n-1 downto 1
dp = k*a(k) + dp*x
Horner's rule was invented by William George Horner (1786--1837) who was born in Bristol, England, and spent much of his life as a schoolmaster in Bristol and, after 1809, in Bath. His paper on solving algebraic equations, which contains the algorithm we know as Horner's rule, was published in 1819, although the Italian mathematician Paolo Ruffini (1765--1822) had previously published a very similar idea.

There is a more efficient algorithm for polynomial evaluation:

A more efficient implementation of Horner's rule for polynomial evaluation
% the polynomial coefficients are stored in the array a(j), j=0,1,..,n
% b(n) = value of polynomial upon completion of the code

b(0) = a(n)
for k = 1 to n
b(k) = x*b(k-1) + a(n-k)

% c(n-1) = value of derivative upon completion of the code

c(0) = b(0)
for k = 1 to n-1
c(k) = x*c(k-1) + b(k)



  1. Bradie, B.A., Friendly Introduction to Numerical Analysis, Pearson Education,2007
  2. Burden, R.L. and Faires, J.D., Numerical Analysis, Cengage Learning, 10 edition, 2015.
  3. Chapra, S.C. and Canale, R.P., Numerical Methods, Online, The Hong Kong University of Science and Technology. 2012.
  4. Chasnov, J.R,, Numerical Methods, McGraw Hill,2009.
  5. Gerald C.F. and Wheatley, P.O., Applied Numerical Analysis, Pearson Education, 7th edition, 2003.
  6. Hauser, John R., Numerical Methods for Nonlinear Engineering Models, ISBN 978-1-4020-9920-5, 2009, Springer Netherlands, doi: 10.1007/978-1-4020-9920-5
  7. Ostrowski, A.M., Solutions of Equations and Systems of Equations, Academic Press, New York, 1960.
  8. Pathan A. and Collyer, T., The Wonder of Horner's Method, The Mathematical Gazette, Vol. 87, No. 509 (Jul., 2003), pp. 230-242
  9. Petković, M.S., Neta, B., Petković, L.D., Džunić, J., Multipoint Methods for Solving Nonlinear Equations, Elsevier, Amsterdam, 2013.
  10. Traub, J.F., Iterative Methods for Solution of Equations, Prentice-Hall, Englewood Cliffs, NJ, 1964.
  11. Wu, X., Wu, H.W., On a class of quadratic convergence iteration formulae without derivatives, Applied Mathematics & Computations, 2000, Vol. 10, Issue 7, pp. 381--389.


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 Problems)