Solving a special rational equation on a very small interval

Question

I need to solve the following equation (for $x$): \begin{equation} \mathcal f(x):=\sum_{i=1}^n b_i \left( \frac{a_i}{1+b_i x}\right)^2-\phi=0, \quad \text{with} \quad -1/b_1< x \le 0. \end{equation} Assume that the equation has at least one solution in $(-1/b_1,0]$ and $b_1\ge b_2\dots\ge b_n>0$. In practice, $b_1$ can be very large, so that the interval of interest $(-1/b_1,0]$ can be very small. Therefore, numerical issues can happen (I think because $1/(1+b_1x)$ can blow up, at least when a numerical method captures a solution that is very close to $-1/b_1$). What is the best practical idea to find a solution of $f(x)=0$ that falls into $(-1/b_1,0]$?

By the way, I am currently using the bisection method and $b_1$ is around 1000 for example. The quantity $n$ could be up to 50 in my case. Further, such an equation must be solved iteratively in my algorithm for different values of $ a_i$'s for $i\in [n]$. Thank you very much in advance for your time and help.

Answer 1

As said in comments, let $x=\frac y{b_1}$ and $c_i=\frac {b_i}{b_1}$ to make $$\sum_{i=1}^n b_i \left( \frac{a_i}{1+b_i x}\right)^2=b_1\Bigg[\frac{a_1^2}{(1+y)^2}+\sum_{i=2}^n \frac{ a_i^2\, c_i}{(1+ c_i\,y)^2}\Bigg]$$ Now, to get rid of the vertical asymptote at $y=-1$, we shall muliply everything by $(1+y)^2$ and consider the function

$$g(y)=b_1\Bigg[{a_1^2}+{(1+y)^2}\sum_{i=2}^n \frac{ a_i^2\, c_i}{(1+ c_i\,y)^2}\Bigg]-\phi(1+y)^2$$ for which $$g(0)=-\phi+b_1 \sum_{i=1}^n a_i^2\qquad \qquad g(-1)=b_1\,a_1^2 >0 \qquad \qquad g'(-1)=0$$ $$g''(-1)=-2\phi+2b_1 \sum_{i=2}^n \frac { a_i^2\, c_i } {(1-c_i)^2 }$$

Assuming $g(0)<0$, using Taylor around $y=-1$, we have $$g(y)=g(-1)+\frac 12 g''(-1) (y+1)^2+ O\left((y+1)^3\right)$$ and then a first estimate $$y_0=-1 +\sqrt {-\frac {b_1 a_1^2} {g''(-1)}}$$

Let us try using $n=5$, $b_i=10^{5-i}$, $a_i=p_{i+5}^2$ and $\phi=12345678987654321$. This gives as an estimate $$y_0=-1+\sqrt{\frac{175914852251400405000}{15208068862859378242920969041} }\sim -0.999892$$ while the exact solution, obtained using Newton method is $-0.999848$

Starting with this guess, Newton iterates would be $$\left( \begin{array}{cc} n & y_n \\ 0 & \color{red}{-0.99989}244905795162394 \\ 1 & \color{red}{-0.9998}3867358692746307 \\ 2 & \color{red}{-0.999847}63616543148912 \\ 3 & \color{red}{-0.999847899}77068160754 \\ 4 & \color{red}{-0.999847899999109}03399 \\ 5 & \color{red}{-0.99984789999910920552} \end{array} \right)$$

Very immodestly, let me precise that this method bears my name.

Edit

If you are coding, for a complete safety, combine bisection and Newton methods . If you go to this place, on page $359$, you will find subroutine RTSAFE (here and here) which does exactly that. It is very robust.

I apologize for giving you a reference to Fortran coding. If you can access the books of Numerical Recipes, you will find the equivalent in C and C++.

Answer 2

I prefer to add a second answer for some crucial points.

The solution you are looking for is not in the range $-\frac 1{b_1} < x <0$ but in the range $$\color{red}{-\frac 1{b_1}+\epsilon \leq x \leq-\frac 1{b_2}-\epsilon}$$ $\epsilon$ being at the level of machine accuracy. So, if your search method use bounds, these are the ones to be used.

Moreover, in this range, there are two possible solutions for $x$ since $$\lim_{x\to {-\frac{1}{b_1}}+\epsilon} \, f(x) = +\infty \qquad \text{and} \qquad \lim_{x\to {-\frac{1}{b_2}}-\epsilon} \, f(x)= +\infty$$but, depending on the value of $\phi$, no root is also possible as well as a double root.

Reusing the example worked in my previous answer, if $\phi < 2.57\times 10^8$, there is no solution.

Using $\phi=10^{10}$, the roots are $-9.085\times 10^{-4}$ and $-1.170\times 10^{-4}$ each of them being close to the bounds.

So, considering before any rescaling $$f(x)=\sum_{i=1}^n b_i \left( \frac{a_i}{1+b_i x}\right)^2- \phi$$ we should consider $$g(x)=(1+b_1x)^2(1+b_2x)^2\sum_{i=1}^n b_i \left( \frac{a_i}{1+b_i x}\right)^2- \phi(1+b_1x)^2(1+b_2x)^2$$ that is to say $$g(x)=a_1^2b_1(1+b_2x)^2+a_2^2b_2(1+b_1x)^2+(1+b_1x)^2(1+b_2x)^2\sum_{i=3}^n b_i \left( \frac{a_i}{1+b_i x}\right)^2-\phi(1+b_1x)^2(1+b_2x)^2$$

So, working around $x=-\frac{1}{b_1}$, we have $$g\left(-\frac{1}{b_1}\right)=a_1^2 b_1 \left(1-\frac{b_2}{b_1}\right)^2$$ $$g'\left(-\frac{1}{b_1}\right)=2 a_1^2 b_1 b_2 \left(1-\frac{b_2}{b_1}\right)$$ $$g''\left(-\frac{1}{b_1}\right)=2 b_1 b_2 \left(a_1^2 b_2+a_2^2 b_1\right)-2 b_1^2 \left(1-\frac{b_2}{b_1}\right)^2 \phi+$$ $$2 b_1^2 \left(1-\frac{b_2}{b_1}\right)^2\sum_{i=3}^n \frac{a_i^2 b_i}{\left(1-\frac{b_i}{b_1}\right)^2}$$

Now, from these informations, build the parabola $$h(x)=A+Bx+Cx^2$$ $$C=\frac 12 g''\left(-\frac{1}{b_1}\right)\qquad B=g'\left(-\frac{1}{b_1}\right)+\frac{2C}{b_1}\qquad A=g\left(-\frac{1}{b_1}\right)+\frac{B}{b_1}-\frac{C}{b_1^2}$$ Solve the quadratic $h(x)=0$ and keep the largest root as the approximate solution.

Applied to the worked example, this would give $$x_0=-0.9998478974\times 10^{-4}$$ while the exact solution is $$x=−0.9998479000\times 10^{-4}$$

Update

In the previous edit, there were $b_1$ and $b_2$ missing as mutliplying factors in the expansion of $g(x)$ but all the remaining was correct.

To allow test with the conditions $n=5$, $b_i=10^{5-i}$, $a_i=p_{i+5}^2$ and $\phi=10^k$, I produced some results $$\left( \begin{array}{ccc} k & \text{estimation} & \text{solution} \\ 10 & -8.26703656651674 \times 10^{-5} & -8.30010445244925 \times 10^{-5} \\ 11 & -9.46205809895313 \times 10^{-5} & -9.46525763471424 \times 10^{-5} \\ 12 & -9.83067188617129 \times 10^{-5} & -9.83098989883711 \times 10^{-5} \\ 13 & -9.94652543341715 \times 10^{-5} & -9.94655718742864 \times 10^{-5} \\ 14 & -9.98309681582244 \times 10^{-5} & -9.98309998987002 \times 10^{-5} \\ 15 & -9.99465543307055 \times 10^{-5} & -9.99465575043391 \times 10^{-5} \\ 16 & -9.99830996825483 \times 10^{-5} & -9.99830999998987 \times 10^{-5} \\ 17 & -9.99946557190165 \times 10^{-5} & -9.99946557507511 \times 10^{-5} \\ 18 & -9.99983099968264 \times 10^{-5} & -9.99983099999999 \times 10^{-5} \\ 19 & -9.99994655747581 \times 10^{-5} & -9.99994655750754 \times 10^{-5} \\ 20 & -9.99998309999683 \times 10^{-5} & -9.99998310000000 \times 10^{-5} \end{array} \right)$$

Edit

If we are sure that $x$ will be very close to $-\frac 1{b_1}$, then, we can write $x= -\frac 1{b_1}+\epsilon$ and expand $f(x)$ as a series around $\epsilon=0$.

Defining $$u_i=\frac{a_i^2 b(i)}{(1+ b_ix)^2}$$ we have $$u_1=\frac{a_1^2}{b_1 \epsilon ^2}$$ and, for $i>1$, $$u_i=b_1^2\sum_{k=0}^\infty (-1)^k (k+1)\frac{a_i^2\, b_1^{k}\, b_i^{k+1} } {(b_1-b_i)^{k+2} }\epsilon^k$$ Neglecting the positive powers of $\epsilon$ then leads to the approximate equation $$\frac{a_1^2}{b_1 \epsilon ^2}+b_1^2 \sum _{i=2}^n \frac{a_i^2 b_i}{(b_1-b_i)^2}=\phi$$

This would give as a first approximation $$ x=-\frac 1{b_1}+\frac{a_1}{\sqrt{b_1}}\left(\phi -b_1^2 \sum _{i=2}^n \frac{a_i^2 b_i}{(b_1-b_i)^2}\right)^{-\frac 12} $$

Using it for the above conditions, for $\phi=10^{10}$ this gives $x=--8.2997747\times 10^{-5}$ and for $\phi=10^{20}$ this gives $x=-9.9999831\times 10^{-5}$; this is much better than the parabolic approximation.

Answer 3

The current framing of the Question suggests that the $b_i$'s are fixed but the problem (root-finding) will need to be solved repeatedly for varying coefficient's $a_i$. My remarks are intended to help with the selection of an initial approximation of the root, to be improved iteratively.

Note that the problem specifies only that the $b_i$'s are ordered by weak inequalities:

$$ b_1 \ge b_2 \ge b_3 \ge \ldots \ge b_n $$

Unless the intention is vary those $b_i$'s, we can eliminate a priori consecutive terms with equal $b_i$'s (by lumping them together) in the summation that defines $f(x)$:

$$ f(x) := \sum_{i=1}^n b_i \left( \frac{a_i}{1+b_i x} \right)^2 - \phi $$

So let's assume the $b_i$'s are distinct, although to be sure there is no requirement that they are widely spaced.

Then all the terms in the summation are continuous and strictly decreasing on the interval $(-1/b_1,0]$ in which the root is sought, and only the first term fails to be continuous at left hand endpoint $x=-1/b_1$ because the first term tends to $+\infty$ there.

It follows that any root $f(r) = 0$ with $r\in (-1/b_1,0]$ is unique, and such a root exists if and only if $f(0) \le 0$. The case $f(0)=0$ is easily checked as it amounts to $\sum a_i^2 b_i = \phi$, so we assume hereafter that the check shows $f(0) \lt 0$.

Locating $r$ is striking the right balance between the contributions of the first term, which are "potentially" unbounded, and the contributions of the rest of the summands, which are at most a finite value we can directly evaluate at $x=-1/b_1$.

My suggestion is that we evaluate those terms:

$$ g(x) := \sum_{i=2}^n \frac{a_i^2 b_i}{(1 + b_i x)^2} $$

at two points and use linear interpolation to approximate their sum in between. Since those terms are concave up (a convex function), that linear interpolation is consistently an overestimate of their true contribution.

Initially we would evaluate $g(x)$ at the two points $x=-1/b_1$ and $x=0$. Then we would solve the simpler rational equation (as equivalent to a cubic):

$$ \frac{a_1^2 b_1}{(1 + b_1 x)^2} + b_1(g(0) - g(-1/b_1)) x + g(0) = \phi $$

The solution $r_1$ of this equation will yield $f(r_1) \lt 0$ because of the error introduced by linear interpolation. Therefore $r_1 \gt r$ overestimates the actual root.

We can then repeat this procedure, replacing evaluation of $g$ at the upper endpoint $x=0$ with evaluation at new upper endpoint $x=r_1$. The attraction of this method is that it conservatively approaches the actual root from above, and therefore does not risk passing the singularity at $x = -1/b_1$.