Как реализовать трапециевидную интеграцию с бесконечными пределами в C?

Я пытаюсь создать программу c для интеграции sin(x)/sqrt(x) между 0 и бесконечностью. Я использую правило трапеции, отрезая конечные точки, когда функция стремится к бесконечности.

Однако общая сумма возврата слишком высока, и я не уверен, почему. Вот код:

#include<math.h>
#include<stdio.h>

double func(double u)
{ double a;
  a = ((sin(u))/(sqrt(u)));
  return a;}

void main()
    {
    int i, N;
    double sum, u, a, b, h, Fa, Fb, F;

    printf("Enter value of N\n");
    scanf("%d" ,&N);

    a=0.01;
    b=1000;

    h=(b-a)/(N-1);

    sum=0;
    F=func(a);
    u=a;

    for(i=0; i<N; i++)
    {
            sum=sum+F;
            u=u+h;
            F=fabs(func(u));
    }

    Fa=func(a);
    Fb=func(b);

    sum=sum-(0.5*Fa)-(0.5*Fb);
    sum=sum*h;

    printf("I: %lf\n", sum);
}

Какие-нибудь мысли?


person CreditChris    schedule 16.03.2015    source источник
comment
Каковы входные данные, каковы ожидаемые и фактические результаты?   -  person Eugene Sh.    schedule 17.03.2015
comment
Ожидаемое значение: 1,2533   -  person CreditChris    schedule 17.03.2015
comment
Входные данные — это количество используемых трапеций, поэтому чем больше входных данных, тем ближе к намеченному результату мы должны быть.   -  person CreditChris    schedule 17.03.2015
comment
Вы не проверяете возвращаемое значение scanf(), поэтому может случиться так, что ваш код вызывает неопределенное поведение.   -  person Iharob Al Asimi    schedule 17.03.2015
comment
Немного странно, что это было опубликовано - stackoverflow.com/questions/29086951/integral-приближение   -  person Ed Heal    schedule 17.03.2015
comment
Для ввода = 100, вывода = 38,5 Что, по-видимому, связано с суммированием цикла for, достигающим слишком высокого значения.   -  person CreditChris    schedule 17.03.2015


Ответы (2)


рабочий пример: http://ideone.com/Xibrov

просто удалите fabs в строке F=fabs(func(u));.

и вы должны использовать int main(void) и return 0; в конце вместо void main().

person mch    schedule 16.03.2015
comment
Привет Мч, N=10000 в твоем примере работает! Но для N = 100 я получаю -0,31, это просто ошибка из-за слишком низких значений N? - person CreditChris; 17.03.2015
comment
Вы должны отметить, что вы интегрируете от 0,01 до 1000, а не от 0 до бесконечности. Это и маленькая N работают вместе, чтобы создать неправильный ответ. - person mch; 17.03.2015

Ваша проблема на линии

F=fabs(func(u));

То, что вы действительно хотите для интеграции func, это

F=func(u);

В этом случае я считаю, что проблема в том, что fabs(func(u)) не интегрируется на бесконечности, поэтому ваш алгоритм будет расходиться.

person Degustaf    schedule 16.03.2015
comment
Привет, Дегустаф! Удаление fabs возвращает -0,31 для N=100. - person CreditChris; 17.03.2015
comment
@CreditChris Это не особенно удивительно, так как Ошибка кубически растет в b-a и уменьшается как квадрат в количестве шагов. - person Degustaf; 17.03.2015
comment
@CreditChris Кроме того, вы используете триггерную функцию с периодом 2\pi (около 6,283) и производите выборку только каждые 10. Вы пропустите большую часть колебаний. Я бы установил N=1000 как минимум. - person Degustaf; 17.03.2015