Нелинейные уравнения одной переменной

Пусть имеем некоторую функцию одной переменной \(f(x)\). Нас интересуют нули этой функции, т.е. задача решения уравнения

\[ f(x)=0 \]

Для простоты пока работаем с функцией одной переменной.

Теоремы о локализации корней алгебраических уравнений.

Теорема (о промежуточном значении)

Если непрерывная функиия \(f(x)\) принимает на концах отрезка значения разного знака, т.е. \(f(a) f(b)<0\), то внутри отрезка имеется, по меньшей мере, один корень уравнения

\[ f(x) = 0 \]

Для более конкретных рекомендаций требуется дополнительная информация об изучаемом классе функций \(f(x)\).

Рассмотрим более подробно алгебраические уравнения вида

\[ f(z)=a_{0} z^{n}+a_{1} z^{n-1}+\cdots+a_{n-1} z+a_{n}=0, a_{0} \neq 0 \]

c действительными коэффициентами \(\left\{a_{i}\right\}_{i=0}^{n}, n>1\).

Приведём без доказательства несколько утверждений, оказывающихся полезными для локализации действительных корней уравнения.

Теорема о локализации корней алгебраического многочлена

Все корни действительного алгебраического уравнения расположены на комплексной плоскости в кольце

\[ \frac{\left|a_{n}\right|}{\left|a_{n}\right|+B} \leq|z| \leq 1+\frac{A}{\left|a_{0}\right|}, \text { где } A=\max _{1 \leq i \leq n}\left\{\left|a_{i}\right|\right\}, B=\max _{0 \leq i \leq n-1}\left\{\left|a_{i}\right|\right\} \text {. } \]

Теорема Декарта

Число положительных действительных корней уравнения либо равно числу перемен знака в последовательности \(\left\{\mathbf{a}_{\mathrm{i}}\right\}_{\mathrm{i}=0}^{\mathrm{n}}\) коэффичиентов, причём коэффициенты, равные нулю, не учитываются, либо меньше этого числа на чётное число.

Пусть \(c\) - действительное число и \(f(c) \neq 0\). Через выражение \(S(c+0)\) обозначим число перемен знаков в последовательности функции и её производных \(\left\{f_{x}^{(i)}(c)\right\}_{i=0}^{n}\) в точке \(c\), где \(f_{x}^{(i)}(c+0)=0\) присваивается знак в точке, сколь угодно близкой к \(c\) справа, а \(f_{x}^{(i)}(c-0)=0\) - слева.

Теорема Бюдана-Фурье

Если \(a\) и \(b\,(a<b)\) не являются корнями алгебраического уравнения указанного выше вида, то число действительных корней этого уравнения с учётом их кратности на \([a ; b]\) равняется разности между числом перемен знаков в последовательностях \(S(a+0)\) и \(S(b-0)\) или меньше этого числа на чётное число

\[ \Delta=S(a+0)-S(b-0) \]

Составим последовательность \(\left(f_{i}\right)_{i=0}^{n}:\) \(f_{0}=f(x)\), \(f_{1}=f_{x}^{\prime}(x)\), \(f_{2}=-\operatorname{Res}\left\{f_{0} / f_{1}\right\}\), \(\ldots\), \(f_{n}=-\operatorname{Res} \left\{f_{n-2} / f_{n-1}\right\}\), где \(\operatorname{Res} \left\{f_{p}(x) / f_{q}(x)\right\}-\) остаток от деления многочлена \(f_{p}(x)\) на многочлен \(f_{q}(x)\), и назовём её рядом Штурма.

Теорема Штурма

Число действительных корней уравнения, если среди них нет кратных, равняется разности между числом перемены знаков в ряде Штурма при \(x=a\) и \(x=b\).

Пример

Локализовать действительные положительные корни многочлена

\[ f(x)=-32 x^{3}+20 x^{2}+11 x+3 . \]

По следствию из основной теоремы алгебры, все действительные корни лежат на пересечении кольца в комплексной плоскости с внутренним радиусом \(\alpha\), внешним радиусом \(\beta\) и положительной частью действительной оси, т.е. в интервале \([\alpha, \beta]\), где

\[ \alpha=\frac{\left|a_{3}\right|}{\left|a_{3}\right|+B}=\frac{3}{35}, \beta=1+\frac{|A|}{\left|a_{0}\right|}=\frac{13}{8} \]

т.к. \(A=\max \{20,11,3\}=20, \quad B=\max \{32,20,11\}=32, \quad\left|a_{0}\right|=32, \quad\left|a_{3}\right|=3\).

Число корней можно определить одним из следующих способов:

a) По теореме Декарта

По теореме Декарта число перемен знаков в последовательности коэффициентов \(\{-32,20,11,3\}\) многочлена \(f(x)\) равно 1 , т.е. \(f(x)\) имеет один положительный корень на отрезке \(\left[\frac{3}{35}, \frac{13}{8}\right]\);

б) По теореме Бюдана-Фурье

По теореме Бюдана-Фурье последовательность функции и ее производных

\[\begin{split} \begin{gathered} f(x)=-32 x^{3}+20 x^{2}+11 x+3, f(\beta) \approx-7.375<0, f(\alpha) \approx 4.07>0, \\ f^{\prime}(x)=-96 x^{2}+40 x+11, f^{\prime}(\beta) \approx-177.5<0, f^{\prime}(\alpha) \approx 13.72>0, \\ f^{\prime \prime}(x)=-192 x+40, f^{\prime \prime}(\beta) \approx-272<0, f^{\prime \prime}(\alpha) \approx 23.54>0, \\ f^{\prime \prime \prime}(x)=-192=f^{\prime \prime \prime}(\beta)=f^{\prime \prime \prime}(\alpha)<0, \\ \Delta=S(\alpha+0)-S(\beta-0)=1-0=1, \end{gathered} \end{split}\]

т.е. \(f(x)\) имеет один положительный корень на отрезке \(\left[\frac{3}{35}, \frac{13}{8}\right]\);

в) По теореме Штурма

По теореме Штурма составим ряд Штурма:

\[\begin{split} \begin{aligned} & f_{0}(x)=f(x)=-32 x^{3}+20 x^{2}+11 x+3, \\ & f_{1}(\beta) \approx-7.375<0, f_{1}(\alpha) \approx 4.07>0, \\ & f_{1}(x)=f^{\prime}(x)=-96 x^{2}+40 x+11, \\ & f_{1}(\beta) \approx-177.5<0, f_{1}(\alpha) \approx 13.72>0, \\ & f_{2}(x)=-\operatorname{Re} s\left\{f_{0}(x) / f_{1}(x)\right\}=-\frac{91}{9} x-\frac{271}{72}, \\ & f_{2}(\beta) \approx-20.19<0, f_{2}(\alpha) \approx-4.63<0, \\ & f_{3}(x)=-\operatorname{Re} s\left\{f_{1}(x) / f_{2}(x)\right\}=-\frac{284751}{16562}=f_{3}(\beta)=f_{3}(\alpha)<0, \end{aligned} \end{split}\]

т.е. \(f(x)\) имеет один положительный корень на отрезке \(\left[\frac{3}{35}, \frac{13}{8}\right]\).

Метод дихотомии (деления пополам)

Пусть функция \(f\) - непрерывная.

Пусть имеем две точки \(x_{\text{left}} < x_{\text{right}}\) такие что функция на них имеет разный знак:

\[ f(x_{\text{left}} )\cdot f(x_{\text{right}}) < 0 \]

Тогда ясно, что где-то в интервале \((x_{\text{left}}, x_{\text{right}})\) будет корень искомого уравнения (по теореме о промежуточных значениях непрерывной функции).

Идея метода заключается в последовательном сужении отрезка до тех пор, пока не достигнем нужной точности корня. Алгоритм:

Алгоритм метода дихотомии

1. Фиксируем малый \(\varepsilon\)

2. Пусть \(x_0=\frac{x_{\text{left}}+x_{\text{right}}}{2}\)

   • Если \(f(x_0) \cdot f(x_{\text{left}}) < 0\), то они имеют разный знак и принимаем \(x_{\text{right}} = x_0\)
     
     

   • Иначе, если \(f(x_0) \cdot f(x_{\text{left}}) > 0\), они имеют одинаковый знак и принимаем \(x_{\text{left}} = x_0\)
     
     

3. Если \(|x_{\text{left}} - x_{\text{right}}| < \varepsilon\), то локализировали корень с точностью \(\varepsilon\) и останавливаемся. Иначе, обратно на 2 шаг.

Метод крайне простой и эффективный. Однако, он не применим для поиска корней чётной кратности (в которой корень - локальный минимум или максимум). Также не обобщается на системы уравнений и на уравнения с комплекснозначной функцией.

Метод Ньютона

Пусть теперь \(f\) - дифференциируема.

Пусть имеем некоторой приближение корня \(x_n\). Смотрим на касательную к \(f(x)\) в этой точке и ищем её пересечение с \(Ox\). Это пересечение есть наше следующее приближение корня. Повторяем, пока не достигнем нужной точности.

Формально,

Шаг метода Ньютона

\[ x_{n+1}=x_n-\frac{f\left(x_n\right)}{f'\left(x_n\right)} \]

Этот метод работает быстрее метода дихотомии и с помощью него можно найти любые корни даже для комплекснозначных функций.

Тем не менее, этот метод имеет существенный недостаток - нестабильная сходимость (особенно при численной аппроксимации производной). Например, на каком-то этапе производная может близко приблизиться к нулю, и мы улетим далеко-далеко от корня. Также возможны замкнутые циклы - ноль может не достигнуться.

Пример с комплекснозначной функцией

Проиллюстрируем метод Ньютона для вычисления корней комплекснозначной функции:

\[ f(z) = z^4 - 1 \]

На итоговой картинке одинаковым цветом обозначены «области притяжения корней» - стартуя алгоритм в одной из этих областей, мы будем призодить к одному и тому же корню. Из этого можно сделать вывод, что метод Ньютона нелокален.

Здесь код на питоне.

Show code cell source Hide code cell source

from scipy import optimize # В этом модуле есть все рассказанные методы оптимизации (нахождения корней или минимумов)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

# A list of colors to distinguish the roots.
colors = ['xkcd:eggshell', 'mediumseagreen', 'c', 'xkcd:sky blue']

TOL = 1.e-8

def newton(z0, f, fprime, MAX_IT=1000):
    """The Newton-Raphson method applied to f(z).

    Returns the root found, starting with an initial guess, z0, or False
    if no convergence to tolerance TOL was reached within MAX_IT iterations.

    """

    z = z0
    for i in range(MAX_IT):

        # Вместо всей этой секции z = optimize.newton(f, z0, fprime)

        dz = f(z)/fprime(z)
        if abs(dz) < TOL:
            return z
        z -= dz
    return False

def plot_newton_fractal(f, fprime, n=200, domain=(-1, 1, -1, 1)):
    """Plot a Newton Fractal by finding the roots of f(z).

    The domain used for the fractal image is the region of the complex plane
    (xmin, xmax, ymin, ymax) where z = x + iy, discretized into n values along
    each axis.

    """

    roots = []
    m = np.zeros((n, n))

    def get_root_index(roots, r):
        """Get the index of r in the list roots.

        If r is not in roots, append it to the list.

        """

        try:
            return np.where(np.isclose(roots, r, atol=TOL))[0][0]
        except IndexError:
            roots.append(r)
            return len(roots) - 1

    xmin, xmax, ymin, ymax = domain
    for ix, x in enumerate(np.linspace(xmin, xmax, n)):
        for iy, y in enumerate(np.linspace(ymin, ymax, n)):
            z0 = x + y*1j
            r = newton(z0, f, fprime)
            if r is not False:
                ir = get_root_index(roots, r)
                m[iy, ix] = ir
    nroots = len(roots)
    if nroots > len(colors):
        # Use a "continuous" colormap if there are too many roots.
        cmap = 'hsv'
    else:
        # Use a list of colors for the colormap: one for each root.
        cmap = ListedColormap(colors[:nroots])
    plt.figure(figsize=(4,4))
    plt.imshow(m, cmap=cmap, origin='lower')
    plt.axis('off')
    plt.show()

    return roots

f = lambda z: z**4 - 1
fprime = lambda z: 4*z**3 

plot_newton_fractal(f, fprime, n = 500)

[(9.075876832133421e-10+1.00000000209645j),
 (0.9999999997118054-8.529533586119149e-10j),
 (9.595113291316482e-17-1.0000000000000002j),
 (-1.0000000000001743-2.3209727148473305e-13j)]

Метод секущих

Пусть \(f\) - непрерывна.

Будем искать нуль функции \(f(x)\). Выберем две начальные точки \(B_0\left(x_0 , f(x_0)\right)\) и \(B_1\left(x_1 , f(x_1)\right)\) и проведем через них прямую. Она пересечет ось абсцисс в точке \(\left(x_2 , 0\right)\).

Теперь найдем значение функции с абсциссой \(x_2\). Зададим точку \(B_2 = (x_2, f(x_2))\), лежащую на графике.

Теперь вместо точек \(B_0\) и \(B_1\) мы возьмём точку \(B_2\) и точку \(B_1\). Теперь с этими двумя точками проделаем ту же операцию и так далее, то есть будем получать две точки \(B_{n+1}\) и \(B_n\) и повторять операцию с ними. Отрезок, соединяющий последние две точки, пересекает ось абсцисс в точке, значение абсциссы которой можно приближённо считать корнем.

Эти действия нужно повторять до тех пор, пока не получим значение корня с нужным приближением.

Формульно,

Шаг метода секущих

\[ x_{n+1}=x_{n-1}-\frac{f\left(x_{n-1}\right) \cdot\left(x_n-x_{n-1}\right)}{f\left(x_n\right)-f\left(x_{n-1}\right)} \]

Скорость у него не такая большая, как у метода Ньютона, но он работает на бОльшем классе функций.

Метод простой итерации (релаксации)

Из математического анализа нам известно, что если некоторая функция \(\Phi(x)\) является сжимающим отображением:

\[ | \Phi(x) - \Phi(y)| \leqslant \alpha \cdot |x - y|, \quad \alpha < 1 \]

то она имеет стационарную точку.

В контексте нашей задачи, пусть свели исследуемое уравнение к виду:

\[ f(x)=0 \quad \iff \quad \Phi(x)=x, \quad \text{где $\Phi$ - cжимающее} \]

Тогда корень будет стационарной точкой \(F\) и его можно найти методом простой итерации:

Шаг метода простой итерации (МПИ)

\[ x_{n+1} = \Phi\left(x_n\right) \]

Если \(f\) - дифференциируема, то проверить сжимаемость \(\Phi\) легко (достаточное условие):

\[ |\Phi(x) - \Phi(y)| = \frac{|\Phi(x) - \Phi(y)|}{|x-y|} \cdot |x-y| \leq \max |\Phi'| \cdot |x-y| \]

\[ \max |\Phi'| < 1 \]

Пример и метод релаксации

Сжимаемость в данном случае критично важна. Например, следующий пример простой итерации не сходится никогда (кроме случае, если мы уже начали в корне):

\[ x^2=a \iff x=\frac{a}{x} \]

\[ \Phi(x)=\frac{a}{x} \]

Как видим, может просто не повезти с функцией \(\Phi\) - она получилась не сжимающей.

Для того, чтобы избежать этого, метод простой итерации некоторым образом обобщают (метод релаксации), при котором мы задаём эквивалентную сжимающую функцию \(\Phi(x)\) в следующем виде:

\[ \Phi(x)=x-\lambda(x) f(x), \quad \lambda(x) \neq 0 \]

и пользуемся прежней итерационной формулой

\[ x_{n+1} = \Phi\left(x_n\right) \]

В контексте примера, \(\lambda\) можно выбрать \(\pm \frac{1}{2a}\) - это создаст «области сжимаемости» около искомых корней.

Примечание. Метод Ньютона является частным случаем метода релаксации (выбором подходящего \(\lambda(x)=\cfrac{1}{f'(x)}\)).

Порядок сходимости метода простых итераций.

Порядок сходимости метода простых итераций

Итерационный процесс \(x_{n+1} = \Phi\left(x_n\right)\) имеет p-порядок сходимости, если

\[ |x_{n+1} - x^*| \leq C |x_{n} - x^*|^p \]

Теорема о порядке сходимости метода простых итераций

Итерационный процесс \(x_{n+1} = \Phi\left(x_n\right)\) имеет p-порядок сходимости, если

\[ \Phi'\left(x^*\right) = \Phi''\left(x^*\right) = ... = \Phi^{(p-1)}\left(x^*\right) = 0, а \ \Phi^{(p)}\left(x^*\right) \neq 0 . \]

Как найти все корни?

Все рассматриваемые ранее методы позволяли найти лишь один корень. А что если нам нужны все? Нетрудно догадаться, что делать.

Пусть мы нашли какое-то количество корней: \(x_1, ..., x_k\) нашей \(f(x)\). Строим функцию \(g(x) = \frac{f(x)}{(x-x_1)...(x-x_k)}\) и находим какой-либо её корень!

Вот только после каждого деления мы повышаем ошибку округления. Таким образом, разумно с помощью данного подхода находить приближение к какому-либо новому корню \(x_{k+1}^*\) и используем его как начальное приближение при поиске корня уже у самой \(f(x)\).

Пример. Нахождение всех корней многочлена

Найдём корни полинома \(x^4-11 x^3+9 x^2+11 x-10\).

Какую (приблизительно) долю площади на комплексной плоскости от круга \(|z|<10\) занимает область притяжения наибольшего корня по модулю?

Решим перебором точек - запустим метод Ньютона из каждой точки на равномерной сетке в квадрате со стороной 10. Код на питоне.

Show code cell source Hide code cell source

# A list of colors to distinguish the roots.
colors = ['xkcd:eggshell', 'mediumseagreen', 'c', 'xkcd:sky blue']

TOL = 1.e-8

def newton(z0, f, fprime, MAX_IT=1000):
    """The Newton-Raphson method applied to f(z).

    Returns the root found, starting with an initial guess, z0, or False
    if no convergence to tolerance TOL was reached within MAX_IT iterations.

    """

    z = z0
    for i in range(MAX_IT):

        # Вместо всей этой секции z = optimize.newton(f, z0, fprime)

        dz = f(z)/fprime(z)
        if abs(dz) < TOL:
            return z
        z -= dz
    return False

def plot_newton_fractal(f, fprime, n=200, domain=(-1, 1, -1, 1)):
    """Plot a Newton Fractal by finding the roots of f(z).

    The domain used for the fractal image is the region of the complex plane
    (xmin, xmax, ymin, ymax) where z = x + iy, discretized into n values along
    each axis.

    """

    roots = []
    m = np.zeros((n, n))

    def get_root_index(roots, r):
        """Get the index of r in the list roots.

        If r is not in roots, append it to the list.

        """

        try:
            return np.where(np.isclose(roots, r, atol=TOL))[0][0]
        except IndexError:
            roots.append(r)
            return len(roots) - 1

    xmin, xmax, ymin, ymax = domain
    for ix, x in enumerate(np.linspace(xmin, xmax, n)):
        for iy, y in enumerate(np.linspace(ymin, ymax, n)):
            z0 = x + y*1j
            r = newton(z0, f, fprime)
            if r is not False:
                ir = get_root_index(roots, r)
                m[iy, ix] = ir


    count_r_a = [0, 0]
    max_ = -1
    i_max = -1

    for i, root in enumerate(roots):
        abs_ = root.real**2 + root.imag**2

        if abs_ > max_:
            max_ = abs_
            i_max = i

    for ix, x in enumerate(np.linspace(xmin, xmax, n)): # Считаем долю площади наиб. корня
        for iy, y in enumerate(np.linspace(ymin, ymax, n)):
            if x**2 + y**2 < 10**2:
                count_r_a[1] += 1

                if (m[iy, ix] == i_max):
                    count_r_a[0] += 1

    nroots = len(roots)
    if nroots > len(colors):
        # Use a "continuous" colormap if there are too many roots.
        cmap = 'hsv'
    else:
        # Use a list of colors for the colormap: one for each root.
        cmap = ListedColormap(colors[:nroots])

    plt.figure(figsize=(4,4))
    plt.imshow(m, extent = [xmin, xmax, ymin, ymax], cmap=cmap)

    for i, root in enumerate(roots):
        plt.scatter(root.real, root.imag)

    plt.show()

    return roots, float(count_r_a[0]/count_r_a[1])*np.pi*100

f = lambda z: z**4 - 11*z**3 + 9*z**2 + 11*z - 10
fprime = lambda z: 4*z**3 - 33*z**2 + 18*z + 11

domain=(-10, 10, -10, 10)
n = 400

roots, m = plot_newton_fractal(f, fprime, n=n, domain=domain)
print(roots)
print(f"S = {m}")

[(9.999999999999996+1.0423244797614476e-15j), (1.000000014266004-8.680596700749564e-09j), (-1.0000000000250509-9.374439725630672e-12j)]
S = 8.792839736163463

# Проверим, что нашли все корни с помощью аналитических вычислений
import sympy as smp # Библиотека для символьных вычислений

x = smp.symbols("x")
f = x**4 - 11*x**3 + 9*x**2 + 11*x - 10
display(f)

\[\displaystyle x^{4} - 11 x^{3} + 9 x^{2} + 11 x - 10\]

f1 = (f/(x-10)).simplify()
display(f1)

\[\displaystyle x^{3} - x^{2} - x + 1\]

f2 = (f1/(x-1)).simplify()
display(f2) 

\[\displaystyle x^{2} - 1\]

f3 = (f2/(x-1)).simplify() # Видим, что это кратный корень
display(f3)

\[\displaystyle x + 1\]

f4 = (f3/(x+1)).simplify()
display(f4)

\[\displaystyle 1\]