Efficient approximation method

Question

We have an equation of the form: $P(x)=(1-xv)(1-xw)(1-x)^{n-2}-u$

$n≥2$ and $n$ belongs to the set of positive integers.

$0<u<1$.

$0<v≤1$.

$0<w≤1$.

$0<x≤1$. ($x>1$, when $n=2$ and $v,w<<1)$ which is very rare.

For this equation it seems, a closed form expression for $x$ is not possible, so what would be the best, most efficient method of approximation, Newton's method requires a lot of iterations and often it doesn't converges, so we need a method that's more robust, requires fewer iterations and would most probably converge, it doesn't matter how rigourous the method is. PS: We can actually make a pretty good guess of value of $x$ but it'll always be smaller than the accurate value.

Answer 1

The equation $(1-xv)(1-xw)(1-x)^{n-2}-u = 0$ can be expressed as

$$x = 1 - \left( \frac{u}{(1-xv)(1-xw)} \right)^{\frac{1}{n-2}}$$

Now, that's not a closed-form solution for $x$, since it still occurs (twice) on the right side, but it might be viable as a fixed-point iteration. Let's try it.

def solve_x(u, v, w, n):
    # Gotta make some initial guess.
    x = 1.0
    # Arbitarily cap the number of iterations.
    for k in range(60):
        old_x = x
        x = 1 - (u / (1 - x * v) * (1 - x * w)) ** (1 / (n - 2))
        # No point in continuing after we get stuck on one root.
        if x == old_x:
            break
        # Tracing intermediate x values for demonstration purposes.
        print('After', k + 1, 'iterations, x =', x)    
    return x

A couple of examples:

>>> solve_x(u=0.1, v=0.2, w=0.3, n=4)
After 1 iterations, x = 0.7041960108450191
After 2 iterations, x = 0.6970087853011551
After 3 iterations, x = 0.6968484196253488
After 4 iterations, x = 0.696844848525994
After 5 iterations, x = 0.6968447690065588
After 6 iterations, x = 0.696844767235862
After 7 iterations, x = 0.6968447671964331
After 8 iterations, x = 0.6968447671955551
After 9 iterations, x = 0.6968447671955355
After 10 iterations, x = 0.696844767195535
0.696844767195535
>>> solve_x(u=0.5, v=0.6, w=0.7, n=8)
After 1 iterations, x = 0.1508093352175236
After 2 iterations, x = 0.11158050977372924
After 3 iterations, x = 0.11088585440922338
After 4 iterations, x = 0.1108738930488482
After 5 iterations, x = 0.11087368718522883
After 6 iterations, x = 0.11087368364219763
After 7 iterations, x = 0.11087368358122007
After 8 iterations, x = 0.11087368358017058
After 9 iterations, x = 0.11087368358015248
After 10 iterations, x = 0.11087368358015215
0.11087368358015215

Convergence after 10 iterations. Not bad. But maybe some tweaking to the algorithm would make it converge even faster.

Answer 2

$$P(x)=(1-x\,v)(1-x\,w)(1-x)^{n-2}-u$$ is bounded by $$Q(x)=(1-x\, 1)(1-x\, 1)(1-x)^{n-2}-u=(1-x)^{n}-u\implies x_Q=1-u^{\frac{1}{n}}$$and $$R(x)=(1-x \,0)(1-x\,0)(1-x)^{n-2}-u=(1-x)^{n-2}-u\implies x_R=1-u^{\frac{1}{n-2}}$$

So, being lazy, draw secant between the two points; this defines $$x_P= \frac { x_R\,P(x_Q)-x_Q\,P(x_R) }{P(x_Q)-P(x_R)}$$ and use $x_P$ as the $x_0$ of Newton method.

Using the same values as @Dan in his/her answer, we have for the first case $(u=0.1,v=0.2,w=0.3,n=4)$

$$\left( \begin{array}{cc} n & x_n \\ 0 & 0.6413440156 \\ 1 & 0.6244679503 \\ 2 & 0.6249584091 \\ 3 & 0.6249588304 \end{array} \right)$$ and for the second case $(u=0.5,v=0.6,w=0.7,n=8)$

$$\left( \begin{array}{cc} n & x_n \\ 0 & 0.0913449582 \\ 1 & 0.0908387018 \\ 2 & 0.0908395909 \end{array} \right)$$

I am almost sure that, being less lazy, we could do better.

Notice that the required derivative is very inexpensive using logarithmic differentiation $$\frac{P'(x)}{P(x)}=-\frac{n-2}{1-x}-\frac{v}{1-x\,v }-\frac{w}{1-x\,w }$$

Edit

Being less lazy, using the fact that we know three points $$(x_Q,f(x_Q)) \qquad \qquad(x_R,f(x_R))\qquad\qquad (x_P,f(x_P))$$compute the equation of the parabola going through them, solve the quadratic and keep as estimate the root which is $\in (x_R,x_Q)$ and use the result as $x_0$.

Fot the first example, it gives $x_0= 0.62510$

Answer 3

I prefer to add another answer because I do not know yet what could be the limitations of such an approach.

Assume that we can write $$(1-x\,v)(1-x\,w)(1-x)^{n-2}\sim (1-x)^m$$

Expanding as series around $x=0$ $$\Delta=(1-x\,v)(1-x\,w)(1-x)^{n-2}-(1-x)^m$$ $$\Delta=(2 + m - n - v - w)\, x+O(x^2)\quad \implies m=n+u+w-2$$

Using this $m$ $$\Delta=\frac 12 [v(1-v)+w(1-w)]\,x^2+O(x^3)$$

Applied to the examples used in the previous answers

• for the first case : $m=2.3$ which leads to $x_0=0.632534$ (the exact solution being $0.624959$)
• for the second case : $m=7.3$ which leads to $x_0=0.090583$ (the exact solution being $0.090840$)

and this did not require the computation of any of the terms $x_R$, $f(x_R)$, $x_Q$, $f(x_Q)$, $x_P$.

Generating random numbers, some results (which look, at least to me, more than promizing)

$$\left( \begin{array}{cccccc} u & v & w & n & \text{first estimate} & \text{exact solution} \\ 0.48 & 0.84 & 0.46 & 43 & 0.0173482 & 0.0172032 \\ 0.19 & 0.44 & 0.96 & 25 & 0.0664561 & 0.0658234 \\ 0.87 & 0.33 & 0.53 & 54 & 0.0026045 & 0.0026311 \\ 0.23 & 0.84 & 0.41 & 74 & 0.0200290 & 0.0198649 \\ 0.06 & 0.90 & 0.71 & 79 & 0.0355295 & 0.0351589 \\ 0.88 & 0.46 & 0.57 & 18 & 0.0072989 & 0.0074791 \\ 0.94 & 0.62 & 0.60 & 31 & 0.0020240 & 0.0020454 \\ 0.17 & 0.49 & 0.59 & 58 & 0.0307361 & 0.0305705 \\ 0.48 & 0.43 & 0.81 & 32 & 0.0231840 & 0.0232241 \\ 0.57 & 0.68 & 0.06 & 37 & 0.0156528 & 0.0156059 \\ 0.03 & 0.71 & 0.55 & 89 & 0.0392474 & 0.0389549 \\ 0.10 & 0.99 & 0.82 & 41 & 0.0560480 & 0.0548658 \\ 0.88 & 0.71 & 0.17 & 12 & 0.0115020 & 0.0116828 \\ 0.64 & 0.33 & 0.30 & 92 & 0.0048955 & 0.0049122 \\ 0.84 & 0.57 & 0.36 & 51 & 0.0034672 & 0.0034859 \\ 0.91 & 0.56 & 0.53 & 38 & 0.0025158 & 0.0025396 \\ 0.01 & 0.93 & 0.59 & 45 & 0.1002360 & 0.0983040 \\ 0.32 & 0.97 & 0.41 & 41 & 0.0282721 & 0.0278260 \\ 0.28 & 0.12 & 0.14 & 100 & 0.0128510 & 0.0128717 \\ 0.51 & 0.17 & 0.41 & 82 & 0.0082866 & 0.0083216 \end{array} \right)$$

Edit

In view of these results, we can even do better.

Compute $$\Phi(m)=\int_0^1\Big[(1-x\,v)(1-x\,w)(1-x)^{n-2}-(1-x)^m\Big]^2 \,dx$$ which does not make any problem. $\Phi'(m)=0$ implies the nasty equation in $m$ $$\frac{2 v w \left(3 (m+n)^2-1\right)-v (m+n+1)^2 (2 m+2 n-1)}{(m+n)^2 (m+n+1)^2}+$$ $$\frac{w (-2 m-2 n+1)}{(m+n)^2}+\frac{(m-n+2) (3 m+n)}{(2 m+1)^2}=0$$

Making $m=n+u+w-2+\epsilon$ and expanding as a series around $\epsilon=0$ to first order gives a longish but fully explicit expression $$\epsilon(v,w,n)=\frac A {2B}$$ where $$A=\frac{1-v}{2 n+v+w-3}+\frac{2 v-1}{2 n+v+w-2}-\frac{v}{2 n+v+w-1}+\frac{(v-1) (2 n+v-2)}{(2 n+v+w-3)^2}-$$ $$\frac{2 \left(n (2 v-1)+(v-1)^2\right)}{(2 n+v+w-2)^2}+\frac{v (2 n+v-1)}{(2 n+v+w-1)^2}+\frac{1}{(2 n+2 v+2 w-3)^2}$$ and $$B=\frac{1-v}{(2 n+v+w-3)^2}+\frac{2 v-1}{(2 n+v+w-2)^2}-\frac{v}{(2 n+v+w-1)^2}+\frac{(v-1) (2 n+v-2)}{(2 n+v+w-3)^3}-\frac{2 \left(n (2 v-1)+(v-1)^2\right)}{(2 n+v+w-2)^3}+\frac{v (2 n+v-1)}{(2 n+v+w-1)^3}+\frac{2}{(2 n+2 v+2 w-3)^3}$$

For the first terms of the above table, the first estimate is now

$$\{0.017205,0.065821,0.002631,0.019865,0.035158,0.007488,0.002046\}$$ which is much better.

Update

Working on a logarithmic scale, there is an explicit expression of $m$.

Minimizing $$\Psi(m)=\int_0^1 \Big[\log(1-x\,v)+\log(1-x\,w)+(n-2-m)\log(1-x) \Big]^2\,dx$$ leads to $$2(m-n)=\frac{1-v}v \left(\text{Li}_2\left(\frac{v}{v-1}\right)+\log (1-v)\right)+$$ $$\frac{1-w} w \left(\text{Li}_2\left(\frac{w}{w-1}\right)+\log (1-w)\right)$$

The function $$g(t)=\frac{1-t} t \left(\text{Li}_2\left(\frac{t}{t-1}\right)+\log (1-t)\right)$$ is very well approximated by its $[2,2]$ Padé approximant built around $t=\frac 12$.

If we make the very crude approximation $$g(t)=-\left(\frac{\pi ^2}{12}+\log (2)\right)+\left(\frac{\pi ^2}{3}-2\right) \left(t-\frac{1}{2}\right)+O\left(\left(t-\frac{1}{2}\right)^2\right)$$

$$m=n +\left(\frac{\pi ^2}{6}-1\right) (v+w)+\left(1-\frac{\pi ^2}{4}-\log (2)\right)$$