Merge branch 'master' of ssh://vo20.nl/git/uva

modsim/ass2/Makefile

@@ -9,7 +9,7 @@ q%: q%.o
 
 q2: bisection.o
 q3: bisection.o regula_falsi.o newton_raphson.o
-q4: newton_raphson.o
+q4: bisection.o regula_falsi.o newton_raphson.o
 q5: integral.o
 q6: integral.o
 

modsim/ass2/bisection.c

 #include "bisection.h"
 
 double bisec(func_ptr f, double left, double right,
-        double epsilon, unsigned int *steps) {
-    double mid, fmid;
+        double epsilon, unsigned int *steps, unsigned int max_steps) {
+    double mid = left, fmid;
 
-    for( *steps = 0; fabs(right - left) > 2 * epsilon; (*steps)++ ) {
+    for( *steps = 0; (fabs(right - left) > 2 * epsilon)
+            && (!max_steps || (*steps < max_steps)); (*steps)++ ) {
         mid = (right + left) / 2;
 
         if( f(left) * (fmid = f(mid)) < 0 )
             right = mid;
         else if( f(right) * fmid < 0 )
             left = mid;
-        else
+        else {
+            (*steps)++;
             break;
+        }
     }
 
     return mid;

modsim/ass2/bisection.h

@@ -5,6 +5,6 @@
 #include "func_ptr.h"
 
 double bisec(func_ptr f, double left, double right,
-        double epsilon, unsigned int *steps);
+        double epsilon, unsigned int *steps, unsigned int max_steps);
 
 #endif

modsim/ass2/q2.c

@@ -20,7 +20,7 @@ int main(int argc, char *argv[]) {
     end = atoi(argv[2]);
 
     for( i = begin; i <= end; i++ ) {
-        bisection = bisec(&func, 0, 2, pow(10, -1.0 * i), &steps);
+        bisection = bisec(&func, 0, 2, pow(10, -1.0 * i), &steps, 100000);
         printf("zero point: %.30f for epsilon = 1e-%d (%d steps)\n", bisection, i, steps);
     }
 

modsim/ass2/q3.c

@@ -27,7 +27,7 @@ int main(int argc, char *argv[]) {
 
     for( i = begin; i <= end; i++) {
         epsilon = pow(10, -1.0 * i);
-        tmp= bisec(&f, 1, 2, epsilon, &steps);
+        tmp= bisec(&f, 1, 2, epsilon, &steps, 100000);
         printf("Sqrt(2) using bisection:      %.20f (%d steps; epsilon=%.0e)\n",
             tmp, steps, epsilon);
 

modsim/ass2/q4.c

@@ -2,9 +2,11 @@
 #include <stdio.h>
 #include <math.h>
 #include "newton_raphson.h"
+#include "bisection.h"
+#include "regula_falsi.h"
 
 #define EPSILON 1e-11
-#define BISEC_STEPS 5
+#define MAX_STEPS 100000
 
 /*
  * Manually define the functions and their derivatives.
@@ -34,14 +36,44 @@ double df3(double x) {
     return (x * x + 1) + 2 * x * (x - 4);
 }
 
-int main(void) {
-    unsigned int steps;
-    double root;
+#define PRINT_ZERO(f, df, x_start) { \
+    if( !isnan(root = newton_raphson(&f, &df, x_start, EPSILON, &steps, MAX_STEPS)) ) \
+        printf(#f ": %.11f (%d steps from start point %.11f)\n", root, steps, (double)x_start); \
+    else \
+        printf(#f ": could not find a root after %d steps\n", steps); \
+}
+
+#define PRINT_ZERO_RANGE(f, df, method, lnbd, ubnd) { \
+    x_start = method(&f, lnbd, ubnd, EPSILON, &start_steps, max_start_steps); \
+    printf("Did first %d steps using " #method "\n", start_steps); \
+    PRINT_ZERO(f, df, x_start); \
+}
+
+int main(int argc, char *argv[]) {
+    unsigned int steps, start_steps, max_start_steps;
+    double root, x_start;
+
+    if( argc != 2 ) {
+        printf("Usage: %s START_STEPS\n", argv[0]);
+        return 1;
+    }
+
+    max_start_steps = atoi(argv[1]);
+
+    PRINT_ZERO(f1, df1, 0);
+    PRINT_ZERO(f2, df2, -10);
+    PRINT_ZERO(f2, df2, 10);
+    PRINT_ZERO(f3, df3, 0);
+
+    printf("\nNow using Bisection to determine the first %d steps.\n\n", max_start_steps);
+
+    PRINT_ZERO_RANGE(f2, df2, bisec, 0, 10);
+    PRINT_ZERO_RANGE(f3, df3, bisec, 0, 10);
+
+    printf("\nNow using Regula Falsi to determine the first %d steps.\n\n", max_start_steps);
 
-    if( !isnan(root = newton_raphson(&f2, &df2, 1000000, EPSILON, &steps, 100000)) )
-        printf("f2: %.11f (%d steps)\n", root, steps);
-    else
-        printf("f2: could not find a root after %d steps\n", steps);
+    PRINT_ZERO_RANGE(f2, df2, regula_falsi, 0, 10);
+    PRINT_ZERO_RANGE(f3, df3, regula_falsi, 0, 10);
 
     return 0;
 }

modsim/ass2/report.tex

@@ -99,7 +99,7 @@ berekening is exponentieel.}
 \noindent Uit de grafiek kunnen we aflezen dat het aantal benodigde stappen voor
 de berekening lineair toeneemt met het aantal decimalen waarop de berekening
 nauwkeurig is. Merk op dat de datapunten niet op een perfect rechte lijn
-liggen, dit komt doordat het aantal stappen wordt afgerond op gehele getallen.
+liggen, dit komt doordat het aantal stappen wordt afgerond op een geheel getal.
 
 % }}}
 
@@ -120,10 +120,10 @@ berekening voor drie verschillende ``root-finding'' methodes (blauw: Bisection,
 rood: Regula Falsi, groen: Newton-Raphson).}
 \end{figure}
 
-We zien, net als bij opgave 2, dat het aantal benodigde stappen lineair toeneemt
-met het aantal decimale waarop het resultaat is afgerond. Merk op dat ook hier
-de grafiekkrommen niet perfect recht zijn, omdat het aantal stappen is afgerond
-op gehele getallen.
+We zien, net als bij opgave 2, dat bij elke methode het aantal benodigde stappen
+lineair toeneemt met het aantal decimalen waarop het resultaat is afgerond. Merk
+op dat ook hier geen perfect rechte lijnen door de datapunten gaan, omdat het
+aantal stappen is afgerond op een geheel getal.
 
 % }}}
 
@@ -213,7 +213,8 @@ sin & 0 & $8\pi$ & gauss       & $-1.797258919631 \cdot 10^{-14}$ & $1.797258919
 \end{tabular}
 \end{table}
 
-% TODO: calculate integral
+We kunnen de integraal $\int_0^2{x^{-0.5}dx}$ als volgt exact berekenen:
+$$ \int_0^2{x^{-0.5}dx} = \left[2\sqrt{x}\right]_0^2 = 2\sqrt{2} $$
 
 % }}}