<< Computer Algorithms | Use Ctrl + or – to enlarge or reduce text size. |

Problem: You are given an invertable function **F** of a single variable over a finite interval. This function might take a long time to compute. Write a program that computes a good approximation to the function’s inverse at any point, and that is always fast. You are allowed one arbitrarily long setup computation.

In other words, write a program that takes a given **y** in the function’s range and quickly computes **x** so that **F(x) = y** to a good approximation. And this even though **F** might take a long time to compute.

The solution is to use function interpolation on the inverse, using some algorithm to compute the inversion table. Any of the usual root finding algorithms (bisection, secant, Brent’s, or – if the derivative is known – Newton’s) could be used for the latter. Since this computation needn’t be very fast we use the bisection method, which is always reliable. (If you use one of the other root finding algorithms, note that **F** is strictly monotonic because it is invertable, which might simplify the computation.)

Using the bisection algorithm (or any other) by itself will not solve the original problem because it takes even longer to compute than the function. The method shown below, after its long initial computation, computes the inverse at once for any given **y**.

Before we get to the method first a description of the two algorithms it depends on.**Bisection algorithm for solving F(x) = 0**

The accuracy is determined by the constant **tol** (tolerance). No upper bound has been placed on the number of iterations. In order to avoid an infinite loop **tol** must not be smaller than the machine precision used for **x**.

Two initial numbers must be given that bracket the solution**a** ≤ **x** ≤ **b**

We assume the function is invertable over the given domain which means there is only one root.

The program below and the ones given later are written for PBCC and can be modified for Visual Basic, FORTRAN and any other BASIC-like compiler.

Function F(x As Single) As Single Function = x*x - 2 'simple example, the solution in the domain 1 < x < 2 is Sqr(2) End Function '---------------------------- Function PBMain Local i As Long Local x, x1, x2 As Single Local tol As Single x1 = 1 'bracket solution x2 = 2 tol = .000001 'fixed i = 0 'count loops, i not used otherwise Do Incr i x = (x1 + x2) * .5 If F(x) = 0 Or Abs(x1 - x2) < tol Then Exit Do If Sgn(F(x)) = Sgn(F(x2)) Then x2 = x Else x1 = x End If Loop Print Print tol; "tolerance" Print i; "iterations" Print " Solution: "; x Print " Actual solution: "; Sqr(2) WaitKey$ End Function

Given evaluations of **G(y)** at evenly spaced **y** over an interval, find a good approximation of **G(y)** for any **y** in the interval by linear interpolation.

This is useful when **G** takes a lot of time to compute and great accuracy is unnecessary. You compute **G** at a finite number of values once and for all, then use this table to approximate **G** at any value any number of times.

The greater the finite number of these values (**n**) or the shorter the interval over which they are spread (**interval**), the more accurate the interpolation.

In order to minimize truncation error we put up with a trivial time penalty by using

... * n / interval

instead of pre-computing the fraction.

Usually **y** is considered a function of **x**, but here the relation is the reverse, because we have in mind an application to the inverse of a function.

In this example the function is **G(y) = Sqr(y)**. You can change this.

The interpolation will be exact at the **y()** and less accurate in between. For example, with the fixed parameters in the program as is, where

n = 50, y0 = 0, yn = 5, G(y) = Sqr(y)

consider y = 4:

Ginterpolate(4) = 2

G(4) = 2

With 4 it’s exact, but at y = 3.95

Ginterpolate(3.95) = 1.987340

G(3.95) = 1.987461

Global y0, yn As Single 'first and last points Global interval As Single 'length of interval from y0 to yn Global n As Long 'number of divisions plus one Global x(), y() As Single 'x = G(y), y() equally spaced from y0 to yn '---------------------------- ' The larger n is (see Sub SetupG) the more accurate the interpolation. Function G(y As Single) As Single 'defined from y0 to yn Function = Sqr(y) 'simple example End Function '---------------------------- Sub SetupG n = 50 'fixed Dim x(n), y(n) As Single Local i As Long ' chose the first and last points so the line through them on the graph ' does not intersect a zero of the function y0 = 0 'first point, fixed yn = 10 'last point, fixed interval = yn - y0 For i = 0 To n - 1 y(i) = y0 + interval * i / n ' y0 + s*i where s = (yn - y0) / n x(i) = G(y(i)) Next y(n) = yn 'possibly more accurate than putting in loop x(n) = G(yn) End Sub '---------------------------- Function Ginterpolate(y As Single) As Single Local k As Long k = Int((y - y0)*n / interval) ' Int((y - y0) / s) where s = (yn - y0) / n Function = x(k) + (y - y(k)) * ( x(k+1) - x(k) ) * n / interval ' ...*n/interval = .../s End Function '---------------------------- Function PBMain Local y As Single Local a As String SetupG Do Line Input "Enter a number: "; a : If Trim$(a) = "" Then End y = Val(a) If y < y0 Or y > yn Then 'check range Print "Number must be between" + Str$(y0) + " and" + Str$(yn) Iterate End If Print "Interpolated: "; Ginterpolate(y); " From function: "; G(y) Print Loop End Function

For the purpose of demonstrating the idea described in the introduction the particular function to invert below is fast, simply **x** squared, and has a known simple solution. You can change the function in **Function F**, then either change or delete the “Actual inverse” in the output.

Global y0, yn As Single 'first and last points Global interval As Single 'length of interval from y0 to yn Global n As Long 'number of interpolation points plus one Global x(), y() As Single 'x = G(y), y() equally spaced from y0 to yn Global tol As Single '---------------------------- Function F(x As Single) As Single Function = x*x 'simple example for the purpose of demonstration End Function '---------------------------- Function G(y As Single) As Single 'the inverse of F Local x, x1, x2 As Single Local y1, y2 As Single x1 = y0 x2 = yn Do x = (x1 + x2) * .5 If F(x) = y Or Abs(x1 - x2) < tol Then Exit Do If Sgn(F(x) - y) = Sgn(F(x2) - y) Then x2 = x Else x1 = x End If Loop Function = x End Function '---------------------------- Sub SetupG 'called only once Dim x(n), y(n) As Global Single Local i As Long interval = yn - y0 For i = 0 To n - 1 y(i) = y0 + interval * i / n x(i) = G(y(i)) Next y(n) = yn x(n) = G(yn) End Sub '---------------------------- Function Ginterpolate(y As Single) As Single Local k As Long k = Int((y - y0)*n / interval) Function = x(k) + (y - y(k)) * ( x(k+1) - x(k) ) * n / interval End Function '---------------------------- Function PBMain Local y As Single Local a As String y0 = 1 '<--- change to suit yn = 10 '<--- change to suit tol = .000001 'must be no smaller than precision of x SetupG Print " Range: "; Format$(y0);" to ";Format$(yn) Print " Number of interpolation points: "; Format$(n) Print " Tolerance: "; Format$(tol,"0.######") Print Do Line Input " Enter a number: "; a : If Trim$(a) = "" Then End y = Val(a) If y < y0 Or y > yn Then 'check range Print " Number must be between" + Str$(y0) + " and" + Str$(yn) Iterate End If Print " Interpolate:"; Ginterpolate(y) Print " Our inverse:"; G(y) '<--- not used in practice Print " Actual inverse:"; Sqr(y) '<--- not used in practice Print Loop End Function