Exam project

exam/Makefile

+CFLAGS = -O1 -std=gnu11
+CFLAGS += $(shell /usr/bin/gsl-config --cflags)
+LDLIBS += $(shell /usr/bin/gsl-config --libs)
+LDLIBS += -lm
+CC = gcc
+
+default: main
+	./main > out.txt
+
+.PHONEY: clean
+clean:
+	rm main *.o plot.png out.txt
+
+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) 
+
+

exam/README.TXT

+Since my student number is Au635140 ive made exam project 18, which consists of implementing a two dimensional integrator.
+The integrator should be able to solve integrals of type 1, where the x limits are constants and the y limits functions of x.
+Thus the integrator has the signature:
+returns integrate2D(double f(double x, double y), double a, double b, double d(double x),\
+                    double u(double x), double acc, double eps)
+Where returns contains the value and an estimate of the error.
+
+The functionality of the integrator is split into 3 functions. An initializer which is the function the user uses, a 2D driver and
+a 1D driver. The initializer takes the minimum required amount of arguments to make it user friendly.
+The 2D driver is a standard recursive adaptive integrator with the higher order weights given by the trapezium rule 
+(2/6,1/6,1/6,2/6) and the lower order by the rectangle rule (1/4,1/4,1/4,1/4). Each "function" evaluation then calls the 1D driver,
+which for a given x integrates f(x,y) from d(x) to u(x). The 1D driver is of course also a recursive adaptive integrator with the
+same weights. An estimate of the error is as always given by the difference between the higher order and lower order value.
+
+For a larger project i would implement variable transformations (e.g. Clenshaw-Curtis) such that the integrator could 
+handle singular endpoint and/or infinite limits.
\ No newline at end of file

exam/functions.c

+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+#include "functions.h"
+
+//FUNCTIONS TO INTEGRATE//
+double f(double x, double y) {
+    return pow(x,2)-10+3*y;
+}
+double RB(double x, double y) {
+    return pow(1-x,2)+10*pow(y-pow(x,2),2);
+}
+double d(double x) {
+    return -sqrt(x);
+}
+double u(double x) {
+    return x;
+}
+//VARIABLE TRANSFORMATIONS//
+double x_t(double t, double a, double b) {
+    return a+(b-a)*t;
+}
+double w(double w, double a, double b) {
+    return (b-a)*w;
+}
+
+//RECURSIVE INTEGRATION FUNCTION 1D//
+double driver1D(double f(double, double), double x, double a, double b, double acc, double eps,\
+              double two, double three, double* high, double* low, double* xs) {
+    double one = f(x,x_t(xs[0],a,b)), four = f(x,x_t(xs[3],a,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)) {
+            return Q;
+    }
+    else {
+            return driver1D(f,x,a,(a+b)/2,acc/sqrt(2.0),eps,one,two,high,low,xs)+
+                   driver1D(f,x,(a+b)/2,b,acc/sqrt(2.0),eps,three,four,high,low,xs);
+    }
+}
+
+//RECURSIVE INTEGRATION FUNCTION 2D//
+double driver2D(double f(double, double), double a, double b, double acc, double eps,\
+              double two, double three, double* high, double* low, double* xs, double* error) {
+    double x1 = x_t(xs[0],a,b), a1 = d(x1), b1 = u(x1);
+    double x4 = x_t(xs[3],a,b), a4 = d(x4), b4 = u(x4);
+    double one = driver1D(f,x1,a1,b1,acc,eps,f(x1,x_t(xs[1],a1,b1)),f(x1,x_t(xs[2],a1,b1)),\
+                        high,low,xs);
+    double four = driver1D(f,x4,a4,b4,acc,eps,f(x4,x_t(xs[1],a4,b4)),f(x4,x_t(xs[2],a4,b4)),\
+                        high,low,xs);
+    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 += err;
+            return Q;
+    }
+    else {
+            return driver2D(f,a,(a+b)/2,acc/sqrt(2.0),eps,one,two,high,low,xs,error)+
+                   driver2D(f,(a+b)/2,b,acc/sqrt(2.0),eps,three,four,high,low,xs,error);
+    }
+}
+
+//INITIALIZING FUNCTION 2D//
+returns integrate2D(double f(double x, double y), double a, double b, double d(double x),\
+                    double u(double x), double acc, double eps) {
+    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;
+    double x2 = x_t(interval_points[1],a,b), a2 = d(x2), b2 = u(x2);
+    double x3 = x_t(interval_points[2],a,b), a3 = d(x3), b3 = u(x3);
+    double start2 = driver1D(f,x2,a2,b2,acc,eps,f(x2,x_t(interval_points[1],a2,b2)),f(x2,x_t(interval_points[2],a2,b2)),\
+                           Higher_order,Lower_order,interval_points);
+    double start3 = driver1D(f,x3,a3,b3,acc,eps,f(x3,x_t(interval_points[1],a3,b3)),f(x3,x_t(interval_points[2],a3,b3)),\
+                           Higher_order,Lower_order,interval_points);
+    estimate = driver2D(f,a,b,acc,eps,start2,start3,Higher_order,Lower_order,interval_points,&error);
+    returns returns = {estimate,sqrt(error)};
+    return returns;
+}
\ No newline at end of file

exam/functions.h

+#ifndef FUNCTIONS_H
+#define FUNCTIONS_H
+
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+typedef struct {double value; double error;} returns;
+
+double f(double x, double y);
+double RB(double x, double y);
+double d(double x);
+double u(double x);
+double x(double t, double a, double b);
+double w(double w, double a, double b);
+
+double driver1D(double f(double, double), double x, double a, double b, double acc, double eps,\
+              double two, double three, double* high, double* low, double* xs);
+double driver2D(double f(double, double), double a, double b, double acc, double eps,\
+              double two, double three, double* high, double* low, double* xs, double* error);
+returns integrate2D(double f(double x, double y), double a, double b, double d(double x),\
+                    double u(double x), double acc, double eps);
+
+#endif
\ No newline at end of file

exam/main.c

+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+#include "functions.h"
+
+int main() {
+    printf("Calculating an easy integral:\n");
+    returns returns = integrate2D(f,0,10.0,d,u,1.0e-8,1.0e-8);
+    printf("integrate (integrate pow(x,2)-10+3*y from -sqrt(x) to x) from 0 to 10: %10g with global error: %2g \n",returns.value,returns.error);
+    printf("Value calculated by CAS: 3117.6894\n\n");
+    printf("Calculating the integral over the Rosenbrock function:\n");
+    returns = integrate2D(RB,0,1.0,d,u,1.0e-8,1.0e-8);
+    printf("integrate (integrate pow(1-x,2)+10*pow(y-pow(x,2),2) from -sqrt(x) to x) from 0 to 1: %10g with global error: %2g \n",returns.value,returns.error);
+    printf("Value calculated by CAS: 6.3872\n");
+}

exam/out.txt

+Calculating an easy integral:
+integrate (integrate pow(x,2)-10+3*y from -sqrt(x) to x) from 0 to 10:    3117.69 with global error: 0.00393203 
+Value calculated by CAS: 3117.6894
+
+Calculating the integral over the Rosenbrock function:
+integrate (integrate pow(1-x,2)+10*pow(y-pow(x,2),2) from -sqrt(x) to x) from 0 to 1:    6.38723 with global error: 0.000472214 
+Value calculated by CAS: 6.3872