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

Задача: Рассмотрим трансцендентное уравнение, которое не имеет явной формулы для корня: $e^x=x+2.$ С помощью компьютера необходимо найти положительный корень этого уравнения с погрешностью по оси $x$ не более $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, то есть $1\cdot 10^{-10}$. Число EPS — это погрешность по оси $x$, с которой мы хотим найти корень.

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

Вопросы и задачи

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