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;
}

I'm required to find all the roots at once, and it was one of the methods suggested on wikipedia http://en.wikipedia.org/wiki/Root-finding_algorithm#Finding_all_roots_at_once I have also used it before for a polynomial (it worked the previous time) so I could reuse the code — Ca01an, Dec 14 '14 at 20:47
No method was suggested, although it is an extension of a previous question which used the newton-raphson method, but I dont think newtons method will work for this. Not familiar with the bisection method but will look it up. — Ca01an, Dec 14 '14 at 21:01
Newton's method would work perfectly fine, especially if you pre-screen to find good first guesses. My posted answer does not use that method, however. EDIT: Newton's method will also has a higher rate of convergence than the method I posted. — druckermanly, Dec 14 '14 at 21:06
Also, please accept answers that have helped you or helped solve your problem. — druckermanly, Dec 14 '14 at 21:18
My problem is that I've used newtons method to find one root at a time, but I'm not sure how to get it to find all the roots at once. — Ca01an, Dec 14 '14 at 21:19
Do you see how I got the values for `relevant_intervals`? You can use the midpoint of each of the values in `relevant_intervals` as your initial guess in Newton's method, and then you would (hopefully) get all the distinct answers. — druckermanly, Dec 14 '14 at 21:21

score 0 · Answer 1 · answered Dec 14 '14 at 20:06

I am not familiar with Durand-Kerner method, but according to Wiki http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method, it is suitable only for solving polynomial equations. Note the line in the wiki page: "The explanation is for equations of degree four. It is easily generalized to other degrees."

Your equation is not polynomial. The numerical solution will probably not converge.

Regardless of that your function f returns wrong formula: return a4 * x * exp(-abs(x)) * cos(x) + a0; (you forgot about complex modulo, i.e. abs)

and your iterations seem also wrong. They should read:

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

but even if you make these changes, the solution will not converge - it will be oscilating. The reason is probably as written above - the method is not suitable for this type of equation.

Thanks, any ideas what other method I could use? – Ca01an Dec 14 '14 at 20:35 — Ca01an, Dec 14 '14 at 20:35

druckermanly · Accepted Answer · 2014-12-14T22:57:55.037

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;
    }
}

Finding the roots of a nonlinear equation using C++

2 Answers2