Как да приложим трапецовидна интеграция с безкрайни граници в C?

Опитвам се да създам c програма за интегриране на sin(x)/sqrt(x) между 0 и Infinity. Използвам правилото на трапеца, като отрязвам крайните точки, тъй като функцията клони към безкрайност.

Общата върната сума обаче е твърде висока и не знам защо. Ето кода:

#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-approximation   -  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
Здравейте Mch, 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
Здравей Degustaf, премахването на 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