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

modsim_sander/ass2/Makefile

+CFLAGS := -Wextra -Wall -Werror -g -std=c99 -D_GNU_SOURCE
+LDFLAGS := -lm
+
+progs := bisection numdiff sqrt newton_raphson
+
+all: $(progs)
+
+bisection: bisection.c
+numdiff: numdiff.c
+sqrt: sqrt.c
+newton_raphson: newton_raphson.c
+
+$(progs): CFLAGS += -DMAIN_$@
+
+clean: ; $(RM) $(progs)

modsim_sander/ass2/bisection.c

+#include <stdio.h>
+#include <math.h>
+
+
+static double bisect(double (*f)(double), double a, double b, double tolerance,
+        unsigned int *steps, unsigned int max_steps, double approximation,
+        int error_per_step)
+{
+    double c = INFINITY;
+
+    for ( *steps = 1; *steps < max_steps; *steps += 1) {
+        c = (a + b) / 2;
+
+        if ( error_per_step )
+            printf("step %2d error: %.3E (%3f%%)\n", *steps,
+                fabs(approximation - c),
+                100 * fabs(c - approximation) / approximation);
+
+        if (f(c) == 0 || (b - a) / 2 < tolerance)
+            return c;
+
+        if (signbit(f(c)) == signbit(f(a)))
+            a = c;
+        else
+            b = c;
+    }
+
+    return c;
+}
+
+
+#ifdef MAIN_bisection
+
+static inline double f(double x) { return x * sin(x) - 1.0; }
+
+int main(void)
+{
+    double approximation = 1.11415714087193008730052518;
+    unsigned int steps;
+
+    bisect(&f, 0.0, 2.0, 10e-12, &steps, 100, approximation, 0);
+
+    printf("bisection of x * sin(x) - 1 for x in [%.1f, %.1f] in %d steps.\n",
+           0.0, 2.0, steps);
+
+    return 0;
+}
+
+#endif

modsim_sander/ass2/newton_raphson.c

+#include <stdio.h>
+#include <math.h>
+
+
+static double newton_raphson(double (*f)(double), double (*df)(double),
+        double start_x, double tolerance, unsigned int *steps,
+        unsigned int max_steps, double approximation, int error_per_step)
+{
+    double x = start_x, last_x = x - 1;
+
+    for( *steps = 1; fabs(x - last_x) >= tolerance; (*steps)++ ) {
+        if ( *steps == max_steps )
+            return NAN;
+
+        last_x = x;
+        x -= f(x) / df(x);
+
+        if ( error_per_step )
+            printf("step %2d error: %.3E (%3f%%)\n", *steps,
+                fabs(approximation - last_x),
+                100 * fabs(last_x - approximation) / approximation);
+    }
+
+    return x;
+}
+
+#ifdef MAIN_newton_raphson
+
+static inline double f1(double x) { return x * x - x + 2; }
+static inline double df1(double x) { return 2 * x - 1; }
+
+static inline double f2(double x) { return x * x * x - 3 * x - 2; }
+static inline double df2(double x) { return 3 * x * x - 3; }
+
+static inline double f3(double x) { return (x * x + 1) * (x - 4); }
+static inline double df3(double x) { return 3 * x * x - 8 * x + 1; }
+
+int main(void)
+{
+    unsigned int steps;
+    double result;
+
+    result = newton_raphson(&f1, &df1, 0, 10e-12, &steps, 100, 0, 0);
+    printf("f1 using newton-raphson = %.20f (%d steps)\n", result, steps);
+
+    result = newton_raphson(&f2, &df2, -10, 10e-12, &steps, 100, 0, 0);
+    printf("f2 using newton-raphson = %.20f (%d steps)\n", result, steps);
+
+    result = newton_raphson(&f2, &df2, 10, 10e-12, &steps, 100, 0, 0);
+    printf("f2 using newton-raphson = %.20f (%d steps)\n", result, steps);
+
+    result = newton_raphson(&f3, &df3, 0, 10e-12, &steps, 100, 0, 0);
+    printf("f3 using newton-raphson = %.20f (%d steps)\n", result, steps);
+
+    return 0;
+}
+
+#endif

modsim_sander/ass2/numdiff.c

+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+
+static double numdiff_central(double (*f)(double), double x, double h) {
+    return (f(x + h) - f(x - h)) / (2 * h);
+}
+
+
+static double numdiff_right(double (*f)(double), double x, double h) {
+    return (f(x + h) - f(x)) / h;
+}
+
+
+#define NUMDIFF(f, x) { \
+    volatile double h = sqrt(DBL_EPSILON) * (x); \
+    \
+    printf("right-hand diff. of " #f "(x) for x = " #x ": %lf\n", \
+           numdiff_right(f, x, h)); \
+    \
+    printf("central diff.    of " #f "(x) for x = " #x ": %lf\n", \
+           numdiff_central(f, x, h)); \
+    \
+    printf("calculated h = √e * (" #x "): %lf\n\n", h); \
+}
+
+
+int main(void)
+{
+    NUMDIFF(sin, M_PI / 3.0);
+    NUMDIFF(sin, 100 * M_PI + M_PI / 3.0);
+    NUMDIFF(sin, 10e12 * M_PI + M_PI / 3.0);
+
+    return 0;
+}

modsim_sander/ass2/sqrt.c

+#include <stdio.h>
+#include <math.h>
+
+#include "bisection.c"
+#include "newton_raphson.c"
+
+static double regula_falsi(double (*f)(double), double a, double b,
+        double tolerance, unsigned int *steps, unsigned int max_steps,
+        double approximation, int error_per_step)
+{
+    int side = 0;
+    double c, fc, fa = f(a), fb = f(b);
+
+    for( *steps = 1; *steps < max_steps; (*steps)++ ) {
+        c = (fa * b - fb * a) / (fa - fb);
+
+        if ( error_per_step )
+            printf("step %2d error: %.3E (%3f%%)\n", *steps,
+                fabs(approximation - c),
+                100 * fabs(c - approximation) / approximation);
+
+        if ( fabs(b - a) < tolerance * fabs(b + a) )
+            break;
+
+        fc = f(c);
+
+        if ( fc * fb > 0 ) {
+            b = c;
+            fb = fc;
+
+            if( side == -1 )
+                fa /= 2;
+
+            side = -1;
+        } else if (fa * fc > 0) {
+            a = c;
+            fa = fc;
+
+            if( side == 1 )
+                fb /= 2;
+
+            side = 1;
+        } else
+            break;
+    }
+
+    return c;
+}
+
+#ifdef MAIN_sqrt
+
+static inline double f(double x) { return x * x - 2; }
+
+static inline double df(double x) { return 2 * x; }
+
+int main(void)
+{
+    unsigned int steps;
+
+    double approx = 1.41421356237309504880168872420, result;
+    printf("√2 is approximate       = %.20f\n", approx);
+
+    result = bisect(&f, 1, 2, 10e-12, &steps, 100, approx, 0);
+    printf("√2 using bisection      = %.20f (%d steps)\n", result, steps);
+
+    result = regula_falsi(&f, 1, 2, 10e-12, &steps, 100, approx, 0);
+    printf("√2 using regula falsi   = %.20f (%d steps)\n", result, steps);
+
+    result = newton_raphson(&f, &df, 1.4, 10e-12, &steps, 100, approx, 0);
+    printf("√2 using newton-raphson = %.20f (%d steps)\n", result, steps);
+
+    return 0;
+}
+
+#endif