Added regula falsi method and two-point gauss.

modsim/ass2/q5.c

@@ -19,12 +19,15 @@ 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) {
-    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;
+    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) {
@@ -48,31 +51,35 @@ 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)))
 
 int main(void) {
-    //PRINT_INTEGRAL(f1, rectangle, 0, 1);
-    //PRINT_INTEGRAL(f1, trapezoidal, 0, 1);
-    //PRINT_INTEGRAL(f1, simpson, 0, 1);
-    //puts("");
-    //PRINT_INTEGRAL(f2, rectangle, 0, 2);
-    //PRINT_INTEGRAL(f2, trapezoidal, 0, 2);
-    //PRINT_INTEGRAL(f2, simpson, 0, 2);
-    //puts("");
-    //PRINT_INTEGRAL(f2, rectangle, 0, 20);
-    //PRINT_INTEGRAL(f2, trapezoidal, 0, 20);
-    //PRINT_INTEGRAL(f2, simpson, 0, 20);
-    //puts("");
-    //PRINT_INTEGRAL(f2, rectangle, 0, 200);
-    //PRINT_INTEGRAL(f2, trapezoidal, 0, 200);
-    //PRINT_INTEGRAL(f2, simpson, 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_INTEGRAL(f1, rectangle, 0, 1);
+    PRINT_INTEGRAL(f1, trapezoidal, 0, 1);
+    PRINT_INTEGRAL(f1, simpson, 0, 1);
     PRINT_GAUSS(f1, 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);
+    puts("");
+    PRINT_INTEGRAL(f2, rectangle, 0, 20);
+    PRINT_INTEGRAL(f2, trapezoidal, 0, 20);
+    PRINT_INTEGRAL(f2, simpson, 0, 20);
+    PRINT_GAUSS(f2, 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);
+    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);
 
     return 0;
 }

modsim/ass2/regula_falsi.h

+
+#include "func_ptr.h"
+
+double regula_falsi(func_ptr f, double s, double t, double e,
+        int *steps, int m) {
+    int n, side = 0;
+    double r, fr, fs = f(s), ft = f(t);
+
+    for( n = 0; n < m; n++ ) {
+        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;
+    }
+
+    *steps = n;
+
+    return r;
+}
+