Algorithms - Fundamentals of Numerical Computation

An idealized mathematical problem $f(x)$ can usually only be approximated using a finite number of steps in finite precision. A complete set of instructions for transforming data into a result is called an algorithm. In most cases it is reasonable to represent an algorithm by another mathematical function, denoted here by $\tilde{f}(x)$ .

Even simple problems can be associated with multiple algorithms.

Example 1.3.1

Suppose we want to find an algorithm that maps a given $x$ to the value of the polynomial $f(x)= 5x^3 + 4x^2 + 3x + 2$ . Representing $x^2$ as $(x)(x)$ , we can find it with one multiplication. We can then find $x^3=(x)(x^2)$ with one more multiplication. We can then apply all the coefficients (three more multiplications) and add all the terms (three additions), for a total of 8 arithmetic operations.

There is a more efficient algorithm, however: organize the polynomial according to Horner’s algorithm,

f(x) = 2 + x \bigl( 3 + x( 4 + 5x) \bigr).

(1.3.1)

In this form you can see that evaluation takes only 3 additions and 3 multiplications. The savings represent 25% of the original computational effort, which could be significant if repeated billions of times.

1.3.1Algorithms as code¶

Descriptions of algorithms may be presented as a mixture of mathematics, words, and computer-style instructions called pseudocode, which varies in syntax and level of formality. In this book we use pseudocode to explain the outline of an algorithm, but the specifics are usually presented as working code.

Of all the desirable traits of code, we emphasize clarity the most after correctness. We do not represent our programs as always being the shortest, fastest, or most elegant. Our primary goal is to illustrate and complement the mathematical underpinnings, while occasionally pointing out key implementation details.

As our first example, Function 1.3.1 implements an algorithm that applies Horner’s algorithm to a general polynomial, using the identity

\begin{split} p(x) &= c_1 + c_2 x + \cdots + c_n x^{n-1} \\ &= \Bigl( \bigl( (c_n x + c_{n-1}) x + c_{n-2} \bigr) x + \cdots +c_{2} \Bigr)x + c_{1}. \end{split}

(1.3.2)

Algorithm 1.3.1 (horner)

Julia

MATLAB

Python

Horner’s algorithm for evaluating a polynomial

1
2
3
4
5
6
7
8
9
10
11
12
13
14
"""
    horner(c,x)

Evaluate a polynomial whose coefficients are given in ascending
order in `c`, at the point `x`, using Horner's rule.
"""
function horner(c,x)
    n = length(c)
    y = c[n]
    for k in n-1:-1:1
        y = x*y + c[k]
    end
    return y
end

Horner’s algorithm for evaluating a polynomial

1
2
3
4
5
6
7
8
9
10
11
12
13
14
"""
    horner(c,x)

Evaluate a polynomial whose coefficients are given in ascending
order in `c`, at the point `x`, using Horner's rule.
"""
function horner(c,x)
    n = length(c)
    y = c[n]
    for k in n-1:-1:1
        y = x*y + c[k]
    end
    return y
end

Horner’s algorithm for evaluating a polynomial

horner.py

1
2
3
4
5
6
7
8
9
10
11
12
def horner(c,x):
    """
    horner(c,x)

    Evaluate a polynomial whose coefficients are given in descending order 
    in `c`, at the point `x`, using Horner's rule.
    """
    n = len(c)
    y = c[0]
    for k in range(1, n):
        y = x * y + c[k]   
    return y

Example 1.3.2

Julia

MATLAB

Python

Using a function

Here we show how to use Function 1.3.1 to evaluate a polynomial. It’s not a part of core Julia, so you need to download and install this text’s package once, and load it for each new Julia session. The download is done by the following lines.

import Pkg
Pkg.add("FundamentalsNumericalComputation");

Once installed, any package can be loaded with the using command, as follows.

Many Julia functions, including the ones in this text, are in packages that must be loaded via using or import in each session. Sometimes a using statement can take a few seconds or even minutes to execute, if packages have been installed or updated.

using FundamentalsNumericalComputation

For convenience, this package also imports many other packages used throughout the book and makes them available as though you had run a using command for each of them.

:::{card}

If you are not sure where a particular function is defined, you can run `methods` on the function name to find all its definitions.

Returning to horner, let us define a vector of the coefficients of $p(x)=(x-1)^3=x^3-3x^2+3x-1$ , in ascending degree order.

c = [-1, 3, -3, 1]

In order to avoid clashes between similarly named functions, Julia has boxed all the book functions into a namespace called FNC. We use this namespace whenever we invoke one of the functions.

:::{card}

You must use the module name when a package is loaded by `import`, but when loaded via `using`, some functions may be available with no prefix.

FNC.horner(c, 1.6)

The above is the value of $p(1.6)$ .

While the namespace does lead to a little extra typing, a nice side effect of using this paradigm is that if you type FNC. (including the period) and hit the Tab key, you will see a list of all the functions known in that namespace.

The multi-line string at the start of Function 1.3.1 is documentation, which we can access using ?FNC.horner.

Using a function

Here we show how to use horner to evaluate a polynomial. First, we have to ensure that the FNC package is imported.

import FNC

Here is the help string for the function:

help(FNC.horner)

We now define a vector of the coefficients of $p(x)=(x−1)^3=x^3−3x^2+3x−1$ , in descending degree order. Note that the textbook’s functions are all in a namespace called FNC, to help distinguish them from other Python commands and modules.

c = array([1, -3, 3, -1])
print(FNC.horner(c, 1.6))

The above is the value of $p(1.6)$ , up to a rounding error.

The quoted lines at the beginning of Function 1.3.1 are a documentation string. The function itself starts off with the keyword function, followed by a list of its input arguments. The first of these is presumed to be a vector, whose length can be obtained and whose individual components are accessed through square bracket notation. After the computation is finished, the return keyword indicates which value or values are to be returned to the caller.

The Polynomials package for Julia provides its own fast methods for polynomial evaluation that supersede our simple Function 1.3.1 function. This will be the case for all the codes in this book because the problems we study are well-known and important. In a more practical setting, you would take implementations of basic methods for granted and build on top of them.

1.3.2Writing your own functions¶

Any collection of statements organized around solving a type of problem should probably be wrapped in a function. One clue is that if you find yourself copying and pasting code, perhaps with small changes in each instance, you should probably be writing a function instead.

Julia

MATLAB

Python

Functions can be defined in text files with the extension .jl, at the command line (called the REPL prompt), or in notebooks.

As seen in Function 1.3.1, one way to start a function definition is with the function keyword, followed by the function name and the input arguments in parentheses. For example, to represent the mathematical function $e^{\sin x}$ , we could use

function myfun(x)
    s = sin(x)
    return exp(s)
end

The return statement is used to end execution of the function and return one or more (comma-separated) values to the caller of the function. If an executing function reaches its end statement without encountering a return statement, then it returns the result of the most recent statement, but this is considered poor style.

For a function with a short definition like the one above, there is a more compact syntax to do the same thing:

myfun(x) = exp(sin(x))

You can also define anonymous functions or lambda functions, which are typically simple functions that are provided as inputs to other functions. This is done with an arrow notation. For example, to plot the function above (in the Plots package) without permanently creating it, you could enter

plot(x -> exp(sin(x)), 0, 6)

As in most languages, input arguments and variables defined within a function have scope limited to the function itself. However, they can access values defined within an enclosing scope. For instance:

mycfun(x) = exp(c*sin(x))
c = 1;  mycfun(3)   # returns exp(1*sin(3))
c = 2;  mycfun(3)   # returns exp(2*sin(3))

There’s a lot more to be said about functions in Julia, but this is enough to get started.

Functions can be defined in text files with the extension .m, at the command line, or in Live Scripts.

As seen in Function 1.3.1, one way to start a function definition is with the function keyword, followed one or more output arguments in brackets, the function name, and the input arguments in parentheses. For example, to represent the mathematical function $e^{\sin x}$ , we could use the following in a file called myfun.m:

function [y] = myfun(x)
    s = sin(x);
    y = exp(s);
end

Whatever value is assigned to y when the function terminates will be returned as the output of the function.

For a function with a short definition like the one above, there is a more compact syntax to do the same thing:

myfun = @(x) exp(sin(x));

The syntax on the right of the = above defines an anonymous function (called a lambda function in computer science), which can be used in place without giving it a name as we did here. We’ll have examples of doing this later on.

c = 1;
mycfun = @(x) exp(c*sin(x));
mycfun(3)   % returns exp(1*sin(3))
c = 2;  
mycfun(3)   % also returns exp(1*sin(3))
mycfun = @(x) exp(c*sin(x));   % redefines mycfun
mycfun(3)   % now returns exp(2*sin(3))

There’s a lot more to be said about functions, but this is enough to get started.

Functions can be defined in text files with the extension .py, at the command line (called the REPL prompt), or in notebooks.

As seen in Function 1.3.1, one way to start a function definition is with the def keyword, followed by the function name and the input arguments in parentheses, ending with a colon. The statements for the body of the function must then all be indented. For example, to represent the mathematical function $e^{\sin x}$ , we could use

def myfun(x):
    s = np.sin(x)
    return np.exp(s)

The return statement is used to end execution of the function and return one or more (comma-separated) values to the caller of the function.

For a function with a short definition like the one above, there is a more compact syntax to do the same thing:

myfun = lambda x : np.exp(np.sin(x))

mycfun = lambda x : np.exp(c * np.sin(x))
c = 1;  print(mycfun(3))   # exp(1*sin(3))
c = 2;  print(mycfun(3))   # exp(2*sin(3))

There’s a lot more to be said about functions in Python, but this is enough to get started.

1.3.3Exercises¶

⌨ Write a function poly1(p) that returns the value of a polynomial $p(x) = c_1 + c_2 x + \cdots + c_n x^{n-1}$ at $x=-1$ . You should do this directly, not by a call to or imitation of Function 1.3.1. Test your function on $r(x)=3x^3-x+1$ and $s(x)=2x^2-x$ .
⌨ In statistics, one defines the variance of sample values $x_1,\ldots,x_n$ by
$s^2 = \frac{1}{n-1} \sum_{i=1}^n (x_i - \overline{x})^2, \qquad \overline{x} = \frac{1}{n} \sum_{i=1}^n x_i.$
(1.3.3)
Write a function samplevar(x) that takes as input a vector x of any length and returns $s^2$ as calculated by the formula. You should test your function on the vectors ones(100) and rand(200). If you enter using Statistics in Julia, then you can compare to the results of the var function.
⌨ Let x and y be vectors whose entries give the coordinates of the $n$ vertices of a polygon, given in counterclockwise order. Write a function polygonarea(x,y) that computes and returns the area of the polygon using this formula based on Green’s theorem:
$A = \frac{1}{2} \left| \sum_{k=1}^n x_k y_{k+1} - x_{k+1}y_k \right|.$
(1.3.4)
Here $n$ is the number of polygon vertices, and it’s understood that $x_{n+1}=x_1$ and $y_{n+1}=y_1$ . (Note: The function abs computes absolute value.) Test your functions on a square and an equilateral triangle.

Preface

Problems and conditioning

Preface

Stability