[Zelepuk 0]

Хочу найти функцию вычисления определённого интеграла для произвольного массива чисел.

В сети полно примеров взятия интеграла, например:

 

 double integral(a, b, f, eps)
           double a, b;    /* концы отрезка */
           double eps;     /* требуемая точность */
           double (*f)();  /* подынтегральная функция */
   {
           register long i;
           double fab = (*f)(a) + (*f)(B); /* сумма на краях */
           double h, h2;   /* шаг и удвоенный шаг */
           long n, n2;     /* число точек разбиения и оно же удвоенное */
           double Sodd, Seven;  /* сумма значений f в нечетных и в
                                   четных точках */
           double S, Sprev;/* значение интеграла на данной
                              и на предыдущей итерациях */
           double x;       /* текущая абсцисса */

           niter = 0;
           n = 10L;        /* четное число */
           n2 = n * 2;

           h = fabs(b - a) / n2;   h2 = h * 2.0;

           /* Вычисляем первое приближение */
           /* Сумма по нечетным точкам: */
           for( Sodd = 0.0, x = a+h, i = 0;
                                     i < n;
                                     i++, x += h2 )
                   Sodd += (*f)(x);

           /* Сумма по четным точкам: */
           for( Seven = 0.0, x = a+h2, i = 0;
                                       i < n-1;
                                       i++, x += h2 )
                   Seven += f(x);

           /* Предварительное значение интеграла: */
           S = h / 3.0 * (fab + 4.0 * Sodd + 2.0 * Seven );
           do{
                   niter++;
                   Sprev = S;

                   /* Вычисляем интеграл с половинным шагом */
                   h2 = h;      h /= 2.0;
                   if( h == 0.0 ) break;   /* потеря значимости */
                   n  = n2;     n2 *= 2;

                   Seven = Seven + Sodd;
                   /* Вычисляем сумму по новым точкам: */
                   for( Sodd = 0.0, x = a+h, i = 0;
                                             i < n;
                                             i++, x += h2 )
                   Sodd += (*f)(x);

                   /* Значение интеграла */
                   S = h / 3.0 * (fab + 4.0 * Sodd + 2.0 * Seven );

           } while( niter < 31 && fabs(S - Sprev) / 15.0 >= eps );
           /* Используем условие Рунге для окончания итераций */

           return ( 16.0 * S - Sprev ) / 15.0 ;
           /* Возвращаем уточненное по Ричардсону значение */
   }

 

Непонятно следущее: f - указывает на подынтегральную функцию. Математически это понятно. Но в Си же это подпрограмма! Я наверное принципиально ен понимаю чего-то но как можно интеграл от подпрограммы взять? или программа должна возвращать указатель на массив?

я же надеялся найти реализацию где в качестве указателя будет указатель на массив, поскольку функция представлена в виде набора отсчётов (массива) от a до b.

 

Заранее спасибо за разъяснения.

Изменено 14 декабря, 2012 пользователем Zelepuk

[MrYuran 16]

double f() - это функция, возвращающая значение в конкретной точке.

Все в соответствии с канонами.

[Zelepuk 0]

double f() - это функция, возвращающая значение в конкретной точке.

Все в соответствии с канонами.

непонятно как производится "перебор" этих точек. Нужно же посчитать для группы точек. Или в функции f() после каждого вызова должен быть инкремент указателя?

Думаю намного проще было бы сделать что-то вроде цикла где перебирались бы элементы массива.

Я полагаю конструкции вида (*f)(a) можно просто заменить на указатель на массив(с тем же индексом)?

Изменено 14 декабря, 2012 пользователем Zelepuk

[xemul 0]

непонятно как производится "перебор" этих точек. Нужно же посчитать для группы точек. Или в функции f() после каждого вызова должен быть инкремент указателя?

Думаю намного проще было бы сделать что-то вроде цикла где перебирались бы элементы массива.

Вы в код, Вами же выложенный заглядывали? Там есть кучка обращений (*f)(x), и именно в циклах.

Я полагаю конструкции вида (*f)(a) можно просто заменить на указатель на массив(с тем же индексом)?

Ваша задача - нарисовать f(x), т.е. собственно подынтегральную функцию. Будет ли функция представлена аналитически или какой-либо интерполяцией - без разницы, но для любого x на [a, b] она должна вернуть осмысленное значение.

[DRUID3 0]

Я полагаю конструкции вида (*f)(a) можно просто заменить на указатель на массив...

Ну а если она не массивом задана?

 

(с тем же индексом)?

??

[Палыч 6]

непонятно как производится "перебор" этих точек. Нужно же посчитать для группы точек.

Для этого в программе присутствуют циклы. Аргумент вызываемой функции в каждом цикле увеличивается на величину шага

x += h2

затем вычисляется значение функции в точке х и вычисляется сумма значений

S... += (*f)(х)

[demiurg_spb 0]

Хочу найти функцию вычисления определённого интеграла для произвольного массива чисел.

Может я чего-то не понимаю?

В чём смысл (достоинство) приведенного вами кода для дискретной функции?

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

[DRUID3 0]

Может я чего-то не понимаю?

В чём смысл (достоинство) приведенного вами кода для дискретной функции?

Точность... Ошибка вычисления... Функция на входе то ожидается аналитическая.

[demiurg_spb 0]

Точность... Ошибка вычисления... Функция на входе то ожидается аналитическая.

Это-то понятно, но у него не аналитическая функция, а массив. И мне думается что этот алгоритм вообще не применим для его случая.

[Xenia 35]

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

 

Интересно, откуда всякий раз выползает предубеждение, что метод прямого суммирования дискрет менее точен, чем сумма трапеций, а последний более трудоёмок?

 

Ведь если точки идут по сетке через равные интервалы, то площадь отдельной трапеции равна полусумме соседних отсчетов (F_i-1 + F_i)/2 или (что тоже самое!) сумме половинок F_i-1/2 + F_i/2. Но когда начнём все эти трапеции друг с дружкой складывать, то соседние половинки объединятся в целое! В результате чего получится простая сумма дискрет.

 

Впрочем, небольшая разница тут всё же есть - на самых краях, там, где половинок не хватает. Поэтому метод трапеций отличается от метода прямого суммирования только величиной (F₀ + F_n)/2, где F₀ - начальная точка, а F_n - конечная. Разница в трудоёмкости, как видим, здесь совершенно ничтожна. А те, у кого она велика, занимаются избыточным суммированием и дележкой на 2 на каждом обороте цикла.

[demiurg_spb 0]

Так я прямо и написал об этом что чуть лучше и чуть трудозатратнее. Никакого заблуждения как и предубеждения нет:-)

 

Разница в трудоёмкости, как видим, здесь совершенно ничтожна.

Всё от задачи зависит. Вот если массив даблов да ещё и на 8-ми битнике.... я бы разницу точно ощутил:)

[_Ivana 0]

Не, ну можно поискать приключений, вылить воду и свести к предыдущей задаче. Проинтерполировать на каждом интервале многочленом n-й степени и брать интеграл на нем, рассчитывая разность первообразных по формуле Ньютона-Лейбница :)

[Xenia 35]

Не, ну можно поискать приключений, вылить воду и свести к предыдущей задаче. Проинтерполировать на каждом интервале многочленом n-й степени и брать интеграл на нем, рассчитывая разность первообразных по формуле Ньютона-Лейбница :)

 

Тогда demiurg_spb на своём 8-битнике просто зашьётся. :)

 

[_Ivana 0]

А кто его знает, с какой частотой у автора идут отсчеты его функции? Он об этом скромно умолчал. Может она близка к частоте Найквиста (при условии ограниченности спектра функции) и остается только использовать суровые методы аппроксимации между отсчетами :) Кстати, если "произвольные отсчеты" идут ещё реже, то задача становится ещё интереснее. Или, если у него отсчеты какой-нибудь функции Дирихле, к примеру :)

[Xenia 35]

А кстати, а чего это народ функцию от указателя так вычурно вычисляет?

S += (*f)(х);

На мой взгляд (а мои взгляды могут расходиться со взглядами компилятора :)), такую функцию следует вызывать как обычную

S += f(х);

Тем более что прецедент уже есть - существует два разных способа объявления массива:

int array[100], *parray=array;

И тем не менее, обоими массивами пользуются одинаково, когда достают из них элементы:

S += array[x];

или

S += parray[x];

а вовсе не (*parray)[x];

И все только потому, что array - тоже указатель, только константный. А, значит, и оператор [] будет их разыменовывать правильно в обоих случаях.

 

По аналогии с этим, и функцию по указателю можно вызывать напрямую, как f(x), поскольку обычная функция (не указатель) тоже является указателем на функцию, только константным. Значит, и оператор () тоже должен работать в обоих случаях одинаково.

Или я не права?