Numerical Analysis: How to Calculate Special Functions

To contact us Click HERE
Your calculator has buttons for all sorts of special functions, like square root, sine, cosine, and logarithms. How does the calculator know how to calculate these functions? In arithmetic, you learned to add, subtract, multiply and divide "by hand". With only those four basic operations available, where would you even start if you had to calculate a square root or transcendental function? Today we are going to explore a little bit of Numerical Analysis, which is the branch of math that focuses on actually calculating things numerically, with high accuracy and efficiency.

To be concrete, let's focus on the base two logarithm function, whichwe will abbreviate L. You can review Logarithmsif you are a little rusty on the concept, but do not worry:you won't have to calculate with them yourself. The whole point is tolook behind the curtain and see how the computer does it. All we care about is that L is a curved function, not astraight line, and that there is no simple formula for calculating it,which makes it challenging.

Let's start with a graph of the L function. We will create it using R,which (as regular readers will recall) is a free, high quality, opensource programming environment well suited to statistics and relatedcalculations. Read Koch Snowflakes for detailed instructions on downloading andinstalling R. Now run the following code:

L <- log2x <- 1:400/10plot(x, L(x), type="l", cex.axis=1.5,      cex.lab=1.5, xlab="x", ylab="L(x)")abline(0,0)

This produces the following graph.

Logarithms in any base are zero at x=1. They go down to negativeinfinity as x approaches zero, as shown by the steep part at the left of thegraph. They grow without limit as x grows, but very slowly, asindicated by the relatively flat part toward the right side of thegraph. Recall that L(2^n) is n, so for instance, L(16) is 4, L(32) is5, and L(64) is 6. Even L(1024) is just 10, so the curve grows everflatter as you move further right.

Now to the question of computing L. Recall that the basic propertyof logarithms is that they translate multiplication into addition. Insymbols:

L(a*b) = L(a) + L(b)

We can use this property of L to our advantage. All we really need to figureout is how to compute L(x) when x is between 1 and 2. That's becauseif x is outside that range, we can double or halve it repeatedly untilwe get into that range; every time we do so, we adjust the answer upor down by one.

More precisely, suppose x =f*(2^n) where n is an integer and f is a number between oneand two. Then L(x) = L(f) + n. So all we really need to figureout is how to compute L(f), where f is between 1 and 2. Adding 'n' iseasy because computers represent real numbers internally usingfloating point, essentially storing both n and f directly.

If we redraw the graph of L over just the range from 1 to 2, itlooks much straighter. For comparison, the red line shows theapproximation y = x-0.95, which seems like a decent fit.

So as a first approximation, we could say that L(x) is roughly x-0.95, atleast on the interval from x=1 to x=2 that is our focus. However, thatis not nearly good enough: computer calculations with real numberstypically are accurate to 16 decimal digits, so we should aim forthat as well.

However, this gives us an idea. If a straight line was a pretty goodfit, might we be able to find a quadratic that is an even better fit?Indeed we can. In fact, calculus provides a useful tool here, calledTaylor series. A Taylor expansion is just a polynomial that is a"good" fit to the function at a particular point. We can ask for alinear fit, a quadratic fit, or even higher order polynomials. Thehigher the order, the better the match, since the Taylor polynomial isdesigned to match the function value, its slope, and as many higherorder derivatives (e.g. curvature) as it can based on its order.

Recall the formula for Taylor series from Calculus:

To use it, we need to know the value of the function and itsderivatives at a specific point x0. To avoid circularity, we assumethat we do not know how to compute logarithms yet, so we need tochoose x0 to be someplace where we know the answers inadvance. Choosing x0=1 works well, since log(1) is just zero.

Since we have to calculate derivatives, it is handy to switch to thenatural log function log(x) instead of the base two log L. (Some people write 'ln' for natural log, instead of 'log'.)Since

log2(x) = log(x)/log(2)

we can easily get back to log base 2 (or any other base), once we knowhow to compute natural logs on values of x between 1 and 2.

Taylor polynomials tend to be most accurate close to the expansionpoint x0, and less accurate further away. Since we want to expandaround x0=1, we will get better results if we shift the range of xvalues to put x0 near the center rather than at an end. We know that

log(x) = log(x/sqrt(2)) + log(2)/2

so it is enough to be able to calculate log(x) for values of x between1/sqrt(2) and sqrt(2), or 0.7 to 1.42.

Putting all these pieces together, the Taylor series for natural logof x around x0=1 is

The first few terms are a good approximation near x=1.The next chart shows log(x) in black, with thequadratic x-1-(x-1)^2/2 in red.

The red line seems to be a very good fit except near the left andright edges of the chart. It is goodenough that we need to switch to plotting the error, rather thanthe curves themselves, in order to really see what is happening.

The next plot therefore shows the error, A(x)-log(x), where A is ourapproximation, here the Taylor quadratic centered at x=1.

x <- 0.7+0.72*(0:400)/400A <- function(x) { x-1-(x-1)^2/2 }plot(x, A(x)-log(x), type="l", xlab="x",      ylab="A(x)-log(x)", lwd=3)abline(0,0)

Now we can see that the error is smaller near x0=1 than faraway. Even at the edges of the graph it is never more than 0.02 insize. Alas, this is still much too big to be practical.

However, we can go further and try a cubic Taylor polynomial. The nextchart adds a red curve to the previous chart. The red curve is theerror when we add the (x-1)^3/3 term to our approximation.

Clearly, the error is smaller, but not enormously so. This suggeststhat to get errors as small as 10^(-16) will require using a lot ofterms in the Taylor series.

An easy way to compute Taylor polynomials is to use Maxima,a free, high quality, open source symbolic mathematics programdescribed in detailin Calculus Calculator.The following code computes the 4th order polynomial around x0=1.

taylor(log(x), x, 1, 4);string(%);

The first line computes the Taylor series. The second line prints the polynomial in "input" notation suitable for graphing in R.The result is

x-1-(x-1)^2/2+(x-1)^3/3-(x-1)^4/4

The problem is apparent: the coefficient to the (x-1)^n term is only1/n, rather than 1/(n!). The symbol n! means n-factorial, i.e. theproduct of the numbers from 1 to n. So for example, 5! is 1*2*3*4*5,or 120. Factorials get big fast: 10! is 3,628,800. If the expansionfor log involved coefficients that went to zero like 1/(n!), it wouldconverge much more quickly. As it stands, 1/n gets small slowly, whichfits with our earlier observation that going from a quadratic to acubic did not help that much.

What if we use the first 10 terms, instead of a quadratic or cubic?We can reduce the error even further, as shown in the next chart.

Now the vertical scale is 4e-06, which means 4 times 10^(-6), or0.000004. That is a big improvement over the 0.02 we had before, butit suggests that we will still need a lot more terms before wefinally get to 10^(-16) accuracy.

Very high order polynomials with large numbers of terms can be numericallyunstable: small amounts of error in the coefficients can get magnifiedinto big errors in the result. Soadding a lot more terms might not actually meet our goal, though it would beworth a try. Also, using lots of terms takes a long time to calculate, so if we couldfind a more efficient process, that would be nice.

Is there something else we can try?

The inverse of the natural log function is the exponential function'exp'. The Taylor series for exp converges much more rapidly, becausethe coefficients scale like 1/(n!) instead of like 1/n. That suggests wecan use Taylor series to efficiently calculate exp. Once we have donethat, can we use it to somehow find natural logs?

Indeed we can. Newton's Method is a powerful algorithm forfinding roots (zeros) of functions. It uses the linear Taylorpolynomial: since f(x) is approximately f(x0)+f'(x0)(x-x0), then iff(x) is zero, we can solve for the root at x = x0 - f(x0)/f'(x0).

If the function f is not exactly linear, this newvalue of x will not be exactly the root, but if x0 was a good initialguess, the new x will be an even better approximation to the true root. Thismeans that if we iterate a few times, using

we should quickly reach a very accurate answer.

Suppose we want to calculate log(c) for some number c.We choose f(x) = exp(x) - c. This function will be zero when x =log(c). Applying Newton's method means using the iteration

xnew = xold - 1 + c*exp(-xold)

Try the following R code:

A <- function(c) {  x <- c - 1 # initial guess  for(i in 1:5) x <- x - 1 + c*exp(-x)  x       }plot(x, A(x)-log(x), type="l", xlab="x",     ylab="A(x)-log(x)", lwd=3)

Starting from c-1 as our initial guess and iterating 5 times yields an answeraccurate to 10^(-16), just as we had wanted, as shown in the following chart:

The jagged, noisy pattern means we have reached the limit of accuracyon this computer: all that is left is random-looking round-off error due tothe finite nature of computer arithmetic.

To achieve this result, we did need to call the exponential functionfive times. That could be expensive, so we should examine how manyterms we need in the Taylor series for 'exp'. Let's check the errorwhen we use a 10th degree polynomial.

The Taylor expansion for 'exp' is simplest if we expand around x0=0,so we will once again shift the range of interest.Since exp(-x) is just 1/exp(x), and since exp(x-1) is exp(x)/e,we can now focus on the range from x = -0.3 to x = 0.42, assuming we knowthe value e = 2.718282... accurately already.

We are pleased to see errors in the 10^(-12) range, a greatimprovement over the 10^(-6) we had with the log function at tenterms. Why the difference? Look at the Taylor series for exp:

Notice the coefficients. The 24 is 4! and the 120 is 5!, so we knowthat by the tenth term, the coefficients are getting verysmall. Therefore, dropping the eleventh and higher order terms doesnot introduce much error. As a result, we can get a highly accurateapproximation to 'exp' just by using a few more terms in this series.

We have used graphs of numerical experiments to assess the accuracy ofour approximations. This is of course somewhat circular: if mycomputer did not already have a working version of the log function, Iwould not have been able to make these graphs. The mathematicalside of Numerical Analysis, therefore, consists of proving inequalitybounds on approximations. For example, using Taylor's theorem in theform that includes a "remainder" term allows us to put an upper boundon how large the error can possibly be, even if we do not yet know howto calculate log or exp. That breaks the circularity and gives us (orcomputer designers, anyhow) peace of mind.

If you liked this article, you may also like System Dynamics: Feedback Models, or check out the Contents page for a complete list of past topics.

Please post questions, comments and other suggestions using the box below, or email me directly at the address given by the clues at the end of the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can follow @ingThruMath on Twitter, where I will tweet about each new post to this blog. The Contents page has a complete list of previous articles in historical order. Now that there are starting to be a lot of articles, you may also want to use the 'Topic', 'Search' or 'Archive' widgets in the side-bar to find other articles of related interest.

I hope you enjoyed this discussion. You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on. Seeyou next time!