@lion137: Chodzi tutaj o to aby scałkować postać funkcji PDF (Probability density function), która ma postać podaną w załączniku. Należy ją scałkować w moim przypadku w przedziale <0, 5633.107>.
Z racji tego tego, że postać funkcji gamma dla zadanych parametrów przyjmuję wartości nieskończone zostało zastosowane przekształcenie logarytmiczne (zmienna s).
@lion137: Część podanej funkcji PDF (którą podałem w poprzednim komentarzu) zawierającej funkcje gamma, została zapisana pod zmienna s (w treści wątku) i dołączona do pozostałej części funkcji PDF ( zmienna x w treści wątku), dlatego zamiast:
s = (gamma((n+1)/2) - gamma(n/2))*(1/sqrt(pi*n))
jest
s = (exp(gammaln((n+1)/2) - gammaln(n/2)))*(1/sqrt(pi*n))
@_13th_Dragon: Póki co nie znam na tyle C, aby móc to zapisać, mam napisany program w C wykorzystujący obliczanie całki dla zadanej funkcji jednak próby zapisania tej funkcji kończą się niepowodzeniem.
C 
#include <stdio.h>
#include <math.h>
int main(){
double a, b, dx, ai, bi, suma_pole, s;
int n, i, c;
c = 1089498;
s = (exp(log(gamma((c+1)/2))) - (log(gamma(c/2))))*(1/(sqrt(M_PI*c)));
double fun (double t) {
return (pow((s*(1+((t*t)/c))),(-(c+1)/2)));
}
/* Wczytanie granic przedziału oraz precyzji całkowania z konsoli. */
printf("Proszę podać dolną granicę przediału całkowania [a]:\n");
scanf("%lf", &a);
printf("Proszę podać górną granicę przediału całkowania [b]:\n");
scanf("%lf", &b);
printf("Proszę podać precyzję całkowania [n]:\n");
scanf("%d", &n);
/* Obliczenie podprzedziału */
dx=(b-a)/n;
/* Pętla sumująca pola kolejnych trapezów wykorzystywanych do obliczenia całki */
for(i=0; i<n-1; i++) {
ai=a+i*dx;
bi=ai+dx;
suma_pole+=dx*((fun(ai)+fun(bi))/2);
}
/* Wydruk wyniku na ekran */
printf("Wartość całki funkcji f(t) dla przedziału (a, b) wynosi %.6lf.\n", suma_pole);
return 0;
}
@_13th_Dragon: Wyniosłem postać funkcji fun() na zewnątrz, program się kompiluje, jednak wynik wychodzi błędny. Wiem, że wynik ma być bliski 1, natomiast wychodzi 0. Wypisałem wynik ze zmiennej s i to tu prawdopodobnie jest jakiś błąd.
@_13th_Dragon: Wciąż ten sam wynik dla s. W porównaniu z matlabem gdzie wychodzi mi liczba po przecinku to tutaj jest to ponad 3 tys.
to sprawdź części składowe:
exp(log(gamma((c+1)/2.0)))
log(gamma(c/2.0))
1/(sqrt(M_PI*c))
Poza tym:
exp(log(x))==x
Zamień
s = (exp(log(gamma((c+1)/2.0))) - (log(gamma(c/2.0))))*(1/(sqrt(M_PI*c)));
na
s = (exp(log(gamma((c+1)/2.0)) - (log(gamma(c/2.0)))))*(1/(sqrt(M_PI*c)));
albo jeszcze lepiej na
s = exp(log(gamma((c+1)/2.0)) - log(gamma(c/2.0))) / sqrt(M_PI*c);
Ogólnie, cały problem wynikał z tego, że w Matlabie zamiast zrobić porządne formatowanie, naciupałeś niestety wzorek bez żadnej spacji, bez żadnego rozbicia na osobne zmienne, żeby to było jakkolwiek czytelne.
To wklej jeszcze raz to co masz w MATLAB'ie a pod spodem to co aktualnie masz w C.
Samo obliczenie s wystarczy.
Dal uproszczenia warto zrobić funkcję gammaln()
s = exp(log(gamma((c+1)/2.0)) - log(gamma(c/2.0))) * (1 / sqrt(M_PI*c));
u = log(gamma((c+1)/2.0));
v = log(gamma(c/2.0));
wyniki C
s = 0.00054
u = 15.71018
v = 15.71018
Matlab
s = (exp(gammaln((n+1)/2) - gammaln(n/2)))*(1/sqrt(pi*n))
u = gammaln((n+1)/2)
v = gammaln(n/2)
wyniki matlab
s = 0.39894
u = 6650340.51807
v = 6650333.91403
@_13th_Dragon: Zgadza się, dziękuje :D
Jednak żeby nie było za dobrze to wciąż całka wychodzi niepoprawna i tu prawdopodobnie problem z funkcja fun() ;(
@przemek.kemezrp: W fun podnosisz ułamki do bardzo dużej liczby ujemnej, co daje nieskończoność.
@lion137: @_13th_Dragon @enedil
Dziękuje za każdą pomoc, było jeszcze kilka błędów, które udało mi się rozwiązać i program działa zgodnie z oczekiwaniami:
#include <stdio.h>
#include <math.h>
double fun (double t, int c) {
double wynik;
double s = exp((gamma((c+1)/2.0)) - (gamma(c/2.0))) * (1 / sqrt(M_PI*c));
double o = (1+((t*t)/c));
double l = (-(c+1)/2);
wynik = s*(pow(o, l));
return wynik;
}
int main(){
double a, b, dx, ai, bi, suma_pole, i, n;
int c = 1089498;
/* Wczytanie granic przedziału oraz precyzji całkowania. */
a = 0;
b = 5633.107;
n = 5633107;
/* Obliczenie podprzedziału */
dx=(b-a)/n;
suma_pole = 0.0;
for(i=0.0; i<n; i++) {
ai=a+i*dx;
bi=ai+dx;
suma_pole+=dx*((fun(ai, c)+fun(bi, c))/2);
}
/* Wydruk wyniku na ekran */
printf("Wartość całki funkcji f(t) dla przedziału (a, b) wynosi %.6lf.\n", suma_pole);
return 0;
}
Tak na przyszłość - jeśli będziesz przenosić kod z Matlaba/Octave do C++ - zerknij na biblioteczkę Armadillo.
pow(), liczba Pi teżM_PI, oraz dołączam do uwagi od @lion137. Nie możesz definiować funkcji wewnątrz innej funkcji.