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

Язык Си в примерах

---

1. Компиляция программ
2. Простейшая программа «Hello World»
3. Учимся складывать
4. Максимум
5. Таблица умножения
6. ASCII-коды символов
7. Верхний регистр
8. Скобочки
9. Факториал
10. Степень числа
11. Треугольник Паскаля
12. Корень уравнения
13. Система счисления
14. Сортировка
15. Библиотека complex
16. Сортировка на основе qsort
17. RPN-калькулятор
18. RPN-калькулятор на Bison
19. Простая грамматика
20. Задача «Расчёт сопротивления схемы»
21. Простая реализация конечного автомата
22. Использование аргументов командной строки
23. Чтение и печать без использования stdio
24. Декодирование звукозаписи в формате ADX
25. Другие примеры

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

Алгоритм вычисления корня основан на методе деления пополам: Предположим, что искомый корень находится между l = 0 и r = 2. Найдем середину c промежутка [l, r). Корень находится на одном из промежутков: либо на [l, c), либо на [с, r), а именно, на том, значение функции на концах которого имеет разные знаки (вспомните теорему Ролля про непрерывную функцию из курса мат. анализа). Выберем нужный из двух промежутков и применим к нему такие же рассуждения. Будем осуществлять деление попалам, пока размер промежутка не станет меньше необходимой точности.

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

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