Язык Си в примерах/Корень уравнения

Задача: Рассмотрим трансцендентное уравнение, которое не имеет явной формулы для корня: ${\displaystyle e^{x}=x+2.\,\!}$ С помощью компьютера необходимо найти положительный корень этого уравнения с погрешностью по оси ${\displaystyle x}$ не более ${\displaystyle 10^{-10}\,\!}$.

Вот решение этой задачи:

 /* root.c: Вычисление корня трансцендентного уравнения */
 #include <stdio.h>
 #include <math.h>
 #define EPS 1e-10            // точность результата
 double f(double x) { 
     return exp(x) - 2 - x;
 }
 int main() {
    double l = 0, r = 2, c;
    while( r - l > EPS ) {
       c = ( l + r ) / 2;      // вычисляем середину отрезка;
       if( f(c) * f(r) < 0 )   // узнаем, в какой из частей
           l = c;              // находится искомый корень;
       else
           r = c; 
    }
    printf ("%.10lf\n", (l + r)/2 ); // выводим результат
 }

При компиляции этой программы с помощью GCC следует указать опцию -lm, которая указывает что при компоновке программы необходимо подключить библиотеку libm с математическими функциями:

bash$ gcc -lm root.c -o root

В библиотеке mlib, в частности, определена функция exp, вычисляющая экспоненту, а также многие другие математические функции: тригонометричекие (sin, cos, asin, acos, tan, ...), корень (sqrt), степень (pow), логарифм (log), ...

Директива #define EPS 1e-10 означает: в тексте программы идентификатор EPS заменить на число 1e-10, то есть ${\displaystyle 1\cdot 10^{-10}\,\!}$. Число EPS — это погрешность по оси ${\displaystyle x}$, с которой мы хотим найти корень.

Алгоритм вычисления корня основан на методе деления пополам: Предположим, что искомый корень находится между l = 0 и r = 2. Найдем середину c отрезка [l, r]. Корень находится на одном из отрезков: либо на [l, c], либо на [с, r], а именно, на том, значение функции на концах которого имеет разные знаки (вспомните теорему Ролля про непрерывную функцию из курса мат. анализа). Выберем нужный из двух отрезков и применим к нему такие же рассуждения. Будем осуществлять деление пополам, пока размер отрезка не станет меньше необходимой точности. (В методе деления отрезка пополам отрезок на оси ${\displaystyle x}$ делится пополам пока ${\displaystyle |f(x_{c})|>\varepsilon _{f}}$, в приведённом же методе отрезок на оси ${\displaystyle x}$ может достичь заданной величины ${\displaystyle \varepsilon _{x}}$, а значения функций ${\displaystyle f(x)}$ (особенно крутых) на оси ${\displaystyle y}$ могут очень далеко отстоять от нуля, при пологих же функциях ${\displaystyle f(x)}$ этот метод приводит к большому числу лишних вычислений.)

• За один шаг длина отрезка [l, r] уменьшается в два раза. Сколько нужно шагов, чтобы уменьшить отрезок более, чем в 1000 раз?
• Сколько требуется шагов, чтобы начиная с отрезка длины ${\displaystyle r-l=2\,\!}$ дойти до отрезка длины меньше ${\displaystyle 10^{-10}\,\!}$? Сколько требуется шагов, чтобы найти корень с точностью до 100 знаков после запятой?
• В случае деления пополам у нас есть нижняя и верхняя граница для значения корня. С каждым шагом эти границы сближаются. В методе Ньютона нахождения корня уравнения у нас имеется одно число x — текущее приближение корня. И следующее приближение получается по следующему алгоритму: находим точку на графике с абсциссой x и проводим из неё касательную к графику; абсцисса точки пересечения касательной с осью абсцисс будет новым значением x. Так делается до тех пор, пока новое x отличается от старого на число меньше, чем ${\displaystyle 10^{-10}/2\,\!}$. Реализуйте этот алгоритм. Для этого вам понадобится определить еще одну функцию, которая возвращает значение производной ${\displaystyle f'(x)=(e^{x}-x-2)'=e^{x}-1\,\!}$.

В рассмотренной программе вычисляется нуль вполне конкретной функции: f(x) = exp(x) - 2 - x. Следуя идеологии Code reuse (повторное использование кода), полезно сделать функцию, которая умеет находить нули произвольных функций. Для этого нужно научиться передавать функцию в качестве аргумента --- это возможно, и совсем несложно:

 /* Универсальная функция вычисления корня уравнения f(x) = 0  */
 #include <stdio.h>
 #include <math.h>
 #define EPS 1e-16
 double root(double l, double r, double (*f)(double)) {
    double c;
    while( r - l > EPS ) {
       c = ( l + r ) / 2;
       if( f(c) * f(r) < 0 )
           l = c;
       else
           r = c;
    }
    return l;
 }
 double f1(double x) {return cos(x) - 3 * x; }
 double f2(double x) {return exp(x) - x - 2; }
 int main() {
     printf("root1 = %lf\n", root(0, 2, f1)); // вычисляем и выводим корень уравнения f1(x) = 0
     printf("root2 = %lf\n", root(0, 2, f2)); // вычисляем и выводим корень уравнения f2(x) = 0
     return 0;
 }

Вывод программы выглядит следующим образом:

root1 = 0.316751
root2 = 1.146193

• Библиотека complex