ModSim: Almost completed assignment 2.5

modsim/ass2/Makefile

-CC=clang
-CFLAGS=-Wall -Wextra -pedantic -std=c99 -D_GNU_SOURCE -g
+CC=gcc
+CFLAGS=-Wall -Wextra -pedantic -std=c99 -D_GNU_SOURCE -g -O0
 LFLAGS=-lm
 
 all: q1 q2 q3 q4 q5 q7 report.pdf

modsim/ass2/integral.c

@@ -20,12 +20,12 @@ double gauss(func_ptr f, double a, double b) {
     return half_dx * (f(x1) + f(x2));
 }
 
-double integral(func_ptr f, method_ptr method, double a, double b, double dx) {
-    int steps = (int)((b - a) / dx), i;
-    double x, surface = 0;
+double integral(func_ptr f, method_ptr method, double a, double b, int steps) {
+    int i;
+    double surface = 0, dx = (b - a) / steps;
 
     for( i = 0; i < steps; i++)
-        surface += method(f, x = a + i * dx, x + dx);
+        surface += method(f, a + i * dx, a + (i + 1) * dx);
 
     return surface;
 }

modsim/ass2/integral.h

@@ -10,6 +10,6 @@ double rectangle(func_ptr f, double a, double b);
 double trapezoidal(func_ptr f, double a, double b);
 double simpson(func_ptr f, double a, double b);
 double gauss(func_ptr f, double a, double b);
-double integral(func_ptr f, method_ptr method, double a, double b, double dx);
+double integral(func_ptr f, method_ptr method, double a, double b, int steps);
 
 #endif

modsim/ass2/q5.c

@@ -4,8 +4,6 @@
 #include "func_ptr.h"
 #include "integral.h"
 
-#define DX 1e-5
-
 double f1(double x) {
     return pow(M_E, -x);
 }
@@ -14,34 +12,35 @@ double f2(double x) {
     return x * f1(x);
 }
 
-#define PRINT_INTEGRAL(func, method, a, b) (printf(#func " from " #a " to " \
-    #b " using %-19s %.11e\n", #method " method:", integral(&func, &method, a, b, DX)))
+#define PRINT_INTEGRALS(func, a, b, real) { \
+    double _i; \
+    PRINT_INTEGRAL(func, a, b, rectangle, (real)); \
+    PRINT_INTEGRAL(func, a, b, trapezoidal, (real)); \
+    PRINT_INTEGRAL(func, a, b, simpson, (real)); \
+    PRINT_INTEGRAL(func, a, b, gauss, (real)); \
+}
+
+#define PRINT_INTEGRAL(func, a, b, method, real) { \
+    _i = integral(&func, &method, a, b, atoi(argv[1])); \
+    printf(#func " from " #a " to " #b " using %-19s %.80e (%.2f%%)\n", \
+        #method " method:", _i, fabs((real - _i) / real * 100)); \
+}
+
+int main(int argc, char **argv) {
+    if( argc != 2 ) {
+        printf("Usage: %s STEPS", argv[0]);
+        return -1;
+    }
 
-int main(void) {
-    PRINT_INTEGRAL(f1, rectangle, 0, 1);
-    PRINT_INTEGRAL(f1, trapezoidal, 0, 1);
-    PRINT_INTEGRAL(f1, simpson, 0, 1);
-    PRINT_INTEGRAL(f1, gauss, 0, 1);
+    PRINT_INTEGRALS(f1, 0, 1, (M_E - 1) / M_E);
     puts("");
-    PRINT_INTEGRAL(f2, rectangle, 0, 2);
-    PRINT_INTEGRAL(f2, trapezoidal, 0, 2);
-    PRINT_INTEGRAL(f2, simpson, 0, 2);
-    PRINT_INTEGRAL(f2, gauss, 0, 2);
+    PRINT_INTEGRALS(f2, 0, 2, 1 - 3/pow(M_E, 2));
     puts("");
-    PRINT_INTEGRAL(f2, rectangle, 0, 20);
-    PRINT_INTEGRAL(f2, trapezoidal, 0, 20);
-    PRINT_INTEGRAL(f2, simpson, 0, 20);
-    PRINT_INTEGRAL(f2, gauss, 0, 20);
+    PRINT_INTEGRALS(f2, 0, 20, 1 - 21/pow(M_E, 20));
     puts("");
-    PRINT_INTEGRAL(f2, rectangle, 0, 200);
-    PRINT_INTEGRAL(f2, trapezoidal, 0, 200);
-    PRINT_INTEGRAL(f2, simpson, 0, 200);
-    PRINT_INTEGRAL(f2, gauss, 0, 200);
+    PRINT_INTEGRALS(f2, 0, 200, 1 - 201/pow(M_E, 200));
     puts("");
-    PRINT_INTEGRAL(sin, rectangle, 0, 8 * M_PI);
-    PRINT_INTEGRAL(sin, trapezoidal, 0, 8 * M_PI);
-    PRINT_INTEGRAL(sin, simpson, 0, 8 * M_PI);
-    PRINT_INTEGRAL(sin, gauss, 0, 8 * M_PI);
+    PRINT_INTEGRALS(sin, 0, 8 * M_PI, 0.0);
 
     return 0;
 }

modsim/ass2/q5.py

@@ -36,7 +36,7 @@ integrals = [(lambda x: exp(-x), 0, 1),
         (lambda x: x * exp(-x), 0, 200),
         (lambda x: sin(x), 0, 8 * pi)]
 
-methods = [rectangle, trapezoidal, simpson]
+methods = [rectangle, trapezoidal, simpson, gauss]
 
 DX = 1e-3
 
@@ -44,7 +44,5 @@ for i in integrals:
     for method in methods:
         print 'integrate from %d to %d using %-11s: %.16f' \
                 % (i[1], i[2], method.func_name, integrate(method, DX, *i))
-    print 'integrate from %d to %d using %-11s: %.16f' \
-            % (i[1], i[2], 'gauss', gauss(*i))
     print '-------------'
 

modsim/ass2/report.tex

@@ -86,11 +86,94 @@ beoogde nauwkeurigheid.
 
 % }}}
 
+\section{Benadering van $\sqrt{2}$} % {{{
+
+\label{sec:benadering van sqrt 2}
+
+% }}}
+
+
 \section{Newton-Raphson} % {{{
 \label{sec:Newton-Raphson}
 
 
 
+% }}}
+
+\section{Vier numerieke integraal routines} % {{{
+\label{sec:vier numerieke integraal routines}
+
+We hebben de volgende vier methodes ge\"implementeerd om een numerieke integraal
+uit te rekenen:
+
+\begin{description}
+    \item[Rectangle rule]
+    $$rectangle(f, a, b) = (b - a) \cdot f(\frac{a + b}{2})$$
+    \item[Trapezoidal rule]
+    $$trapezoidal(f, a, b) = \frac{(b - a) \cdot (f(a) + f(b))}{2}$$
+    \item[Simpson's rule]
+    Simpson's rule kan worden berekend door middel van de bovenstaande methodes
+    ``rectangle'' en ``trapezoidal'':
+    $$simpson(f, a, b) = \frac{2 \cdot rectangle(f, a, b) + trapezoidal(f, a, b)}{3}$$
+    \item[Two-point Gauss integratie]
+    $$gauss(f, a, b) = (b - a) \frac{f(a + \frac{b - a}{2} \cdot (1 -
+    \frac{1}{\sqrt{3}})) + f(a + \frac{b - a}{2} \cdot (1 +
+    \frac{1}{\sqrt{3}}))}{2}$$
+\end{description}
+
+Naast de bovengenoemde methodes hebben we ook een integraal functie gemaakt die
+de methodes aanroept:
+
+\begin{verbatim}
+double integral(f, method, a, b, dx):
+    steps = (int) ((b - a) / dx)
+    surface = 0
+    for(i = 0; i < steps; i++):
+        surface += method(f, a + i * dx, a + (i + 1) * dx)
+    return surface
+\end{verbatim}
+
+\noindent Daarnaast hebben we drie functies gemaakt (waarbij \texttt{sin(x)}
+onderdeel is van glibc) die de drie verschillende integranten voorstellen:
+
+\begin{verbatim}
+double f1(x):
+    return pow(M_E, -x)
+double f2(x):
+    return x * f1(x)
+double sin(x);
+\end{verbatim}
+
+\noindent Dit heeft het volgende resultaat gegeven als we het programma met
+1.000.000 steps uitvoeren:
+
+\begin{verbatim}
+f1 from 0 to 1 using rectangle method:           6.32120558829e-01
+f1 from 0 to 1 using trapezoidal method:         6.32120558829e-01
+f1 from 0 to 1 using simpson method:             6.32120558829e-01
+f1 from 0 to 1 using gauss method:               6.32120558829e-01
+
+f2 from 0 to 2 using rectangle method:           5.93994150290e-01
+f2 from 0 to 2 using trapezoidal method:         5.93994150290e-01
+f2 from 0 to 2 using simpson method:             5.93994150290e-01
+f2 from 0 to 2 using gauss method:               5.93994150290e-01
+
+f2 from 0 to 20 using rectangle method:          9.99999956716e-01
+f2 from 0 to 20 using trapezoidal method:        9.99999956715e-01
+f2 from 0 to 20 using simpson method:            9.99999956716e-01
+f2 from 0 to 20 using gauss method:              9.99999956716e-01
+
+f2 from 0 to 200 using rectangle method:         1.00000000001e+00
+f2 from 0 to 200 using trapezoidal method:       9.99999999965e-01
+f2 from 0 to 200 using simpson method:           9.99999999998e-01
+f2 from 0 to 200 using gauss method:             9.99999999998e-01
+
+sin from 0 to 8 * M_PI using rectangle method:   -1.04730878509e-13
+sin from 0 to 8 * M_PI using trapezoidal method: -1.08168370335e-13
+sin from 0 to 8 * M_PI using simpson method:     -8.63149135973e-14
+sin from 0 to 8 * M_PI using gauss method:       -8.55165510887e-14
+\end{verbatim}
+
 % }}}
 
 \end{document}