Как сделать метод деления пополам?

Question

[Михаил Егоров]

Нужно решить уравнение методом деления пополам. Нарисовал график в Desmos - функция пересекает 0 в 3 местах. Моя программа ищет только 1 корень. Вопрос: как найти еще и 2 других?

Уравнение:
x^3((e^x+e^-x)/2)+pi-9pix=0 [-10.0; 10.0]

Код:

double f2(double a, double b) {
    double t = (b + a) / 2;
    //printf("%.10lf\n", t);
    while (abs(b - a) > 0.0000001) {
        if (riv2(t) > t)
        {
            b = t;
            printf("%.10lf\n", t);
        }
        else
        {
            a = t;
            //printf("%.10lf\n", riv2(t));
        }
        t = (b + a) / 2;
        //printf("%.10lf\n", t);
    }
    
    return t;
}

Функция riv2() подставляет t вместо x, считает и возвращает результат

Comments

[Akina]

Метод половинного деления предназначен только для случаев, когда в исследуемом диапазоне имеется строго один корень. Если их несколько, диапазон следует раздробить (например, сканированием с неким шагом), а в каждом выполнять поиск независимо.

Общего алгоритма сканирования/дробления нет...

Теоретически МПД даже может не обнаружить факта наличия нескольких корней - например, если все корни, кроме одного, находятся в первой половине диапазона.

[Михаил Егоров]

Akina, Спасибо за ответ

[Rsa97]

Можно предварительно пройти по отрезку с маленьким шагом, выделяя интервалы, где значение функции меняет знак, но этот метод тоже не гарантирует, что внутри даже маленького интервала функция не сменит знак дважды, трижды и т.п.

[mayton2019]

Akina, надо это отписать в ответы.

Answer 1

[Игорь Кузьмич]

1) прошагать с мелким шагом найти где функция меняет знак
2) в каждом диапазоне где функция меняет знак найти корень
2.1) если функция при изменении знака очень разко меняет значение (оценка производной слишком большая), то это не пересечения нуля, а разрыв.

такой алгоритм не гарантирует корректное решение.

#define _USE_MATH_DEFINES
#include <iostream>
#include <vector>
#include <optional>
#include <cmath>
#include <cassert>
using namespace std;

double f(double x)
{
	return x * x * x * (cosh(x)) + M_PI - 9 * M_PI * x;// (e^x+e^-x)/2 это "гиперболический косинус" aka "чосинус"
}

bool has_different_signs(double a, double b)
{
	return (a > 0) != (b > 0);
}

vector<pair<double, double>> find_sign_shanges(double from, double to, double step, double function(double))
{
	assert(from < to);
	assert(step > 0);

	vector<pair<double, double>> result;

	double last_val = function(from);
	double last_x = from;
	from += step;
	for (; from < to; last_x = from, from += step)
	{
		double current_val = function(from);
		if (has_different_signs(current_val, last_val))
		{
			result.push_back({last_x, from});
		}
		last_val = current_val;
		last_x = from;
	}
	return result;
}

optional<double> find_root_binary(double from, double to, double function(double), double precision, double slope_limit = +INFINITY)
{
	double y1 = function(from);
	double y2 = function(to);
	assert(from <= to); // не уверен что это нужно
	assert(has_different_signs(y1, y2));
	assert(slope_limit > 0);
	assert(precision > 0);

	double slope = (y2 - y1) / (to - from);
	if (abs(slope) >= slope_limit)
	{
		return nullopt; // если угол наклона граффика очень боьшой, кторее всего знак меняется не через 0, а через разрыв
	}

	double middle_x = (from + to) / 2;
	double middle_y = function(middle_x);
	if (abs(middle_y) < precision)
	{
		return middle_x;
	}
	if (has_different_signs(middle_y, y1))
	{
		return find_root_binary(from, middle_x, function, slope_limit);
	}
	else
	{
		return (find_root_binary(middle_x, to, function, slope_limit));
	}
}

vector<double> find_roots(double from, double to, double function(double), double precision, double step = 0.1)
{
	assert(precision > 0);
	auto sign_shanges = find_sign_shanges(from, to, step, function);
	vector<double> result;
	for (auto [from_mini, to_mini] : sign_shanges) // нужен c++ 17 или более поздняя версия вроде бы как для такого синтаксиса
	{
		auto root = find_root_binary(from_mini, to_mini, function, precision);
		if (root)
		{
			result.push_back(*root);
		}
	}
	return result;
}

int main()
{
	auto roots = find_roots(-10, +10, f, 0.0000001);
	for (double root : roots)
	{
		cout << root << endl;
	}
	return 0;
}