Задача. Используйте метод хорд для того, чтобы отыскать с точностью [latex]\varepsilon[/latex] все действительные корни уравнения [latex]\frac{x}{2 \cdot \sin x +1}=\tan(\ln(x^2+1))[/latex]. Для подготовки необходимых графиков воспользуйтесь этим ресурсом.
Тесты(найдено с помощью математической системы WolframAlpha):
[latex]A[/latex] | [latex]B[/latex] | [latex]x\approx[/latex] |
-20 | 20 | -11.6945378230838209122818536587051434153…
0
10.9948442206261587135425985750810372810… |
Код программы
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 |
#include <stdio.h> #include <math.h> #include <iostream> #include<cfloat> using namespace std; const int max_iterations=100; const double interval = 0.001; // interval length for roots subdivision bool is_root = false; // is false if function has a second type break in root candidate point double delta; //precision of root search double A, B; double f(double x)//our function { double r = DBL_MAX; try { r = x / (2*sin(x)+1) - tan(log(x*x+1)); } catch(...) { //nothing to do for arithmetic exception } return r; } bool is_change_sign(double fa, double fb) { return fa*fb<0; } double find_root(double a, double b, double delta, bool *yes) { double fa; double fb; double fc; double c; int l=0; fa=f(a); fb=f(b); double dy,dy1; //old and new OY projections of chord dy=fabs(fa-fb); dy1=DBL_MAX; while(true) { l++; if(l > max_iterations) { *yes=false; return (a+b)/2; } c = a - (b - a) * fa/(fb - fa); fc = f(c); if(fa*fc>0) { a=c; fa=fc; } else { b=c; fb=fc; } if(fabs(fa)<delta) { *yes=true; return a; } if(fabs(fb)<delta) { *yes=true; return b; } if(l>10) //there's a magic number, we must check conditions of break or stop after some iterations(some = 10) { if(dy1>dy) //checking if there's break of second kind(pole) { *yes=false; return (a+b)/2; } if(fabs(b-a)<delta) //else we are near the root { *yes=true; return (a+b)/2; } } dy=dy1; //change old chord height dy1=fabs(fa-fb); //calculate new chord height } // we can't get here *yes=false; return c; } int main() { scanf("%lf %lf %lf",&A,&B,&delta); double a,b; double fa, fb; /* test interval for finding roots int k=2; A= sqrt(exp(M_PI*(1/2.0+k))-1)-M_PI;//70010637.8772-10; B= sqrt(exp(M_PI*(1/2.0+k))-1)+M_PI;//70010637.8772+10; delta = 1e-12; printf("A=%20.19g A=%20.19g delta=%20.19g \n",A,B,delta); */ fb = f(A); for(a = A, b = A + interval; b<=B; a += interval, b += interval) { fa = fb; fb = f(b); if(is_change_sign(fa,fb)) { double r=find_root(a,b,delta, &is_root); if(is_root) { printf("%14.15lf \n",r); } } } return 0; } |
Алгоритм
Для начала запишем данное нам уравнение в виде функции [latex]y=f(x)[/latex] и построим ее график:
[latex]y=\frac{x}{2 \cdot sin x +1}-tan(ln(x^2+1))[/latex]
Задача о нахождении приближённых значений действительных корней уравнения [latex]f(x)=0[/latex] предусматривает предварительное отделение корня, то есть установление интервала, в котором других корней данного уравнения нет.
Метод хорд предусматривает замену функции на интервале на секущую, и поиск ее пересечения с осью [latex]OX[/latex]. На заданном интервале [latex][a,b][/latex] с точностью [latex]\varepsilon[/latex] корень будет вычисляться согласно итерационному выражению, которое имеет вид:
[latex]x_{i+1}=x_{i-1}-\frac{f(x_{i-1}) \cdot (x_{i}-x_{i-1})}{f(x_{i})-f(x_{i-1}) ) }[/latex]Данный метод имеет свои недостатки. В первую очередь видно, что он не учитывает возможную разрывность функции, вследствие чего могут возникать ложные корни, или пропадать имеющиеся. Как же выяснить, есть ли корень на данном отрезке? Рассмотрим и проанализируем случаи, в которых метод хорд может выдать ошибочный результат. Возможны следующие варианты:
- В точке, где находится предполагаемый корень, имеется разрыв второго рода. Здесь метод хорд обнаруживает перемену знака и начинает сужать отрезок. Однако расстояние между крайними точками отрезка не уменьшается, а увеличивается. А именно увеличивается проекция отрезка на ось [latex]OY[/latex]. И чем ближе мы находимся к точке разрыва, тем она больше, а в самой точке стремится к бесконечности. В качестве примера приведем функцию [latex]\frac{(x-5)^2}{x-4}[/latex], имеющую разрыв второго рода в точке [latex]x=4[/latex].
-
Аналогичный случай — точка разрыва первого рода, где наша хорда стремится к некоторой константе — величине разрыва. Пример — функция [latex]\frac{\sin x\cdot(x-2.5)}{\left | x-2,5 \right |}[/latex], имеющая разрыв первого рода в точке [latex]x=2,5[/latex].
-
Функция не является разрывной и даже имеет корень, однако в точке корня производная стремится к бесконечности. Например, функция [latex]\sqrt[3]{x-1,2}[/latex], имеющая корень в точке [latex]x=1,2[/latex]. Для функций подобного рода длина проекции отрезка на ось [latex]OY[/latex] будет очень медленно меняться, вследствие чего для разумного числа итераций она будет превосходить заранее выбранную точность [latex]\varepsilon[/latex].
- Функция равна нулю или очень близка к нулю на некотором интервале(например, функция [latex]y=rect(x-1,5)\cdot(x-1)\cdot(x-2)^2[/latex]). Здесь метод хорд найдет пересечение с осью [latex]OX[/latex] в интервале, где находятся корни(одна из сторон отрезка будет корнем), но все последующие итерации будут выдавать эту же точку. Поэтому хорда не будет уменьшаться, и даже этот один корень не будет найден(если будет использоваться стандартная [latex]\delta[/latex]-оценка точности по оси [latex]OX[/latex]).
- Функция вида [latex]y=x^{2k}[/latex] или ей подобная, например, [latex]y=1+\sin x[/latex], к которой метод хорд вообще не применим, так как нарушается начальное условие применимости этого метода. Здесь нужно ввести дополнительную проверку. Изменим значение функции на небольшую константу — нашу точность [latex]\varepsilon[/latex] и повторим процедуру поиска корней. В результате мы получим [latex]2k[/latex] корней, каждую пару из которых мы можем считать краями интервала, в котором лежит настоящий корень.
К счастью, в данной нам функции присутствует только один из этих случаев, а именно разрыв второго рода. Аналитически рассмотрев нашу функцию, мы обнаружили, что корни следует искать в окрестности точек[latex]\sqrt{e^{\frac{\pi}{2}+\pi \cdot k}-1}[/latex] с отклонением [latex]\pm\pi[/latex]. В корнях функции ее производная быстро растет с ростом [latex]k[/latex].
Критерием отбрасывания кандидата на корень будет рост длины хорды при сужении интервала. Критерием останова будет сужение интервала до заданной точности [latex]\delta[/latex].
Код программы
Молодец, проделана большая работа и получен хороший результат.
И как с каждой большой работой, есть о чём поговорить.
По оформлению.
Пожалуйста, задавайте стандартные функции в latex правильно: \sin x, а не sin x и т.п.
По коду.
— Если нужна бесконечность типа double, то резонно выбрать DBL_MAX. Задавая своё число «с потолка» Вы рискуете.
— Отступы. Я писал Вам подробный комментарий по поводу GNU-style, который Вы выбрали. Пожалуйста, придерживайтесь стандартов. Что касается дополнительного пробела для демонстрации зависимости по данным, то либо Вы к ней относитесь не очень строго, либо я не все Ваши идеи понял:
— что означает дополнительный сдвиг в строке 64?
— почему 46-я сдвинута на два пробела, а не один?
Предлагаю решить так. В процессе программирования Вы используете любые личные наработки и оригинальный стиль. Но перед публикацией делаете свою работу понятной читателю. Т.е. оформляете код в одном из существующих стилей. В этом случае Вы не потеряете своей оригинальности, но и приобретёте навыки социализации, необходимые в общении с другими программистами. Кстати, я пользуюсь системой автоматической расстановки отступов, чтобы читать Ваш код. Может и Вам она пригодится.
Согласны?
По содержанию.
Вы пишите «В общем случае функция имеет ряд промежутков, на которых это условие не выполняется». Что означает «в общем случае»? У нас одна конкретная функция, мы ничего не обобщаем.
Дальше. Нам не нужно, чтобы функция была монотонной везде. Нужно только, чтобы существовала окрестность корня в которой функция была монотонна.
По тестам.
— Тесты, это то, что ДОЛЖНА выдать программа, а не то, что она выдаёт. Вы пишите, что точка 0 является корнем функции при точности 10-6, но при точности 10-12 она уже не корень. Объясните, почему?
— Правильно было бы найти корни с помощью какой-либо математической системы:
x~ -11.6945378230838209122818536587051434153…
x~ -1.25741503276862309237205903178504130394…
x = 0
x~ 0.547316310185252929580383582338832450320…
x~ 10.9948442206261587135425985750810372810…
и сделать их эталонными тестами для проверки программы.
Что ожидалось?
Как видите количество корней бесконечно, и из-за логарифма функция не является периодической. Как же численно найти «все». Понятное дело — никак! Значит стоит «научить» программу находить отделять интервал очередного (то положительного, то отрицательного) корня и находить его без остановки программы. Т.е. она должна выдавать последовательные корни без остановки.
Вы спросите — а можно оставить как есть? Можно. Только поправьте предыдущие замечания.
Спасибо, Игорь Евгеньевич. Я исправила.
Отлично. Принято.