CODEWITHSUNDEEP
c++calgorithmmath
shybovycha
shybovycha

Reputation: 12265

Simple iteration algorithm

If we are given with an array of non-linear equation coefficients and some range, how can we find that equation's root within the range given?

E.g.: the equation is

enter image description here

So coefficient array will be the array of a's. Let's say the equation is

enter image description here

Then the coefficient array is { 1, -5, -9, 16 }.

As Google says, first we need to morph function given (the equation actually) to some other function. E.g. if the given equation is y = f(x), we should define other function, x = g(x) and then do the algorithm:

while (fabs(f(x)) > etha)
  x = g(x);

To find out the root.

The question is: how to define that g(x) using coefficient array and the range given only?

The problem is: when i define g(x) like this

enter image description here

or

enter image description here

for the equation given, any start value for x will lead me to the second equation's root. And no one of 'em would give me the other two (roots are { -2.5, 1.18, 6.05 } and my code gives 1.18 only).

My code is something like this:

float a[] = { 1.f, -5.f, -9.f, 16.f }, etha = 0.001f;

float f(float x)
{
    return (a[0] * x * x * x) + (a[1] * x * x) + (a[2] * x) + a[3];
}

float phi(float x)
{
    return (a[3] * -1.f) / ((a[0] * x * x) + (a[1] * x) + a[2]);
}

float iterationMethod(float a, float b)
{
    float x = (a + b) / 2.f;

    while (fabs(f(x)) > etha)
    {
        x = phi(x);
    }

    return x;
}

So, calling the iterationMethod() passing ranges { -3, 0 }, { 0, 3 } and { 3, 10 } will provide 1.18 number three times along.

Where am i wrong and how should i act to get it work right?

UPD1: i do not need any third-party libraries.

UPD2: i need "Simple Iteration" algorithm exactly.

Upvotes: 0

Views: 5696

Answers (2)

Mike Seymour
Mike Seymour

Reputation: 254621

The link you posted in your comment explains why you can't find all the roots with this algorithm - it only converges to a root if |phi'(x)| < 1 around the root. That's not the case with any of the roots of your polynomial; for most starting points, the iteration will end up bouncing around the middle root, and eventually get close to it by accident; it will almost certainly never get close enough to the other roots, wherever it starts.

To find all three roots, you need a more stable algorithm such as Newton's method (which is also described in the tutorial you linked to). This is also an iterative method; you can find a root of f(x) using the iteration x -> x - f(x)/f'(x). This is still not guaranteed to converge, but the convergence condition is much more lenient. For your polynomial, it might look a bit like this:

#include <iostream>
#include <cmath>

float a[] = { 1.f, -5.f, -9.f, 16.f }, etha = 0.001f;

float f(float x)
{
    return (a[0] * x * x * x) + (a[1] * x * x) + (a[2] * x) + a[3];
}

float df(float x)
{
    return (3 * a[0] * x * x) + (2 * a[1] * x) + a[2];
}

float newtonMethod(float a, float b)
{
    float x = (a + b) / 2.f;
    while (fabs(f(x)) > etha)
    {
        x -= f(x)/df(x);
    }

    return x;
}

int main()
{
    std::cout << newtonMethod(-5,0) << '\n'; // prints -2.2341
    std::cout << newtonMethod(0,5) << '\n';  // prints  1.18367
    std::cout << newtonMethod(5,10) << '\n'; // prints  6.05043
}

There are many other algorithms for finding roots; here is a good place to start learning.

Upvotes: 3

hugomg
hugomg

Reputation: 69954

One of the more traditional root-finding algorithms is Newton's method. The iteration step involves finding the root of the first order approximation of the function

So if we have a function 'f' and are at a point x0, the linear fisrt order approximation will be

f_(x) = f'(x0)*(x - x0) + f(x0)

and the corresponding approximate root x' is

x' = phi(x0) = x0 - f(x0)/f'(x0)

(Note that you need to have the derivative function handy but it should be very easy to obtain it for polynomials)

The good thing about Newton's method is simple to implement and is often very fast. The bad thing is that sometimes it doesn't behave well: the method fails on points that have f'(x) = 0 and some inputs in some functions it can diverge (so you need to check for that and restart if needed).

Upvotes: 4

Related Questions