In most areas of numerical analysis the first step in the solving a particular problem is to find a technique that can be used to obtain the solution to the problem. The numerical solution of ordinary differential equations is somewhat different. More often the first step is to choose that technique among the many available that will serve the purpose best. There seems to be no method that is the best in all situations.


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.

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.

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. matlab has built-in functions fzero, roots (for polynomial roots), and chebfun 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.

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

  • If every step value xk+1 in an iteration method \( \displaystyle x_{k+1} = \phi \left( x_k \right) \) depends only on the one previous step xk, then the method is called a one-point iteration. Newton's method is an example of such a method. In differential equations, such an iteration method is called a one-step method.
  • If the iteration function \( \displaystyle \phi \left( x_k , x_{k-1} , \ldots , x_{k-n+1} \right) \) depends on n previous steps, the iteration method is called a one-point method with memory. An example of an iteration function with memory is the well-known secant method. In differential equations, such a method is called a multistep method.
  • Another type of iteration functions is derived by using the expressions w1(xk), w2(xk), … , wm(xk), where xk is the common argument. The iteration function ϕ, defined as
    \[ x_{k+1} = \phi \left( x_k , w_1 (x_k ), w_2 (x_k ) , \ldots , w_m (x_k ) \right) , \qquad k=0,1,2,\ldots , \]
    is called a multipoint iteration function without memory. The simplest examples are Steffensen's method and Traub--Steffensen's method. In differential equations, such a method is usually referred to as the Runge--Kutta method.
  • Assume that an iteration function ϕ has arguments zj, where each argument represents m + 1 quantities xj, w1(xj), w2(xj), … , wm(xj), m ≥ 1. The iteration function
    \[ x_{k+1} = \phi \left( z_k , w_1 (z_k ), w_2 (z_k ) , \ldots , w_m (z_k ) \right) , \qquad k=0,1,2,\ldots , \]
    is called a multipoint iteration function with memory.
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:

  • an understanding of computational tools;
  • an understanding of the problem to be solved;
  • construction of an algorithm that will solve the given mathematical problem.

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 , \]
then
\[ | 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, matlab). Since you don't have control on 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 matlab has a dedicated command: deconve for division of two polynomials and polynomialReduce, we implement Synthetic Division directly.

For example, suppose we want to divide the polynomial p(x) = 2x³ + 3x² - x + 4 by u(x) = x - 5. We enter coefficients of the polynomial p(x) in decreasing order and execute the following matlab code

In the code above, q is the quotient polynomial whose coefficients are written in decreasing order as an array of dimension 1 less than the polynomail p(x), and r is the the array that contains the remainder in its last entry. The remainder can be obtained with the aid of polynomialReduce. In our case, just type in 

polynomialReduce(2*x^3 +3*x^2 -x+4,x-5)

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

\[ \text{dividend} = \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 q2 … qn = 𝑎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

Conclusions

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

Actually, matlab has a dedicated command for Horner's rule: horner. For example,

matlab also has a built-in function polyval that uses Horner's rule to evaluate polynomials. The list of coefficients is in descending order of power, where as to task spec specifies ascending order. For example, the polynomial \( p(x) = -7*x^4 + 5*x^3 - 3*x^2 + 8 \) can be evaluated at x = 2 as 

polyval(fliplr([8, 0, -3, 5,-7]),2)


ans = -76

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 matlab 

polynomialReduce(2*x^5 -3*x^4 + 5*x^2 -7*x+9,x-3)
which gives the remainder to be 276. The coefficients of the corresponding quotient polynomial are
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 


[q1,r]=deconv(q,[1 -3])
with output: 

q1 = [ 2     9    36   140 ]
r = 0     0     0     0   509
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 


[q2,r]=deconv(q1,[1 -3])
with output: 

q2 = [2    15    81]
r = 0     0     0   383
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 


[q3,r]=deconv(q2,[1 -3])
with output: 

q2 = [ 2    21 ]
r = 0     0   144
The last step is easy:
Coefficients:     2     21
x0 = 3     ↡     3*2 = 6
Quotient:     2     Remainder = 27 


[q4,r]=deconv(q3,[1 -3])
with output: 

q4 = 2
r = 0    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
endfor

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

% 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)
endfor

 

 

  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.

=================== to be checked ==================

So far we have found an efficient procedure for evaluating a polynomial P(x) at any x = x0. Next we develop an algorithm for getting its derivative P'(x0) at the same point x0. Note that we can divide P(x) by x - x0 for any point x0 ro get quotient and remainder:
\[ P(x) = Q(x) \left( x- x_0 \right) + P\left( x_0 \right) , \]
where the quotient polynomial has degree n-1:
\[ P(x) = \sum_{i=0}^n a_i x^i , \qquad Q(x) = \sum_{j=0}^{n-1} b_j x^j . \]
Once we find coefficients bj of the quotient polynomial Q(x), we can use Horner's algorithm to evaluate it at any point x0. Upon substituting Maclaurin formulas for P(x) and Q(x) and equating coefficients of like power of x, we obtain
\begin{align*} b_{i-1} &= i\,a_i , \qquad i=1,2,\ldots , n. \end{align*}
However, application of Horner's rule yields \begin{align*} c_{n-1} &= a_n , \\ c_{i-1} &= a_i + x_0 c_i , \qquad i = 1,2,\ldots , n-1; \\ P\left( x_0 \right) &= a_0 + x_0 c_0 . \end{align*}

The pseudocode actually holds the coefficients ci of the quotient polynomial Q(x):

Algorithm Horner'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 P(x0) is p and the value of its derivative P'(x0) is q

Horner's rule with derivatives:

Example: Let us evaluate the polynomial of degree 6


   ■