Задача состоит в численном поиске нулей функции.

О функции известно только то, что она непрерывна и достаточно гладка. Вид и поведение функциив окрестности нуля неизвестно.

Пока делаю следующее - прохожу функцию с большим шагом и для интервала поиска нахожу окрестности нуля. Дальше для каждого найденного интервала применяю метод деления отрезка пополам (ньютоновские методы циклят или вылетают за пределы окрестности ).

Плюс анализ пиков на нули - из-за большого шага приходится заниматься и этим.

Может есть более достойный алгоритм?

ничего достойного для столь широкого класса функций, думаю, найти не удастся

Цитата

Пока делаю следующее - прохожу функцию с большим шагом и для интервала поиска нахожу окрестности нуля.

ну представь себе функцию (x-sqrt(2))^2 (любая производная существует и непрерывна )
так как компьютер оперирует только рациональными числами, то при проходе с любым шагом функция даст положительное значение, а значит, промежуток с нулем найден не будет
чтобы можно было разработать нормальный алгоритм, придется значительно сузить класс функций...

Да ну нет. Я же не ставлю перед собой задачу найти ТОЧНОЕ решение.

Как раз для указанного примера и занимаюсь отдельной обработкой пиков\впадин.

Да и про функции тоже не так уж мало известно.

Немного расскажу о задаче.

Программа расчитывает углы между планетами на заданный момент времени для заданной точки на Земле. Это прямой расчет. Он запрограммирован и работает нормально (быстро и достаточно точно).

А обратная задача заключается в поиске всех моментов времени, для которых угол между двумя заданными планетами принимает указанную величину. Для данной точки на Земле.

Реально углы изменяются не резко (строил графики и смотрел). Аналитически задачу решить наверно можно, но в алгоритме прямого расчета уравнения движения планет аппроксимируются рядами, которые к тому же учитывают какие-то возмущения (я не силен в астрофизике). Поэтому копаться с этими уравнениями мне кажется бессмысленно.

Так вот, сначала прохожу весь интервал поиска с большим шагом - скажем, 1 день (максимальный интервал - 100 лет, 36500 итераций максимум). А затем методом деления отрезка пополам уточняю время, скажем, до минуты.

А угол рассчитывается с точностью, скажем, до секунды.

Так вопрос остается - а можно ли как-то оптимизировать алгоритм? И вообще как численно находятся все нули функции (в смысле - если их на интервале несколько)?

Цитата

уравнения движения планет аппроксимируются рядами

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

кстати, приведи приблизительный вид функций, посмотрим, подумаем...

Собрался было написать и назадавать вопросов, но вроде все понял по ходу.

Действительно, алгоритм сработает для полиномов.

Но возникает вопрос - а не слишком ли все сложно? Расчет производной на каждом шаге - достаточно сложная операция.

Хотя, с другой стороны, точно найдем все корни.

Но неважно. В моей задаче - тригонометрические ряды. Откуда они берутся и какой имеют вид - сейчас точно не скажу (книжки нет под рукой ). Да и это не важно. Сам угол для двух планет рассчитывается следующим образом:

Gamma(alpha1, alpha2, dec1, dec2) = arccos(cos(alpha1 - alpha2) * cos(dec1) * cos(dec2) + sin(dec1) * sin (dec2)),

где alpha1, alpha2 - так называемые "склонения" планет -
dec1, dec2 - "прямые восхождения".

(Это если рассматривать сферическую систему координат, то положение планет определяется этми двумя углами и расстоянием до планеты.)

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

А по обратному алгоритму, если задан угол Gamma0, то я ищу нули функции

F(T) = Gamma0 - Gamma(T)

описанным выше способом.

Смотри литературу по Численным методам, главы, посвященные численному решению уравнений.
К примеру, книга Демидовича "Численные методы".

Насколько я помню, смотреть нужно раздел "Алгоритмы оптимизации". Там рассматриваются алгоритмы поиска экстремумов и нулей функции...

А как на счет метода касательных и/или хорд? Или не в тему? У меня находил достаточно быстро. Могу кинуть код. Вспомнил: эта прога ищет для полинома 3-ей степени. ее нужно переделать под конкретную функцию.

Цитата

#include <iostream.h> #include <math.h> double a, b, c, d; //Coefficients for equation int SN; //Solution number int RF=0; //Roots found int NR; //Number of roots int nd; //Number of true digits in result double e1, e2; /* This function returns maximum of two numbers x and y */ double min(double x, double y) { if (y<x) {   return y; } else {   return x; } } /* This function shows found root and increases found roots number */ void declare_root(double root) { cout<<"X"<<RF++<<"="<<root<<endl; } /* The result of this function is polinom of 3rd power */ double f(float x) { return /*f(x)=*/a*pow(x,3)+b*pow(x,2)+c*x+d; } /* Derivative of function */ double ft(double x) { return /*f'(x)=*/3*a*pow(x,2)+2*b*x+c; } /* Second derivative of the function */ double ftt(double x) { return /*f''(x)=*/6*a*x+2*b; } /* This function returns a absolute error for calculating root between e1 and e2 with 10^-nd digits */ double error(double r) { // r - approximate root for function double n, M; n=-b*b/3+c; if ((e1<n)&&(n<e2)) {   M=n; } else {   M=min(ft(e1),ft(e2)); } n=M; M=ftt(e2)/(n*n*n); return e=(f®*f®/2)*M; // return the error } /* Function searches for root from point-Xo using tangent method */ double root(double Xo) { while (error(Xo)<pow(10, -1*nd) {  //while the root not found   Xo=Xo-f(Xo)/ft(Xo); //get next approach } return Xo; } /* This function gets Xo for tangent method. ONLY FOR 6th solution!!! */ double getXo_6(double max, double min) { // Let's find twist-point // From the second derivative f''(x)=6*a*x+2*b=0 double Xt=-b/(3*a); // twist-point if (f(Xt)!=0) { // twist-point is not a root?   e1=Xt-1; // left border for root calculation   e2=Xt+1; // right border for root calculation   return Xt;  // if so, return Xo for tangent method } else {   declare_root(Xt); // So, we found root! Just print it   return max; // and return max as flag, to show, that the root is already found } } /* This function gets Xo for tangent method. ONLY FOR 2nd, 4th and 6th solutions!!! */ double getXo_246(double max, double min) { while (f(max-1)>0) { // while sign of function value is positive   max-=1;        // try next value } e2=max; // right border for root calculation e1=max-1;// left border for root calculation return max-1; //return Xo for tangent method } /* This function gets Xo for tangent method. ONLY FOR 1st, 3rd and 6th solutions!!! */ double getXo_136(double max, double min) { while (f(min+1)<0) { // while sign of function value is negative   min+=1;        // try next value } e1=min; // left border for root calculation e2=min+1; // right border for root calculation return min+1; // return Xo } void main(void) { double x1, x2, des, y1, y2, Xo, res; //Get the coefficients from user cout<<"\nEnter coefficients for function ax^3+bx^2+cx+d: "; cin>>a>>b>>c>>d; cout<<"\nEnter number of digits after point: "; cin>>nd; //If a=0, this is a square equation. if (a==0) {   cout<<"\nThis is a square equation";   return; } //Divide equation by a: b/=a; c/=a; d/=a; a=1; //Calculating number of roots for the equation: // (1) - discremenant for derivative of function des=(b*b-3*a*c)/(a*a); // (2) - calculating roots for derivative (points of maximum and minimum): x1=(-b/a-sqrt(des))/3; x2=(-b/a+sqrt(des))/3; // (3) - calculating function value in these points: y1=f(x1); y2=f(x2); //Some decisions about number of roots and searching for some of them. if (y1<0) { //If y1<0 then there is one root, which greater than x2   Xo=getXo_136(x1, x2); // Calculate Xo for tangent method   res=root(Xo); // Calculate 1st root   declare_root(res); // Show 1st root   SN=1;  // 1st solution   NR=1;  // There is only 1 root } else {   if (y2>0) { //If y2>0 then there is one root, which smaller than x1   Xo=getXo_246(x1, x2); // Calculate Xo for tangent method   res=root(Xo); // Calculate 1st root   declare_root(res); // Show 1st root   SN=2;  // 2nd solution   NR=1;  // There is only 1 root   } else {   if ((y1==0)&&(y2<0)) { // If condition is true, we found 2 of 3 roots     declare_root(x1); //Show the 1st root     declare_root(x1); //Show the 2nd root     Xo=getXo_136(x1, x2); // Calculate Xo for tangent method     res=root(Xo); // Calculate 3rd root     declare_root(res); // Show 3rd root     SN=3; // 3rd solution     NR=3; // But there is one more... mother fucker!   } else {     if ((y1>0)&&(y2==0)) { // If this condition is true, we found 2 of 3 roots too     declare_root(x2); //Show the 1st root     declare_root(x2); //Show the 2nd root     Xo=getXo_246(x1, x2); // Calculate Xo for tangent method     res=root(Xo); // Calculate 3rd root     declare_root(res); // Show 3rd root     SN=4; // 4th solution     NR=3;     } else {     if ((y1==y2)&&(y2==0)) { // If this condition is true, we found...       declare_root(x1); //Show the 1st root       declare_root(x1); //Show the 2nd root       declare_root(x1); //Show the 3rd root       SN=5; // 5th solution       NR=3; // There are three roots.     } else {       if ((y1>0)&&(y2<0)) { // Otherwise there are three roots:                         // first smaller than X1, second between X1 and X2, third greater than X2       Xo=getXo_246(x1, x2); // Calculate Xo for tangent method       res=root(Xo); // Calculate 1st root       declare_root(res); // Show 1st root       Xo=getXo_6(x1, x2); // Calculate Xo for tangent method       if (Xo!=x1) { // if Xo is already root, then we havn't print it, it printed by getXo_6 function         res=root(Xo); // otherwise calculate 2nd root         declare_root(res); // Show the found 2nd root       }       Xo=getXo_136(x1, x2); // Calculate Xo for tangent method       res=root(Xo); // Calculate 3rd root       declare_root(res); // Show 3rd root       SN=6; // Solution 6th       NR=3; // There are three roots. This is worst accident       }     }     }   }   } } // Just some information about calculating cout<<endl<<RF<<" of "<<NR<<" roots found by "<<SN<<" solution."; }

Можно делить не методом дихтомии а методом золотого сечения, сходимость растет быстрее вроде.

Пардон за оффтопик.... смахивает на разработку астрологической проги (вариант - астрономической). Я права?