В Python собираю функцию решения уравнения методом дихотомии

Question

Ищу решение уравнения. Функция состоит из двух частей:

отделение корня;
уточнение корня.

На первом этапе формируем сетку и строим словарь: значение Х: значение Y. Если от одного значения к другому Y меняет знак на противоположный, мы знаем: на этом участке между соседними Х имеется корень.

На втором этапе мы его уточняем: делим отрезок пополам. Корень может быть серединой отрезка либо же находиться левее/правее от него. Производим действия, и когда путем приближения к корню преодолеваем границу погрешности (новое значение становится меньше погрешности), считаем, что корень найден. Выполнение кода показывает только один корень (x = 92.005), хотя корней у функции много больше (ее график многократно пересекает ось Х). Первый из корней примерно равен, по прикидкам, 0,065. Сейчас выход такой:

Корень по методу дихотомии находится в точке x = 92.005

В чем оплошность?

import numpy as np

def f(x):
    return 1.2-np.log(x)-4*np.cos(2*x)

def dichotomy (a,b,n, eps): # отрезок от a до b делим на n частей, погрешность eps
    """
    Функция отделения и уточнения корня
    """
    assert a!=0,  'a равно 0'
    assert b!=0, 'b равно 0'
    # сначала отделим корни
    h=(b-a)/n
    setka={}
    for i in np.arange(a, b, 0.1):
        x=a+i*h
        i=i+0.1
        setka.update({x:f(x)})
    
    # далее уточним корни
    root = None
    for value in setka.values():
        while ( abs(f(b)-f(a)) )>eps:     # пока отрезок больше заданной погрешности, выполняем нижестоящие операции:
            mid = (a+b)/2                   # получаем середину отрезка
            if f(mid) == 0 or f(mid)<eps:    # если функция в середине отрезка равну нулю или меньше погрешности:
                root = mid                  # корень равень серединному значению
                break
            elif (f(mid) * f(a)) < 0:       # иначе если произведение функции в середине отрезка на функцию в т. а <0
                b = mid                     # серединой становится точка b
            else:
                a = mid                     #в другом случае - точка а
    if root is None:
        print('Корень не найден')
    else:
        print(f'Корень по методу дихотомии находится в точке x = {root}')

res = dichotomy(0.01,184,100, 0.001)
res

Answer 1

Пришлось сильно переписать код. Вообще не я не очень знаю постановку этой задачи, возможно, первоначально отрезки должны подбираться другим способом, но мой способ по крайней мере работает.

import numpy as np

def f(x):
    return 1.2-np.log(x)-4*np.cos(2*x)

def dichotomy (a,b,n, eps): # отрезок от a до b делим на n частей, погрешность eps
    """
    Функция отделения и уточнения корня
    """
    assert a!=0,  'a равно 0'
    assert b!=0, 'b равно 0'
    # сначала отделим корни
    setka=np.linspace(a, b, n)
    # далее уточним корни
    for x,y in zip(setka, setka[1:]):
        if f(x) * f(y) > 0: # если на отрезке нет корня, смотрим следующий
            continue
        root = None
        while ( abs(f(y)-f(x)) )>eps:     # пока отрезок больше заданной погрешности, выполняем нижестоящие операции:
            mid = (y+x)/2                   # получаем середину отрезка
            if f(mid) == 0 or f(mid)<eps:    # если функция в середине отрезка равну нулю или меньше погрешности:
                root = mid                  # корень равень серединному значению
                break
            elif (f(mid) * f(x)) < 0:       # иначе если произведение функции в середине отрезка на функцию в т. а <0
                y = mid                     # серединой становится точка b
            else:
                x = mid                     #в другом случае - точка а
        if root:
            yield root

res = list(dichotomy(0.01,10,1000, 0.00001))
print('Корни по методу дихотомии находятся в точках:')
print(', '.join(map(lambda x: f'{x:.4f}', res)))

Вывод:

Корни по методу дихотомии находятся в точках:
0.0650, 0.5525, 2.3975, 3.9450, 5.4375, 7.1650, 8.5250

Судя по графику, действительно на участке [0,10] должно быть 7 корней примерно в этих местах.

Пояснения:

вы делали сетку как-то странно. Её можно сделать гораздо проще.
вы эту сетку потом вообще никак не использовали, а искали корень по всему промежутку [a,b]
Чтобы получить корень в окрестности числа 0,065 как вы думаете, на сколько частей нужно разделить отрезок [0, 184], если на нём очень-очень много корней? А если разделить отрезок [0, 184] всего на 100 частей, то как вы думаете - найдёте ли вы корень 0.065 на первом отрезке или скорее что-то посередине между 0 и 1.84, т.е. примерно где-то возле числа 0.92, если там тоже есть корни?
Вы вообще не делали проверку на то, что на отрезке есть корень, я добавил такую проверку

Answer 2

Добавлю свои пять копеек. Сравнил с исходными решениями, вроде сошлось.

Разделил функции на основную findroot() и бисекции bisect(). Бисекции в условиях задачи упростил, принимая на последнем шаге xm = 0.5 * (x1 + x2).

Ну и в findroot() сразу подставляется функция для поиска корней.

Число разбиений поменял на начальную точность.

import numpy as np


def f(x):
    return 1.2-np.log(x)-4*np.cos(2*x)


def bisect(func, xs: float, xe: float, eps: float):
    if func(xs) * func(xe) < 0.0:
        while abs(xe - xs) > eps:
            xm = 0.5 * (xs + xe)
            xs, xe = (xs, xm) if func(xs) * func(xm) < 0.0 else (xm, xe)
        return 0.5 * (xs + xe)
    return None


def findroot(func, a, b, eps_init, eps_bisect):
    res = []
    setka = np.arange(a, b, eps_init)
    size = len(setka)
    points = list(zip(setka, func(setka)))
    for idx, v in enumerate(points):
        next_idx = idx + 1 if idx + 1 < size else idx
        if v[1] * points[next_idx][1] < 0.0:
            res.append(bisect(f, v[0], points[next_idx][0], eps_bisect))
    return res


res = findroot(f, 0.0001, 10, 0.001, 0.0000001)
print(('{:.5} ' * len(res)).format(*res))

Результат

0.062754 0.55311 2.397 3.9487 5.436 7.1653 8.5205

Сравнение с Wolfram Alpha

БЛОГ НА HUSL

В Python собираю функцию решения уравнения методом дихотомии

Ответы (2 шт):