Root delimitation method for nonlinear equations (C++)

Article directory

  • Graphical method
    • General steps
    • Implement code
    • Case Analysis
  • dichotomy
    • General steps
    • Implement code
    • Case Analysis
  • Trial position method
    • Implement code
    • Case Analysis

Solving the roots of nonlinear equations is another very important problem in scientific research and engineering calculations. Whether it is solving differential equations in mathematical theory, optimization problems, numerical simulations, or in engineering sciences such as fluid dynamics and structural mechanics. , fields such as electromagnetic fields, and in model research in some economic disciplines and other social science disciplines, solving nonlinear equations is often encountered.

Numerical solution methods for nonlinear equations can be roughly divided into two categories: delimitation method and iterative method. Each type of method has many classic methods, each with its own characteristics. Delimitation methods mainly include graphical method, bisection method and trial position method; and Iteration methods mainly include fixed point iteration method, Newton’s method, secant method, parabola method and inverse quadratic interpolation method. In addition, there is also Brent method which combines the two.

The preprocessing required for this code is as follows:

#include <QStandardPaths>
#include <functional>
#include <QFile>
#include <math.h>
using namespace std;

Since the code in this article is written in the Qt project and uses QFile, readers can also change it to use the file stream implementation of the standard library.

Graphic method

General steps

  • Draw curves separately

    y

    =

    f

    (

    x

    )

    y=f(x)

    y=f(x) and straight line

    y

    =

    0

    y=0

    The graph of y=0.

  • observe

    y

    =

    f

    (

    x

    )

    y=f(x)

    y=f(x)sum

    y

    =

    0

    y=0

    The intersection point of y=0, the abscissa of the intersection point is the root of the equation, this value can only be read on the graph (if necessary, you can use the magnification tool to partially enlarge the graph), so the accuracy is generally not very high.

Implementation code

/*
 * Export drawing data (unary function)
 *f: function
 * a: left endpoint of the drawing interval
 * b: Right end point of the drawing interval
 * n: number of interval divisions
 * path: file saving path
 * title: header
 *
 * Return (bool):
 * true: export failed
 * false: Export successful
 */
bool plot(function<double(double)> f, double a = 0, double b = 1, unsigned n = 1000, QString path = QStandardPaths::writableLocation(QStandardPaths::DesktopLocation), const char *title = "x, f(x)")
{<!-- -->
    QFile file(path);
    if (!file.open(QIODevice::WriteOnly))
        return true;
    if (b < a)
    {<!-- -->
        doublet(a);
        a = b;
        b = t;
    }
    file.write(title);
    double step((b - a) / n + + );
    do
    {<!-- -->
        double f_a(f(a));
        if (isnan(f_a))
        {<!-- -->
            file.close();
            return true;
        }
        file.write(('\\
' + QString::number(a) + ',' + QString::number(f_a)).toUtf8());
        a + = step;
    } while (--n);
    return false;
}

Example analysis

Find function

f

(

x

)

=

e

x

+

x

?

2

f(x)=e^x + x-2

f(x)=ex + x?2 in

[

0

,

2

]

[0,2]

Root on [0,2].

Plot the function graph as follows:

It can be seen from the figure that the roots of the equation are approximately

0.45

0.45

0.45, and the exact solution obtained using MATLAB is

0.4428544010

?

0.4428544010\cdots

0.4428544010?, the error is within the allowable range, but if you want to seek higher accuracy, you need other numerical methods.

Dichotomy

The bisection method is a numerical method used to solve nonlinear equations. Its basic idea is to divide the domain of the function into two intervals, and then use linear approximation in each interval to find the roots of the equation. By continuously narrowing the search interval until a solution that meets the accuracy requirements is found.

In the bisection method, we first determine an initial interval that contains the solution to the equation. We then calculate the midpoint of the interval and calculate the value of the equation at the midpoint. If the value of the equation at the midpoint is less than zero, then the solution lies in the subinterval to the right of the midpoint; if the value of the equation at the midpoint is greater than zero, the solution lies in the subinterval to the left of the midpoint. Based on this information, we will divide the original interval into two new subintervals, One of them contains the solution, the other does not.

Repeat the above steps, reducing the subinterval containing the solution by half each time, until the subinterval containing the solution is small enough that we can determine an approximation of the solution.

General steps

Input: Real range

[

a

,

b

]

[a,b]

[a,b], function

f

(

x

)

f(x)

f(x), calculation accuracy

ε

\varepsilon

ε

Output: Approximate solution

x

x

x

  1. If

    f

    (

    a

    )

    f

    (

    b

    )

    >

    0

    f(a)f(b)>0

    f(a)f(b)>0 then

    • Return No solution
  2. Loop until

    b

    ?

    a

    < ε b-a<\varepsilon b?a<ε Until

    • m

      a

      +

      b

      2

      m\gets\dfrac{a + b}2

      m←2a + b?

    • If

      f

      (

      m

      )

      f

      (

      a

      )

      >

      0

      f(m)f(a)>0

      f(m)f(a)>0 then

      • a

        m

        a\gets m

        a←m

    • Otherwise
      • b

        m

        b\gets m

        b←m

Implementation code

/*
 * Dichotomy
 *f: function
 * a: Left endpoint of the interval
 * b: Right endpoint of the interval
 * eps: termination condition
 *
 * Return (double): approximate solution
 */
double dichotomy(const function<double(double)> & amp;f, double a, double b, double eps)
{<!-- -->
    if (a > b)
    {<!-- -->
        double t = a;
        a = b;
        b = t;
    }
    if (f(a) > 0)
    {<!-- -->
        if (f(b) > 0)
            return NAN;
        while (b - a > eps)
        {<!-- -->
            double m((a + b) / 2);
            if (isnan(f(m)))
                return NAN;
            if (f(m) > 0)
                if (a == m)
                    break;
                else
                    a = m;
            else
                if (b == m)
                    break;
                else
                    b = m;
        }
    }
    else
    {<!-- -->
        if (f(b) < 0)
            return NAN;
        while (b - a > eps)
        {<!-- -->
            double m((a + b) / 2);
            if (isnan(f(m)))
                return NAN;
            if (f(m) > 0)
                if (b == m)
                    break;
                else
                    b = m;
            else
                if (a == m)
                    break;
                else
                    a = m;
        }
    }
    return (a + b) / 2;
}

Example analysis

Find functions using the bisection method

f

(

x

)

=

sin

?

x

?

x

?

x

2

+

0.5

f(x)=\sin x-x-x^2 + 0.5

f(x)=sinx?x?x2 + 0.5 is in the interval

[

0

,

1

]

[0,1]

Zero point within [0,1].

Substitute into the program to get the changes in the left and right endpoints of the interval as shown in the figure:

Test position method

The trial position method is a modification of the bisection method. In the bisection method, if the midpoint is changed to the point

(

a

,

f

(

a

)

)

(a,f(a))

(a,f(a)) and

(

b

,

f

(

b

)

)

(b,f(b))

(b,f(b)) Same as

x

x

Linear interpolation of the x-axis, that is

m

=

b

?

b

?

a

f

(

b

)

?

f

(

a

)

f

(

b

)

m=b-\dfrac{b-a}{f(b)-f(a)}f(b)

m=b?f(b)?f(a)b?a?f(b)
This gives us the trial position method.

Implementation code

/*
 * Dichotomy
 *f: function
 * a: Left endpoint of the interval
 * b: Right endpoint of the interval
 * eps: termination condition
 *
 * Return (double): approximate solution
 */
double dichotomy(const function<double(double)> & amp;f, double a, double b, double eps)
{<!-- -->
    if (a > b)
    {<!-- -->
        double t = a;
        a = b;
        b = t;
    }
    if (f(a) > 0)
    {<!-- -->
        if (f(b) > 0)
            return NAN;
        while (b - a > eps)
        {<!-- -->
            double m(b - (b - a) * f(b) / (f(b) - f(a)));
            if (isnan(f(m)))
                return NAN;
            if (f(m) > 0)
                if (a == m)
                    break;
                else
                    a = m;
            else
                if (b == m)
                    break;
                else
                    b = m;
        }
    }
    else
    {<!-- -->
        if (f(b) < 0)
            return NAN;
        while (b - a > eps)
        {<!-- -->
            double m(b - (b - a) * f(b) / (f(b) - f(a)));
            if (isnan(f(m)))
                return NAN;
            if (f(m) > 0)
                if (b == m)
                    break;
                else
                    b = m;
            else
                if (a == m)
                    break;
                else
                    a = m;
        }
    }
    return b - (b - a) * f(b) / (f(b) - f(a));
}

Example analysis

Use the trial position method to solve the above equation.

The changes in the left and right endpoints of the available interval are as shown in the figure: