0

I was asked to code a root finding program which uses the secant method. Although the exercise itself isn't any hard I am passing only 2 of 5 declared examples and I have no idea why. My function looks like this: 

template <class T >
bool secant ( double &x , T f , double x0 , double x1 , double eps , int n = 1000){
    x = x0;
    double c, xp, xInter, xm;
    if(eps<=0){cout << "Wrong epsilon" << endl; return false;}
    if(n<=0){cout << "Wrong n" << endl; return false;}
    if(f==nullptr){cout <<"No function" << endl; return false;}
    if (f(x0)*f(x1)<0){
        do{
            //intermediate value
            xInter = (x0*f(x1)-x1*f(x0))/(f(x1)-f(x0));
            cout << n << endl;
            c = f(xInter)*f(x0);
            x0 = x1;
            x1 = xInter;
            n--;
            if (c==0){
                break;
            }
            if (n==0){
                cout << "Number of iterations exceeded.";
                break;
            }
            xm = (x0*f(x1)-x1*f(x0))/(f(x1)-f(x0));
        } while(abs(xm - xInter)>= eps);
    }
    x = xInter;
    return true;
}

And is called in main like so:

double x;
    for ( auto fx : { f1 , f2 , f3 , f4 , f5 }) {
        for ( auto eps : {0.1 , 0.01 , 0.001 , 0.0001 , 0.0000001}) {
            if ( secant (x , fx , 0.1 , 3 , eps )) {
                cout << setprecision (8) << "eps = " << eps
                << "\t root = " << x << endl;
            } else {
                cout << " Unable to find root " <<endl;
            }
        }

        cout << endl ;
    }

There are 5 functions to test:

  1. x^2
  2. x^2-2
  3. e^x+x-1
  4. log(x)-x+1
  5. sin(x+0.5)

Both 2nd and 5th functions get good results (1.4142136 and 2.6415926, as expected) but every other function returns the same numbers being:

  1. (for the epsilon = 0.1): 8.0...e-307
  2. (for every other epsilon): 1.5...e+231

What am I doing wrong?

Jakub Sapko
  • 304
  • 2
  • 15
  • Two things that jumped out at me: you use `abs` in your `while`, which is an integer function. Try `fabs` for floating point numbers. Smaller and probably unimportant, you have a check of `c==0` where c is a double. I'd change that to `==0.0` – possum May 06 '20 at 13:32
  • Why do you compute an additional midpoint `xm`? What exactly is the role of `c`? A sensible break condition based on the function value would be `fabs(f(xInter)) < eps*fabs(f(x0))`. – Lutz Lehmann May 06 '20 at 15:04
  • Why do you think that the method should converge in all examples, esp. in examples 1 and 3 that are positive except at the double root `x=0`? The secant method is rather fragile, it is fast **if** it converges, but goes astray in more circumstances than the Newton method. Are you sure your method is secant and not regula falsi in one of its variants or Dekker's method? – Lutz Lehmann May 07 '20 at 09:57

0 Answers0