start dn4

dn4/quad.c

@@ -0,0 +1,60 @@
+// metoda za adaptivno numericno integracijo
+// prevedi s: gcc -O2 quad.c -fopenmp -lm -o quad
+
+#include <stdio.h>
+#include <math.h>
+#include <omp.h>
+
+#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;
+}