Newton method for computing an inverse

Question

This is following the question I asked in this thread : Link error missing vtable

I defined a class 'function' and two others classes 'polynomial' and 'affine' that inherit from 'function'.

class function {

        public:
            function(){};
            virtual function* clone()const=0;
            virtual float operator()(float x)const=0; //gives the image of a number by the function
            virtual function* derivative()const=0;
            virtual float inverse(float y)const=0; 
            virtual ~function(){}

        };

        class polynomial : public function {
        protected:
            int degree;
        private:
            float *coefficient;
        public:

            polynomial(int d);
            virtual~polynomial();
            virtual function* clone()const;
            int get_degree()const;
            float operator[](int i)const; //reads coefficient number i
            float& operator[](int i); //updates coefficient number i
            virtual float operator()(float x)const; 
            virtual function* derivative()const;
            virtual float inverse(float y)const; 

        };

        class affine : public polynomial {
            int a; 
            int b; 
            //ax+b
    public:
            affine(int d,float a_, float b_);
            function* clone()const;
            float operator()(float x)const;
            function* derivative()const;
            float inverse(float y)const;
            ~affine(){}
        };

Method inverse in polyomial does not seem to work fine. It is based on the Newton method applied to the function x->f(x)-y for fixed y (the element for which we're computing the inverse) and the current polynomial f.

float polynomial::inverse(float y)const
{
    int i=0;
    float x0=1;
    function* deriv=derivative();
    float x1=x0+(y-operator()(x0))/(deriv->operator()(x0));
    while(i<=100 && abs(x1-x0)>1e-5)
    {
        x0=x1;
        x1=x0+(y-operator()(x0))/(deriv->operator()(x0));
        i++;
    }

    if(abs(x1-x0)<=1e-5) 
    {
        //delete deriv; //I get memory problems when I uncomment this line
        return x1;   
    }

    else    
    {
        cout<<"Maximum iteration reached in polynomial method 'inverse'"<<endl;
        //delete deriv; //same here
        return -1;
    }
}

double polynomial::operator()(double x)const
{
    double value=0;
    for(int i=0;i<=degree;i++) value+=coefficient[i]*pow(x,i);
    return value;
}
polynomial* polynomial::derivative()const
{
    if(degree==0)
    {
        return new affine(0,0,0);
    }
    polynomial* deriv=new polynomial(degree-1);
    for(int i=0;i<degree;i++)
        deriv[i]=(i+1)*coefficient[i+1];   
    return deriv;
}

I test this method with p:x->x^3 :

#include "function.h"

int main(int argc, const char * argv[])
{

        polynomial p(3);
        for(int i=0;i<=2;i++) p[i]=0;
        p[3]=1;
        cout<<"27^(1/3)="<<p.inverse(27);



        return 0;
}

This script outputs 27^(1/3)=Maximum iteration reached in polynomial method 'inverse' -1 even if I put 10,000 instead of 100. I've read some articles on the internet and it seems that it's a common way to compute the inverse.

Answer 1

Well, the problem was in method 'derivative'. Instead of using the 'operator[]' that I redefined, I used '->coefficient[]' and the main script worked fine for p.inverse(27) (only 14 iterations). I just replaced deriv[i]=(i+1)*coefficient[i+1]; with deriv->coefficient[i]=(i+1)*coefficient[i+1];

Answer 2

the abs function prototype is: int abs(int) So a test like abs(x1-x0)<=1e-5 won't behave as you expect; you compare a int with a float. In this case the float will be converted to int so it the same as abs(x1-x0)<=0

This is probably why you don't get the expected result - I suggest adding a few more printouts to get to the bottom of things.

Answer 3

Check This Code :

#include<iostream>
#include<cmath>
#include<math.h>

using namespace std;

void c_equation(int choose, double x);
void Processes(double x, double fx1, double fdx1, int choose);

void main()
{
    int choose,choose2;
    double x;
    system("color A");
    cout << "                                                             " << endl;
    cout << "=============================================================" << endl;
    cout << "Choose Equation : " << endl;
    cout << "_____________________________________" << endl;
    cout << "1- x-2sin(x)" << endl;
    cout << "2- x^2 + 10 cos(x)" << endl;
    cout << "3- e^x - 3x^2" << endl;
    cout << "                                                    " << endl;
    cin >> choose;
    cout << "If you have values  press 1/ random press 2 :" << endl;
    cin >> choose2;
    if (choose2 == 1)
    {
        cout << "                                         " << endl;
        cout << "Enter Xo : " << endl;
        cin >> x;
        c_equation(choose, x);
    }
    else if (choose2 == 2)
    {
        x = rand() % 20;
        cout << "Xo = " << x << endl;
        c_equation(choose, x);
        choose2 = NULL;
    }
    else
    {
        cout << "Worng Choice !! " << endl;
        choose = NULL;
        choose2 = NULL;
        main();
    }
}

void c_equation(int choose, double x)
{
    double fx;
    double fdx;
    double fddx;
    double result;
    if (choose == 1)
    {
        fx = x - 2 * sin(x);
        fdx = 1 - 2 * cos(x);
        fddx = 2 * sin(x);
        result = abs((fx * fddx) / pow(fdx, 2));
    }

    else if (choose == 2)
    {

        fx = pow(x, 2) + 10 * cos(x);
        fdx = 2 * x - 10 * sin(x);
        fddx = 2 - 10 * cos(x);
        result = abs((fx * fddx) / pow(fdx, 2));
    }

    else if (choose == 3)
    {
        fx = exp(x) - 3 * pow(x, 2);
        fdx = exp(x) - 6 * x;
        fddx = exp(x) - 6;
        result = abs((fx * fddx) / pow(fdx, 2));
    }
    else
    {
        cout << " " << endl;
    }
    //------------------------------------------------------------
    if (result < 1)
    {
        cout << "True Equation :) " << endl;
        Processes(x, fx, fdx , choose);
    }
    else
    {
        system("cls");
        cout << "False Equation !!" << endl;
        choose = NULL;
        x = NULL;
        main();
    }

}
void Processes(double x, double fx, double fdx , int choose)
{
    double xic;
    for (int i = 0; i < 3; i++)
    {
        xic = x - (fx / fdx);
        cout << "                                          " << endl;
        cout << "Xi = " << x << "   " << "F(Xi) = " << fx   << "   " << "  F'(Xi) = " << fdx   << "  " << "  Xi+1 = " << xic << endl;
        x = xic;
        if (choose == 1)
        {
            fx = xic - 2 * sin(xic);
            fdx = 1 - 2 * cos(xic);
        }
        else if (choose == 2)
        {
            fx = pow(xic, 2) + 10 * cos(xic);
            fdx = 2 * xic - 10 * sin(xic);
        }
        else if (choose == 3)
        {
            fx = exp(xic) - 3 * pow(xic, 2);
            fdx = exp(xic) - 6 * xic;
        }

    }
}