Solving non linear equations 1.1

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <cmath>
#include <vector>
#include <stdexcept>

// Для возведения в квадрат
#define sqr(a) ((a)*(a))



// Отрезок [a;b] на котором находится один корень уравнения
struct Section {
    double leftBorder; // a
    double rightBorder;// b
};



// Точность эпсилон
const double EPS = 0.001;



// Прототипы функций для уравнений
double f1(double x); // [0 ; 1] root is 0.1
double f2(double x); // just y = x
double f3(double x); // [1.5 ; 3.14] root is 2
double f4(double x); // just y = sin(x)
double f5(double x); // just x*x

// Тип указателя на функции выше
typedef double (*pointFunc)(double x);



// Прототипы функций для вычисления производных
double derivative(pointFunc f, double x);
double doubleDerivative(pointFunc f, double x);

// Прототип функции sign() - возвращает знак (+1/-1)
int sign(double x);



// Прототипы функций для решения уравнения
double halfDivisionMethod(pointFunc f, const Section & section);
double chordMethod(pointFunc f, const Section & section);
double newtonMethod(pointFunc f, const Section & section);
double interationMethod(pointFunc f, const Section & section);

// Тип указателя на функции выше
typedef double(*pointMethod)(pointFunc f, const Section & section);



int main() {
    // Здесь будут хранится отрезки, содержащие один корень уравнения
    std::vector<Section> sections;

    // Здесь будет храниться найденные корни уравнения
    std::vector<double> roots;



    //############################################################################
    //# Ввод исходных данных и разбиение функции на интервалы (отделение корней) #
    //############################################################################

    const pointFunc f = f1; // Здесь задаём уравнение через указатель на функцию
    int n = 0;
    double leftBorder = 0.0;
    double rightBorder = 0.0;

    // Вводим область определения, на которой хотим искать корни
    // Также вводим кол-во отрезков, необходимое для поиска всех корней
    initInput:
    printf("Input left border: ");
    scanf("%lf", &leftBorder);
    printf("Input right border: ");
    scanf("%lf", &rightBorder);
    printf("Input amount of sections: ");
    scanf("%d", &n);
    if (rightBorder <= leftBorder) {
        puts("Error: incorrect input of section borders! Try again!\n");
        goto initInput;
    }
    if (n <= 0) {
        puts("Error: incorrect input amount of sections! Try again!\n");
        goto initInput;
    }

    // Запоминаем те отрезки [a; b], где есть корень
    // Таковыми будем считать все отрезки, знак функции на концах которых различается
    // Т.е. f(a) * f(b) = -1
    double x = leftBorder;
    double h = (rightBorder - leftBorder) / double(n);
    while (x < rightBorder)
    {
        if (f(x) * f(x + h) <= 0)
            sections.push_back({ x, x + h });
        x += h;
    }
    if (sections.empty()) { // Если таких отрезков нет
        printf("Error: something went wrong trying to get sections.");
        printf("Ensure equation have roots or increase N.");
        getchar();
        getchar();
        exit(1);
    }
    puts("");
    //sections.push_back({ -1.0, 1.0 });



    //############################################
    //# Выбор метода решения и вычисление корней #
    //############################################

    // Запрашиваем у пользователя номер метода
    pointMethod method = NULL;
    backToChoice:
    puts("Choose method: ");
    puts("Print 1 to use half division method.");
    puts("Print 2 to use chord method.");
    puts("Print 3 to use Newton method.");
    puts("Print 4 to use iteration method.");
    int choice = 0;
    scanf("%d", &choice);

    // Инициализируем указатель на функцию выбранным методом
    switch (choice)
    {
    case 1:
        method = halfDivisionMethod;
        break;
    case 2:
        method = chordMethod;
        break;
    case 3:
        method = newtonMethod;
        break;
    case 4:
        method = interationMethod;
        break;
    default:
        puts("There is no such method! Choose again!\n");
        goto backToChoice;
    }

    // Решаем уравнение выбранном методом на заданных отрезках
    while (!sections.empty()) {
        try {
            // Запоминаем вычисленные корни
            roots.push_back(method(f, sections.back()));
        }
        catch (std::runtime_error& e) { // Если что-то пошло не так
            printf("Error: %s\n\n", e.what());
        }
        sections.pop_back();
    }
    puts("");



    //####################################
    //# Вывод результатов и тестирование #
    //####################################

    // Выводим все корни на экран
    printf("Roots: ");
    for (int i = 0; i != roots.size(); ++i)
    printf("%f   ", roots.at(i));
    puts("\n");

    // Проверяем корни на правильность
    double currentRoot = 0.0;
    double currentAnswer = 0.0;
    while (!roots.empty()) {
        currentRoot = roots.back();
        currentAnswer = f(currentRoot);
        if (sqr(currentAnswer) > EPS*2)
            printf("Test failed: %f is not root of the equation! f(%f) = %f\n", currentRoot, currentRoot, f(currentRoot));
        else
            printf("Test passed.\n");
        roots.pop_back();
    }
    puts("");



    //###################
    //# Конец программы #
    //###################

    puts("End of program. Press enter to exit.");
    getchar();
    getchar();
    return 0;
}



//############################################
//# Функции f(x) для уравнений вида f(x) = 0 #
//############################################

double f1(double x) {
    return exp(x) - 10 * x;
}

double f2(double x) {
    return 2 * x;
}

double f3(double x) {
    return x * x - 5 * sin(x);
}

double f4(double x)
{
    return sin(x);
}

double f5(double x)
{
    return x*x;
}


//######################################
//# Функции для нахождения производных #
//######################################

// Первая производная функции в точке
double derivative(pointFunc f, double x)
{
    //const double h = 0.0005;
    //return (f(x) - f(x + h)) / h;
    return exp(x) - 10; // For f1
    //return 2 * x - 5 * cos(x); // For f3
    //return 2 * x; // For f5
}

// Вторая производная функции в точке
double doubleDerivative(pointFunc f, double x)
{
    //const double h = 0.0005;
    //return (f(x + h) - 2 * f(x) + f(x - h)) / (h * h);
    return exp(x); // For f1
    //return 2 + 5 * sin(x); // For f3
    //return 2.0; // For f5
}

// Знак (проверяет знак числа, возвращает +1 или -1)
int sign(double x)
{
    return (x > 0) ? (1) : (-1);
}

//######################################################################
//# Функции для нахождения одного корня уравнения на указанном отрезке #
//######################################################################

// Метод половинного деления
double halfDivisionMethod(pointFunc f, const Section & section) {
    int iterationsCounter = 0;
    double a = section.leftBorder;
    double b = section.rightBorder;
    double x = (a + b) / 2.0;
    double yLeft = f(a), yRight = f(b), yCenter = f(x);

    while (fabs(yRight - yLeft) > EPS  && iterationsCounter != 10000) {
        if (yCenter * yRight < 0) {
            a = x;
            yLeft = f(a);
        }
        else {
            b = x;
            yRight = f(b);
        }

        x = (a + b) / 2.0;
        yCenter = f(x);
        ++iterationsCounter;
    }

    if (iterationsCounter == 10000)
        throw std::runtime_error("Half division method does not respond.");

    return x;
}

// Метод хорд
double chordMethod(pointFunc f, const Section & section) {
    int iterationsCounter = 0;
    double a = section.leftBorder;
    double b = section.rightBorder;
    double result = 0.0;
    double yLeft = f(a), yRight = f(b);


    if ((yLeft*doubleDerivative(f ,a)) > 0)
    {
        while (fabs(b - ((b - a) / (yRight - yLeft))*yRight - b) > EPS && iterationsCounter != 100000)
        {
            if (sqr(yRight - yLeft) <= EPS)
                throw std::runtime_error("Chord method does not converge.");

            b = b - ((b - a) / (yRight - yLeft))*yRight;
            ++iterationsCounter;
            yRight = f(b);
        }
        result = b;
    }
    else
    {
        if (sqr(yRight - yLeft) <= EPS)
            throw std::runtime_error("Chord method does not converge.");

        while (fabs(a - ((b - a) / (yRight - yLeft))*yLeft - a) > EPS && iterationsCounter != 100000)
        {
            a = a - ((b - a) / (yRight - yLeft))*yLeft;
            ++iterationsCounter;
            yLeft = f(a);
        }
        result = a;
    }


    if (iterationsCounter == 100000)
        throw std::runtime_error("Chord method does not respond.");

    return result;
}

// Метод Ньютона
double newtonMethod(pointFunc f, const Section & section)
{
    int iterationsCounter = 0;
    double a = section.leftBorder;
    double b = section.rightBorder;
    double yLeft = f(a), yRight = f(b);


    //if(sign(derivative(f, a)) != sign(derivative(f, b)) || sign(doubleDerivative(f, a)) != sign(doubleDerivative(f, b)))
    //  throw std::runtime_error("Newton method does not converge.");


    double x = a;
    if (!(yLeft * doubleDerivative(f, x) > 0))
    {
        x = b;
        if (!(yRight * doubleDerivative(f, x) > 0))
        {
            x = (a + b) / 2;
            if (!(f(x) * doubleDerivative(f, x) > 0))
                puts("Warning: lost of accuracy is possible using Newton Method!\n");
        }
    }


    double next = x;
    do
    {
        x = next;
        if(sqr(derivative(f, x)) <= EPS)
            throw std::runtime_error("Newton method does not converge.");

        next = x - f(x) / derivative(f, x);
    } while (fabs(next - x) > EPS && iterationsCounter != 100000);


    if (iterationsCounter == 100000)
        throw std::runtime_error("Newton method does not respond.");


    return next;
}

// Метод итераций
double interationMethod(pointFunc f, const Section & section)
{
    int iterationsCounter = 0;
    double x = section.leftBorder;
    double k = 0.1;


    double next = x;
    do
    {
        k = 0.1;
        while (fabs(1 + k * derivative(f, x)) > 0.5)
            k /= 10;

        x = next;
        next = x + k * f(x);
    } while (fabs(next - x) > EPS && iterationsCounter != 100000);


    if (iterationsCounter == 100000)
        throw std::runtime_error("Iteration method does not respond.");


    return next;
}