Numerical algorithm for finding the inverse of a function

Question

Is there a numerical method to approximate the inverse of a function for a given interval?

Thank you

Answer 1

Suppose that you have $y=f(x)$ and you want to have an approximation $x=g(y)$ valid over a rnage $ a \leq x \leq b$.

What you can do is to expand $f(x)$ as a Taylor series built around $x=c=\frac{a+b}2$ and get $$y=f(c)+ f'(c)(x-c)+\frac{1}{2} f''(c)(x-c)^2+\frac{1}{6} f'''(c) (x-c)^3+O\left((x-c)^4\right)$$ and then use series reversion to get $$x=c+\frac{1}{f'(c)}(y-f(c))-\frac{ f''(c)}{2 f'(c)^3}(y-f(c))^2+\frac{ 3 f''(c)^2-f'''(c) f'(c)}{6 f'(c)^5}(y-f(c))^3+O\left((y-f(c))^4\right)$$

Let us take the basic $f(x)=e^x$ with $a=-\frac 12$ and $b=\frac 12$ $$x=(y-1)-\frac{1}{2} (y-1)^2+\frac{1}{3} (y-1)^3+O\left((y-1)^4\right)$$ Now, for testing, assign a value to $x$ in order to get $y$ and recompute $x$ by the approximation formula. This will give $$\left( \begin{array}{cc} -0.5 & -0.491184 \\ -0.4 & -0.395969 \\ -0.3 & -0.298573 \\ -0.2 & -0.199684 \\ -0.1 & -0.099978 \\ +0.0 & +0.000000 \\ +0.1 & +0.100028 \\ +0.2 & +0.200511 \\ +0.3 & +0.302933 \\ +0.4 & +0.410535 \\ +0.5 & +0.529304 \end{array} \right)$$ For sure, with more terms in the first expansion the comparison will be batter.

Answer 2

You can also use simple reverse interpolation. You start with a table of points $(x_i, y_i)_{0\leq i \leq n}$, where $y_i=f(x_i)$, then you reverse the table and interpolate. For example, if you start with a table of values for $e^x$ $$ \left( \begin{array}{cc} x_i & y_i \\ -0.5 & 0.606531 \\ -0.4 & 0.67032 \\ -0.3 & 0.740818 \\ -0.2 & 0.818731 \\ -0.1 & 0.904837 \\ 0. & 1. \\ 0.1 & 1.10517 \\ 0.2 & 1.2214 \\ 0.3 & 1.34986 \\ 0.4 & 1.49182 \\ 0.5 & 1.64872 \\ \end{array} \right) \to \left( \begin{array}{cc} y_i & x_i \\ 0.606531 & -0.5 \\ 0.67032 & -0.4 \\ 0.740818 & -0.3 \\ 0.818731 & -0.2 \\ 0.904837 & -0.1 \\ 1. & 0. \\ 1.10517 & 0.1 \\ 1.2214 & 0.2 \\ 1.34986 & 0.3 \\ 1.49182 & 0.4 \\ 1.64872 & 0.5 \\ \end{array} \right), $$ The interpolating polynomial for the second table is $$ p(y)=-0.0959517 y^{10}+1.11934 y^9-5.89122 y^8+18.4819 y^7-38.4598 y^6+55.8896 y^5-58.1704 y^4+43.8458 y^3-24.158 y^2+10.4132 y-2.97439, $$ which provides an approximation for the inverse function in the region [0.606531, 1.64872]. Below you can see the graph of the interpolation error for the inverse.

In practice, you should probably use piecewise lower order interpolation, Chebyshev nodes or a combination of both.