Finding the roots of a nonlinear equation using C++

Question

What method would I use to find the roots of f(x) = 5x(e^-mod(x))cos(x) + 1 ? I have being trying the durand-kerner method but I can't get it to work. Are there any easier ways of doing it?

Here is the my code using the durand-kerner method

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

using namespace std;

typedef complex<double> dcmplx;

dcmplx f(dcmplx x)
{
    // the function we are interested in
    double a4 = 5;

    double a0 = 1;

    return (a4 * x * exp(-x) * cos(x) )+ a0;
}


int main()
{

dcmplx p(.9,2);
dcmplx q(.1, .5);
dcmplx r(.7,1);
dcmplx s(.3, .5);

dcmplx p0, q0, r0, s0;

int max_iterations = 100;
bool done = false;
int i=0;

while (i<max_iterations && done == false)
{
    p0 = p;
    q0 = q;
    r0 = r;
    s0 = s;


p = p0 - f(p0)/((p0-q)*(p0-r)*(p0-s));
q = q0 - f(q0)/((q0-p)*(q0-r)*(q0-s));
r = r0 - f(r0)/((r0-p)*(r0-q)*(r0-s0));
s = s0 - f(s0)/((s0-p)*(s0-q)*(s0-r));

    // if convergence within small epsilon, declare done
    if (abs(p-p0)<1e-5 && abs(q-q0)<1e-5 && abs(r-r0)<1e-5 && abs(s-s0)<1e-5)
        done = true;

    i++;
}

cout<<"roots are :\n";
cout << p << "\n";
cout << q << "\n";
cout << r << "\n";
cout << s << "\n";
cout << "number steps taken: "<< i << endl;

return 0;
}

Answer 1

This approach makes use of the bisection method, and the fact that you can do a little math to find an upper bound for the highest zero, respectively.

Reproduced at http://ideone.com/fFLjsh

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <utility>

#define MINX (-20.0f)
//MAXX Happens to be an upper bound for all zeroes of the function...
#define MAXX (1.0f)
#define NUM_INTERVALS (1000000)
#define NUM_BISECTION_ITERATIONS (30)

using namespace std;

double f(double x){
    return 5 * x * exp(-x) * cos(x) + 1;
}

double bisection_method(double x0, double x1){
    for (unsigned int i = 0; i < NUM_BISECTION_ITERATIONS; i++){
        double midpoint = 0.5*x0 + 0.5*x1;
        f(x0) * f(midpoint) < 0 ? x1 = midpoint : x0 = midpoint;
    }
    return 0.5*x0 + 0.5*x1;
}

int main(int argc, char** argv){
    vector<pair<double, double>> relevant_intervals;
    for (unsigned int i = 0; i < NUM_INTERVALS - 1; i++){
        double x0 = MINX + (MAXX - MINX) / NUM_INTERVALS * (i);
        double x1 = MINX + (MAXX - MINX) / NUM_INTERVALS * (i + 1);
        if (f(x0) * f(x1) < 0)
            relevant_intervals.push_back(make_pair(x0, x1));
    }
    cout << fixed << setprecision(20);
    for (const auto & x : relevant_intervals){
        cout << "One solution is approximately " << bisection_method(x.first, x.second) << endl;
    }
}