Implemented Gauss correctly and updated all header files.

modsim/ass2/Makefile

@@ -7,14 +7,15 @@ all: q1 q2 q3 q4 q5 q7 report.pdf
 q%: q%.o
 	$(CC) $(CFLAGS) $(LFLAGS) -o $@ $^
 
-%.o: %.c
-	$(CC) $(CFLAGS) -o $@ -c $^
+q2: bisection.o
+q3: bisection.o regula_falsi.o newton_raphson.o
+q4: newton_raphson.o
+q5: integral.o
 
 %.pdf: %.tex
 	pdflatex $^
 	pdflatex $^
 
-.PHONY: clean
 clean:
 	for i in `seq 7`; do \
 		rm -vf q$$i; \

modsim/ass2/bisec.h → modsim/ass2/bisection.c

-#include "func_ptr.h"
+#include "bisection.h"
 
 double bisec(func_ptr f, double left, double right,
         double epsilon, unsigned int *steps) {

modsim/ass2/bisection.h

+#ifndef BISECTION_H
+#define BISECTION_H
+
+#include <math.h>
+#include "func_ptr.h"
+
+double bisec(func_ptr f, double left, double right,
+        double epsilon, unsigned int *steps);
+
+#endif

modsim/ass2/integral.c

+#include "integral.h"
+
+double rectangle(func_ptr f, double a, double b) {
+    return (b - a) * f((a + b) / 2);
+}
+
+double trapezoidal(func_ptr f, double a, double b) {
+    return .5 * (b - a) * (f(a) + f(b));
+}
+
+double simpson(func_ptr f, double a, double b) {
+    return (2 * rectangle(f, a, b) + trapezoidal(f, a, b)) / 3;
+}
+
+double gauss(func_ptr f, double a, double b) {
+    double half_dx = .5 * (b - a),
+           x1 = a + half_dx * (1 - 1 / sqrt(3)),
+           x2 = a + half_dx * (1 + 1 / sqrt(3));
+
+    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;
+
+    for( i = 0; i < steps; i++)
+        surface += method(f, x = a + i * dx, x + dx);
+
+    return surface;
+}

modsim/ass2/integral.h

+#ifndef INTEGRAL_H
+#define INTEGRAL_H
+
+#include <math.h>
+#include "func_ptr.h"
+
+typedef double (*method_ptr)(func_ptr, double, double);
+
+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);
+
+#endif

modsim/ass2/newton_raphson.c

+#include "newton_raphson.h"
+
+double newton_raphson(func_ptr f, func_ptr df, double x_start,
+        double epsilon, unsigned int *steps, unsigned int max_steps) {
+    double x = x_start, last_x = x - 1;
+
+    for( *steps = 0; fabs(x - last_x) >= epsilon; (*steps)++ ) {
+        if( *steps == max_steps )
+            return NAN;
+
+        last_x = x;
+        x -= f(x) / df(x);
+    }
+
+    return x;
+}

modsim/ass2/newton_raphson.h

+#ifndef NEWTON_RAPHSON_H
+#define NEWTON_RAPHSON_H
+
+#include <math.h>
 #include "func_ptr.h"
 
 double newton_raphson(func_ptr f, func_ptr df, double x_start,
-        double epsilon, unsigned int *steps, unsigned int max_steps) {
-    double x = x_start, last_x = x - 1;
-
-    for( *steps = 0; fabs(x - last_x) >= epsilon; (*steps)++ ) {
-        if( *steps == max_steps )
-            return NAN;
-
-        last_x = x;
-        x -= f(x) / df(x);
-    }
+        double epsilon, unsigned int *steps, unsigned int max_steps);
 
-    return x;
-}
+#endif

modsim/ass2/q2.c

 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
-#include "bisec.h"
+#include "bisection.h"
 
 double func(double x) {
     return x * sin(x) - 1;
@@ -13,7 +13,7 @@ int main(int argc, char *argv[]) {
     unsigned int steps, i, begin, end;
 
     if( argc != 3 ) {
-        printf("Usage: %s BEGIN END", argv[0]);
+        printf("Usage: %s BEGIN END\n", argv[0]);
         return 1;
     }
 

modsim/ass2/q3.c

 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
-#include "bisec.h"
+#include "bisection.h"
 #include "newton_raphson.h"
 #include "regula_falsi.h"
 

modsim/ass2/q5.c

@@ -2,44 +2,10 @@
 #include <stdio.h>
 #include <math.h>
 #include "func_ptr.h"
+#include "integral.h"
 
 #define DX 1e-3
 
-typedef double (*method_ptr)(func_ptr, double, double);
-
-double rectangle(func_ptr f, double a, double b) {
-    return (b - a) * f((a + b) / 2);
-}
-
-double trapezoidal(func_ptr f, double a, double b) {
-    return 0.5 * (b - a) * (f(a) + f(b));
-}
-
-double simpson(func_ptr f, double a, double b) {
-    return (2 * rectangle(f, a, b) + trapezoidal(f, a, b)) / 3;
-}
-
-#define GAUSS_F(x) (f(x / sqrt(3) * (b - a) / 2 + (a + b) / 2))
-
-double gauss(func_ptr f, double a, double b) {
-    return (b - a) / 2 * (GAUSS_F(1) + GAUSS_F(-1));
-    //double s = (b - a) / sqrt(3),
-    //    // calculate abscissae
-    //    left  = (b + a - s) / 2.0,
-    //    right = (b + a + s) / 2.0;
-    //return (b - a) * (f(left) + f(right)) / 2.0;
-}
-
-double integral(func_ptr f, method_ptr method, double a, double b) {
-    int steps = (int)((b - a) / DX), i;
-    double x, surface = 0;
-
-    for( i = 0; i < steps; i++)
-        surface += method(f, x = a + i * DX, x + DX);
-
-    return surface;
-}
-
 double f1(double x) {
     return pow(M_E, -x);
 }
@@ -49,37 +15,33 @@ double f2(double x) {
 }
 
 #define PRINT_INTEGRAL(func, method, a, b) (printf(#func " from " #a " to " \
-    #b " using %-19s %.11f\n", #method " method:", \
-            integral(&func, &method, a, b)))
-
-#define PRINT_GAUSS(func, a, b) (printf(#func " from " #a " to " \
-    #b " using %-19s %.11f\n", "gauss method:", gauss(&func, a, b)))
+    #b " using %-19s %.11f\n", #method " method:", integral(&func, &method, a, b, DX)))
 
 int main(void) {
     PRINT_INTEGRAL(f1, rectangle, 0, 1);
     PRINT_INTEGRAL(f1, trapezoidal, 0, 1);
     PRINT_INTEGRAL(f1, simpson, 0, 1);
-    PRINT_GAUSS(f1, 0, 1);
+    PRINT_INTEGRAL(f1, gauss, 0, 1);
     puts("");
     PRINT_INTEGRAL(f2, rectangle, 0, 2);
     PRINT_INTEGRAL(f2, trapezoidal, 0, 2);
     PRINT_INTEGRAL(f2, simpson, 0, 2);
-    PRINT_GAUSS(f2, 0, 2);
+    PRINT_INTEGRAL(f2, gauss, 0, 2);
     puts("");
     PRINT_INTEGRAL(f2, rectangle, 0, 20);
     PRINT_INTEGRAL(f2, trapezoidal, 0, 20);
     PRINT_INTEGRAL(f2, simpson, 0, 20);
-    PRINT_GAUSS(f2, 0, 20);
+    PRINT_INTEGRAL(f2, gauss, 0, 20);
     puts("");
     PRINT_INTEGRAL(f2, rectangle, 0, 200);
     PRINT_INTEGRAL(f2, trapezoidal, 0, 200);
     PRINT_INTEGRAL(f2, simpson, 0, 200);
-    PRINT_GAUSS(f2, 0, 200);
+    PRINT_INTEGRAL(f2, gauss, 0, 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_GAUSS(sin, 0, 8 * M_PI);
+    PRINT_INTEGRAL(sin, gauss, 0, 8 * M_PI);
 
     return 0;
 }

modsim/ass2/q7.c

@@ -92,8 +92,10 @@ int main(int argc, char *argv[]) {
 
     generations = atoi(argv[1]);
 
-    if( argc > 2 )
-        max_age = atoi(argv[2]);
+    if( argc > 2 && (max_age = atoi(argv[2])) == 1 ) {
+        puts("MAX_AGE should be at least 2 (or 0 for infinite).");
+        return 1;
+    }
 
     sequence(generations, max_age);
 }

modsim/ass2/regula_falsi.c

+#include "regula_falsi.h"
+
+double regula_falsi(func_ptr f, double s, double t, double e,
+        unsigned int *steps, unsigned int m) {
+    int side = 0;
+    double r, fr, fs = f(s), ft = f(t);
+
+    for( *steps = 0; *steps < m; (*steps)++ ) {
+        r = (fs * t - ft * s) / (fs - ft);
+
+        if( fabs(t - s) < e * fabs(t + s) )
+            break;
+
+        fr = f(r);
+
+        if( fr * ft > 0 ) {
+            t = r;
+            ft = fr;
+
+            if( side == -1 )
+                fs /= 2;
+
+            side = -1;
+        } else if (fs * fr > 0) {
+            s = r;
+            fs = fr;
+
+            if( side == 1 )
+                ft /= 2;
+
+            side = 1;
+        } else
+            break;
+    }
+
+    return r;
+}
+

modsim/ass2/regula_falsi.h

+#ifndef REGULA_FALSI_H
+#define REGULA_FALSI_H
 
+#include <math.h>
 #include "func_ptr.h"
 
 double regula_falsi(func_ptr f, double s, double t, double e,
-        unsigned int *steps, unsigned int m) {
-    int side = 0;
-    double r, fr, fs = f(s), ft = f(t);
-
-    for( *steps = 0; *steps < m; (*steps)++ ) {
-        r = (fs * t - ft * s) / (fs - ft);
-
-        if( fabs(t - s) < e * fabs(t + s) )
-            break;
-
-        fr = f(r);
-
-        if( fr * ft > 0 ) {
-            t = r;
-            ft = fr;
-
-            if( side == -1 )
-                fs /= 2;
-
-            side = -1;
-        } else if (fs * fr > 0) {
-            s = r;
-            fs = fr;
-
-            if( side == 1 )
-                ft /= 2;
-
-            side = 1;
-        } else
-            break;
-    }
-
-    return r;
-}
+        unsigned int *steps, unsigned int m);
 
+#endif