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

Материал из Викиучебника

Задача: Рассмотрим трансцендентное уравнение, которое не имеет явной формулы для корня: С помощью компьютера необходимо найти положительный корень этого уравнения с погрешностью не более .

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

 /* 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, то есть . Число EPS — это погрешность, с которой мы хотим найти корень.

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

[править] Вопросы и задачи

• За один шаг длина промежутка [l, r) уменьшается в два раза. Сколько нужно шагов, чтобы уменьшить отрезок более, чем в 1000 раз?
• Сколько требуется шагов, чтобы начиная с отрезка длины дойти до отрезка длины меньше ? Сколько требуется шагов, чтобы найти корень с точностью до 100 знаков после запятой?
• В случае деления попалам у нас есть нижняя и верхняя граница для значения корня. С каждым шагом эти границы сближаются. В методе Ньютона нахождения корня уравнения у нас имеется одно число x — текущее приближение корня. И следующее приближение получается по следующему алгоритму: находим точку на графике с абциссой x и проводим из неё касательную к графику; абcцисса точки пересечения касательной с осью абцисс будет новым значением x. Так делается до тех пор, пока новое x отличается от старого на число меньше, чем . Реализуйте этот алгоритм. Для этого вам понадобится определить еще одну функцию, которая возвращает значение производной .

[править] Универсальная функция вычисления корня

В рассмотренной программе вычисляется нуль вполне конкретной функции: 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