diff --git a/dn4/quad.c b/dn4/quad.c new file mode 100644 index 0000000..f24519e --- /dev/null +++ b/dn4/quad.c @@ -0,0 +1,60 @@ +// metoda za adaptivno numericno integracijo +// prevedi s: gcc -O2 quad.c -fopenmp -lm -o quad + +#include +#include +#include + +#define TOL 1e-8 + +double func(double x) +{ + return sin(x*x); +} + +double quad(double (*f)(double), double lower, double upper, double tol) +// funkcija za integracijo, spodnja meja, zgornja meja, dovoljena napaka +{ + double quad_res; // rezultat + double h; // dolzina intervala + double middle; // sredina intervala + double quad_coarse; // groba aproksimacija + double quad_fine; // fina aproksimacija (two trapezoids) + double quad_lower; // rezultat na spodnjem intervalu + double quad_upper; // rezultat na zgornjem intervalu + double eps; // razlika + + h = upper - lower; + middle = (lower + upper) / 2; + + // izracunaj integral z obema aproksimacijama trapezoidnega pravila + quad_coarse = h * (f(lower) + f(upper)) / 2.0; + quad_fine = h/2 * (f(lower) + f(middle)) / 2.0 + h/2 * (f(middle) + f(upper)) / 2.0; + eps = fabs(quad_coarse - quad_fine); + + //ce se nismo dosegli zelene natancnosti, razdelimo interval na pol in ponovimo + if (eps > tol) + { + quad_lower = quad(f, lower, middle, tol / 2); + quad_upper = quad(f, middle, upper, tol / 2); + quad_res = quad_lower + quad_upper; + } + else + quad_res = quad_fine; + + return quad_res; +} + +int main(int argc, char* argv[]) +{ + double quadrature; + double dt = omp_get_wtime(); + + quadrature = quad(func, 0.0, 50.0, TOL); + + dt = omp_get_wtime() - dt; + + printf("Integral: %lf\nCas: %lf\n", quadrature, dt); + + return 0; +}