Numerical integration homework

homework/numerical-integration/Makefile

+CFLAGS = -O1 -std=gnu11
+CFLAGS += $(shell /usr/bin/gsl-config --cflags)
+LDLIBS += $(shell /usr/bin/gsl-config --libs)
+LDLIBS += -lm
+CC = gcc
+
+.PHONEY: default
+default: main run 
+.PHONEY: run
+run:
+	./main > out.txt
+
+.PHONEY: clean
+clean:
+	rm main *.o *.txt plot.png
+
+main: main.o functions.o
+
+functions.o: functions.c functions.h
+	$(CC) $(CFLAGS) -c functions.c $(LDLIBS) 
+
+main.o: main.c functions.h
+	$(CC) $(CFLAGS) -c main.c $(LDLIBS) 
+
+

homework/numerical-integration/functions.c

+#include<stdio.h>
+#include<string.h>
+#include<math.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_matrix.h>
+#include<stdlib.h>
+#include<gsl/gsl_linalg.h>
+#include<gsl/gsl_blas.h>
+#include<assert.h>
+#include "functions.h"
+
+//FUNCTIONS TO INTEGRATE//
+double f(double x)
+{
+    return 1/sqrt(x);
+}
+double g(double x)
+{
+    return 4*sqrt(1-pow(x,2));
+}
+double h(double x)
+{
+    return exp(-pow(x,2));
+}
+//VARIABLE TRANSFORMATIONS//
+double x(double t, double a, double b)
+{
+    return a+(b-a)*t;
+}
+double w(double w, double a, double b)
+{
+    return (b-a)*w;
+}
+double CC(double f(double), double x, double a, double b)
+{
+    return f((a+b)/2+(b-a)/2*cos(x))*sin(x)*(b-a)/2;
+}
+double a_inf(double f(double), double x, double b)
+{
+    return f(b-(1-x)/x)/pow(x,2);
+}
+double b_inf(double f(double), double x, double a)
+{
+    return f(a+(1-x)/x)/pow(x,2);
+}
+double both_inf(double f(double), double x)
+{
+    return (f((1-x)/x)+f(-(1-x)/x))/pow(x,2);
+}
+//RECURSIVE INTEGRATION FUNCTION//
+double driver(double f(double), double org_a, double org_b, double a, double b, double acc, double eps,\
+              double two, double three, double* high, double* low, double* xs, int* evals, double* error, char* type)
+{
+    double one, four;
+    if(type=="normal") {
+        one = f(x(xs[0],a,b));
+        four = f(x(xs[3],a,b));
+    }
+    if(type=="a_inf") {
+        one = a_inf(f,x(xs[0],a,b),org_b);
+        four = a_inf(f,x(xs[3],a,b),org_b);
+    }
+    if(type=="b_inf") {
+        one = b_inf(f,x(xs[0],a,b),org_a);
+        four = b_inf(f,x(xs[3],a,b),org_a);
+    }
+    if(type=="both_inf") {
+        one = both_inf(f,x(xs[0],a,b));
+        four = both_inf(f,x(xs[3],a,b));
+    }
+    if(type=="CC") {
+        one = CC(f,x(xs[0],a,b),org_a,org_b);
+        four = CC(f,x(xs[3],a,b),org_a,org_b);
+    }
+    double Q = w(high[0],a,b)*one+w(high[1],a,b)*two+w(high[2],a,b)*three+w(high[3],a,b)*four; 
+    double q = w(low[0],a,b)*one+w(low[1],a,b)*two+w(low[2],a,b)*three+w(low[3],a,b)*four;
+    double err=fabs(Q-q);
+    if (err < acc+eps*fabs(Q))
+        {
+            *error += pow(err,2);
+            return Q;
+        }
+    else 
+        {
+            *evals += 2;
+            return driver(f,org_a,org_b,a,(a+b)/2,acc/sqrt(2.0),eps,one,two,high,low,xs,evals,error,type)+
+                    driver(f,org_a,org_b,(a+b)/2,b,acc/sqrt(2.0),eps,three,four,high,low,xs,evals,error,type);
+        }
+}
+//INITIALIZING FUNCTION//
+returns integrate(double f(double), double a, double b, double acc, double eps, char* is_CC) 
+{
+    assert(is_CC=="CC" || is_CC=="no_CC" && "is_CC param must be CC or no_CC");
+    if(is_CC=="CC") {assert(a!=-INFINITY && b!=INFINITY && "CC not available with infinite boundaries");}
+    double Higher_order[4] = {2.0/6,1.0/6,1.0/6,2.0/6};
+    double Lower_order[4] = {1.0/4,1.0/4,1.0/4,1.0/4};
+    double interval_points[4] = {1.0/6,2.0/6,4.0/6,5.0/6};
+    double estimate = 0;
+    double error = 0;
+    int evals = 0;
+
+    if(is_CC=="no_CC" && a!=-INFINITY && b!=INFINITY) {
+        estimate = driver(f,a,b,a,b,acc,eps,f(x(interval_points[1],a,b)),f(x(interval_points[2],a,b)),\
+                          Higher_order,Lower_order,interval_points,&evals,&error,"normal");
+    }
+    if(is_CC=="no_CC" && a==-INFINITY && b!=INFINITY) {
+        estimate = driver(f,a,b,0,1,acc,eps,a_inf(f,x(interval_points[1],0,1),b),a_inf(f,x(interval_points[2],0,1),b),\
+                          Higher_order,Lower_order,interval_points,&evals,&error,"a_inf");
+    }
+    if(is_CC=="no_CC" && a!=-INFINITY && b==INFINITY) {
+        estimate = driver(f,a,b,0,1,acc,eps,b_inf(f,x(interval_points[1],0,1),a),b_inf(f,x(interval_points[2],0,1),a),\
+                          Higher_order,Lower_order,interval_points,&evals,&error,"b_inf");
+    }
+    if(is_CC=="no_CC" && a==-INFINITY && b==INFINITY) {
+        estimate = driver(f,a,b,0,1,acc,eps,both_inf(f,x(interval_points[1],0,1)),both_inf(f,x(interval_points[2],0,1)),\
+                          Higher_order,Lower_order,interval_points,&evals,&error,"both_inf");
+    }
+    if(is_CC=="CC") {
+        estimate = driver(f,a,b,0,M_PI,acc,eps,CC(f,x(interval_points[1],0,M_PI),a,b),CC(f,x(interval_points[2],0,M_PI),a,b),\
+                             Higher_order,Lower_order,interval_points,&evals,&error,"CC");
+    }
+    
+    returns returns = {estimate,sqrt(error)};
+    printf("Integration(%s) took %i evals\n",is_CC,evals);
+    return returns;
+}

homework/numerical-integration/functions.h

+#ifndef FUNCTIONS_H
+#define FUNCTIONS_H
+
+#include<stdio.h>
+#include<math.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_matrix.h>
+#include<stdlib.h>
+#include<gsl/gsl_linalg.h>
+#include<gsl/gsl_blas.h>
+#include<assert.h>
+
+typedef struct {double value; double error;} returns;
+
+double f(double x);
+double g(double x);
+double h(double x);
+double x(double t, double a, double b);
+double w(double w, double a, double b);
+double CC(double f(double), double x, double a, double b);
+double a_inf(double f(double), double x, double b);
+double b_inf(double f(double), double x, double a);
+double both_inf(double f(double), double x);
+double driver(double f(double), double org_a, double org_b, double a, double b, double acc, double eps,\
+              double two, double three, double* high, double* low, double* xs, int* evals, double* error, char* type);
+returns integrate(double f(double), double a, double b, double acc, double eps, char* is_CC);
+
+
+
+
+
+#endif
\ No newline at end of file

homework/numerical-integration/main.c

+#include<stdio.h>
+#include<math.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_matrix.h>
+#include<stdlib.h>
+#include<gsl/gsl_linalg.h>
+#include<gsl/gsl_blas.h>
+#include<assert.h>
+#include "functions.h"
+
+int main()
+{
+    printf("To investigate the effect of CC-transformation on function with singularity at endpoint:\n");
+    returns returns = integrate(f,0,1,1.0e-8,1.0e-8,"no_CC");
+    printf("1/sqrt(x) from 0 to 1 without CC: %34g    with global error: %2g \n",returns.value,returns.error);
+    printf("\n");
+    returns = integrate(f,0,1,1.0e-8,1.0e-8,"CC");
+    printf("1/sqrt(x) from 0 to 1 with CC: %37g    with global error: %2g \n",returns.value,returns.error);
+    printf("\n");printf("\n");
+
+    printf("To investigate the effect of CC-transformation on PI precision:\n");
+    returns = integrate(g,0,1,1.0e-8,1.0e-8,"no_CC");
+    printf("4*sqrt(1-pow(x,2)) from 0 to 1 without CC: %25.20g    with global error: %2g \n",returns.value,returns.error);
+    printf("\n");
+    returns = integrate(g,0,1,1.0e-8,1.0e-8,"CC");
+    printf("4*sqrt(1-pow(x,2)) from 0 to 1 with CC: %28.20g    with global error: %2g \n",returns.value,returns.error);
+    printf("\n");printf("\n");
+
+    printf("To demonstrate converging integrals with infinite boundaries:\n");
+    returns = integrate(h,-INFINITY,0,1.0e-8,1.0e-8,"no_CC");
+    printf("exp(-pow(x,2)) from -INFINITY to 0 without CC: %21g    with global error: %2g \n",returns.value,returns.error);
+    printf("\n");
+    returns = integrate(h,0,INFINITY,1.0e-8,1.0e-8,"no_CC");
+    printf("exp(-pow(x,2)) from 0 to INFINITY without CC: %22g    with global error: %2g \n",returns.value,returns.error);
+    printf("\n");
+    returns = integrate(h,-INFINITY,INFINITY,1.0e-8,1.0e-8,"no_CC");
+    printf("exp(-pow(x,2)) from -INFINITY to INFINITY without CC: %14g    with global error: %2g \n",returns.value,returns.error);
+    return 0;
+}

homework/numerical-integration/out.txt

+To investigate the effect of CC-transformation on function with singularity at endpoint:
+Integration(no_CC) took 400670 evals
+1/sqrt(x) from 0 to 1 without CC:                                  2    with global error: 5.9092e-09 
+
+Integration(CC) took 950 evals
+1/sqrt(x) from 0 to 1 with CC:                                     2    with global error: 6.81757e-09 
+
+
+To investigate the effect of CC-transformation on PI precision:
+Integration(no_CC) took 1892 evals
+4*sqrt(1-pow(x,2)) from 0 to 1 without CC:     3.1415926535960836397    with global error: 5.97213e-09 
+
+Integration(CC) took 3006 evals
+4*sqrt(1-pow(x,2)) from 0 to 1 with CC:        3.1415926535897393812    with global error: 5.26236e-09 
+
+
+To demonstrate converging integrals with infinite boundaries:
+Integration(no_CC) took 1646 evals
+exp(-pow(x,2)) from -INFINITY to 0 without CC:              0.886227    with global error: 4.88183e-09 
+
+Integration(no_CC) took 1646 evals
+exp(-pow(x,2)) from 0 to INFINITY without CC:               0.886227    with global error: 4.88183e-09 
+
+Integration(no_CC) took 2026 evals
+exp(-pow(x,2)) from -INFINITY to INFINITY without CC:        1.77245    with global error: 5.24721e-09