Increasing Accuracy of Solution of Transcendental Equation

Approximate logarithm uses inaccurate approximation for log_2(1+x)-x but it still works

With some help from the comments I realized what my mistake was and adjusted the formula accordingly which leads to the expected result:

How approximation search works

Approximation search

This is analogy to binary search but without its restrictions that searched function/value/parameter must be strictly monotonic function while sharing the O(log(n)) complexity.

For example Let assume following problem

We have known function y=f(x) and want to find x0 such that y0=f(x0). This can be basically done by inverse function to f but there are many functions that we do not know how to compute inverse to it. So how to compute this in such case?

knowns

• y=f(x) - input function
• y0 - wanted point y value
• a0,a1 - solution x interval range

Unknowns

• x0 - wanted point x value must be in range x0=<a0,a1>

Algorithm

1. probe some points x(i)=<a0,a1> evenly dispersed along the range with some step da

   So for example x(i)=a0+i*da where i={ 0,1,2,3... }

2. for each x(i) compute the distance/error ee of the y=f(x(i))

   This can be computed for example like this: ee=fabs(f(x(i))-y0) but any other metrics can be used too.

3. remember point aa=x(i) with minimal distance/error ee

4. stop when x(i)>a1

5. recursively increase accuracy

   so first restrict the range to search only around found solution for example:

   a0'=aa-da;
   a1'=aa+da;
   

   then increase precision of search by lowering search step:

   da'=0.1*da;
   

   if da' is not too small or if max recursions count is not reached then go to #1

6. found solution is in aa

This is what I have in mind:

On the left side is the initial search illustrated (bullets #1,#2,#3,#4). On the right side next recursive search (bullet #5). This will recursively loop until desired accuracy is reached (number of recursions). Each recursion increase the accuracy 10 times (0.1*da). The gray vertical lines represent probed x(i) points.

Here the C++ source code for this:

//---------------------------------------------------------------------------
//--- approx ver: 1.01 ------------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _approx_h
#define _approx_h
#include <math.h>
//---------------------------------------------------------------------------
class approx
    {
public:
    double a,aa,a0,a1,da,*e,e0;
    int i,n;
    bool done,stop;

    approx()            { a=0.0; aa=0.0; a0=0.0; a1=1.0; da=0.1; e=NULL; e0=NULL; i=0; n=5; done=true; }
    approx(approx& a)   { *this=a; }
    ~approx()           {}
    approx* operator = (const approx *a) { *this=*a; return this; }
    //approx* operator = (const approx &a) { ...copy... return this; }

    void init(double _a0,double _a1,double _da,int _n,double *_e)
        {
        if (_a0<=_a1) { a0=_a0; a1=_a1; }
        else          { a0=_a1; a1=_a0; }
        da=fabs(_da);
        n =_n ;
        e =_e ;
        e0=-1.0;
        i=0; a=a0; aa=a0;
        done=false; stop=false;
        }
    void step()
        {
        if ((e0<0.0)||(e0>*e)) { e0=*e; aa=a; }         // better solution
        if (stop)                                       // increase accuracy
            {
            i++; if (i>=n) { done=true; a=aa; return; } // final solution
            a0=aa-fabs(da);
            a1=aa+fabs(da);
            a=a0; da*=0.1;
            a0+=da; a1-=da;
            stop=false;
            }
        else{
            a+=da; if (a>a1) { a=a1; stop=true; }       // next point
            }
        }
    };
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

This is how to use it:

approx aa;
double ee,x,y,x0,y0=here_your_known_value;
//            a0,  a1, da,n, ee  
for (aa.init(0.0,10.0,0.1,6,&ee); !aa.done; aa.step())
    {
    x = aa.a;        // this is x(i)
    y = f(x)         // here compute the y value for whatever you want to fit
    ee = fabs(y-y0); // compute error of solution for the approximation search
    }

in the rem above for (aa.init(... are the operand named. The a0,a1 is the interval on which the x(i) is probed, da is initial step between x(i) and n is the number of recursions. so if n=6 and da=0.1 the final max error of x fit will be ~0.1/10^6=0.0000001. The &ee is pointer to variable where the actual error will be computed. I choose pointer so there are not collisions when nesting this and also for speed as passing parameter to heavily used function creates heap trashing.

[notes]

This approximation search can be nested to any dimensionality (but of coarse you need to be careful about the speed) see some examples

• Approximation of n points to the curve with the best fit
• Curve fitting with y points on repeated x positions (Galaxy Spiral arms)
• Increasing accuracy of solution of transcendental equation
• Find Minimum area ellipse enclosing a set of points in c++
• 2D TDoA Time Difference of Arrival
• 3D TDoA Time Difference of Arrival

In case of non-function fit and the need of getting "all" the solutions you can use recursive subdivision of search interval after solution found to check for another solution. See example:

• Given an X co-ordinate, how do I calculate the Y co-ordinate for a point so that it rests on a Bezier Curve

What you should be aware of?

you have to carefully choose the search interval <a0,a1> so it contains the solution but is not too wide (or it would be slow). Also initial step da is very important if it is too big you can miss local min/max solutions or if too small the thing will got too slow (especially for nested multidimensional fits).

Related Topics