ODE homework

homework/ODE/Makefile

+CFLAGS = -Wall -O3 -std=gnu11 $(shell gsl-config --cflags)
+LDLIBS = -lm $(shell gsl-config --libs)
+
+Three_body.png: path.txt Makefile
+	echo '\
+	set terminal png;\
+	set output "$@";\
+	set key top left;\
+	set tics out;\
+	set xlabel "Time (days)";\
+	set ylabel "Nr. of people";\
+	set title "Basic SIR-model";\
+	set yrange [-1:1];\
+	plot\
+	"$<" using 2:3 with line title "Body 1",\
+	"$<" using 4:5 with line title "Body 2",\
+	"$<" using 6:7 with line title "Body 3",\
+	' | gnuplot
+
+path.txt: main
+	./$<
+
+.PHONEY: clean
+clean:
+	$(RM) main Three_body.png path.txt
\ No newline at end of file

homework/ODE/main.c

+#include<math.h>
+#include<stdio.h>
+#include<assert.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_matrix.h>
+#include<gsl/gsl_blas.h>
+
+void rkstep12(void f(double t,gsl_vector* y,gsl_vector* dydt),double t,gsl_vector* yt,double h,gsl_vector* yh,gsl_vector* err) {
+    gsl_vector* k_zero = gsl_vector_alloc(yt->size);
+    f(t,yt,k_zero);
+    gsl_vector_memcpy(err,k_zero);
+    gsl_vector_scale(err,-1.0);
+    gsl_vector_scale(k_zero,h/2);
+    gsl_vector_add(k_zero,yt);
+    f(t+h/2,k_zero,yh);
+    gsl_vector_add(err,yh);
+    gsl_vector_scale(err,h);
+    gsl_vector_scale(yh,h);
+    gsl_vector_add(yh,yt);
+    gsl_vector_free(k_zero);
+}
+
+int driver(void f(double,gsl_vector*,gsl_vector*),double a,gsl_vector* ya,double b,gsl_vector* yb,double h,double acc,double eps) {
+    FILE* path_stream = fopen("path.txt","w");
+    fprintf(path_stream,"%g ",a);
+    for(int i=0;i<ya->size;i++) {fprintf(path_stream,"%g ",gsl_vector_get(ya,i));} fprintf(path_stream,"\n");
+    gsl_vector_memcpy(yb,ya);
+    int n = 1;
+    double x = a;
+    double step_size = h;
+    while(x<b) {
+        gsl_vector* yx = gsl_vector_alloc(ya->size);
+        gsl_vector* err = gsl_vector_alloc(ya->size);
+        rkstep12(f,x,yb,step_size,yx,err);
+        double loc_tol = (eps*gsl_blas_dnrm2(yb)+acc)*sqrt(step_size/(b-a));
+        double loc_err = gsl_blas_dnrm2(err);
+        if(loc_tol>loc_err) {
+            n++;
+            x += step_size;
+            gsl_vector_memcpy(yb,yx);
+            fprintf(path_stream,"%g ",x);
+            for(int i=0;i<ya->size;i++) {fprintf(path_stream,"%g ",gsl_vector_get(yb,i));} fprintf(path_stream,"\n");
+        }
+        step_size *= pow(loc_tol/loc_err,0.25)*0.95;
+        gsl_vector_free(yx);
+        gsl_vector_free(err);
+    }
+    fclose(path_stream);
+return n;
+}
+
+void load_path(gsl_matrix* path_matrix) {
+    FILE* path_stream = fopen("path.txt","r");
+    for(int i=0;i<path_matrix->size1;i++) {
+        for(int j=0;j<path_matrix->size2;j++) {
+            double x;
+            int i=fscanf(path_stream,"%lg",&x);
+            gsl_matrix_set(path_matrix,i,j,x);
+        }
+    }
+    fclose(path_stream);
+}
+
+void func(double t,gsl_vector* y, gsl_vector* dydt) {
+    double x1 = gsl_vector_get(y,0);
+    double y1 = gsl_vector_get(y,1);
+    double x2 = gsl_vector_get(y,2);
+    double y2 = gsl_vector_get(y,3);
+    double x3 = gsl_vector_get(y,4);
+    double y3 = gsl_vector_get(y,5);
+    double dist_12 = sqrt(pow(x2-x1,2)+pow(y2-y1,2));
+    double dist_13 = sqrt(pow(x3-x1,2)+pow(y3-y1,2));
+    double dist_23 = sqrt(pow(x3-x2,2)+pow(y3-y2,2));
+    gsl_vector_set(dydt,0,gsl_vector_get(y,6));
+    gsl_vector_set(dydt,1,gsl_vector_get(y,7));
+    gsl_vector_set(dydt,2,gsl_vector_get(y,8));
+    gsl_vector_set(dydt,3,gsl_vector_get(y,9));
+    gsl_vector_set(dydt,4,gsl_vector_get(y,10));
+    gsl_vector_set(dydt,5,gsl_vector_get(y,11));
+    gsl_vector_set(dydt,6,(x2-x1)/pow(dist_12,3)+(x3-x1)/pow(dist_13,3));
+    gsl_vector_set(dydt,7,(y2-y1)/pow(dist_12,3)+(y3-y1)/pow(dist_13,3));
+    gsl_vector_set(dydt,8,(x1-x2)/pow(dist_12,3)+(x3-x2)/pow(dist_23,3));
+    gsl_vector_set(dydt,9,(y1-y2)/pow(dist_12,3)+(y3-y2)/pow(dist_23,3));
+    gsl_vector_set(dydt,10,(x1-x3)/pow(dist_13,3)+(x2-x3)/pow(dist_23,3));
+    gsl_vector_set(dydt,11,(y1-y3)/pow(dist_13,3)+(y2-y3)/pow(dist_23,3));
+}
+
+int main() {
+    gsl_vector* ya = gsl_vector_alloc(12);
+    gsl_vector* yb = gsl_vector_alloc(12);
+    gsl_vector_set(ya,0,0.9700436);
+    gsl_vector_set(ya,1,-0.24308753);
+    gsl_vector_set(ya,2,-0.9700436);
+    gsl_vector_set(ya,3,0.24308753);
+    gsl_vector_set(ya,4,0);
+    gsl_vector_set(ya,5,0);
+    gsl_vector_set(ya,6,0.466203685);
+    gsl_vector_set(ya,7,0.43236573);
+    gsl_vector_set(ya,8,0.466203685);
+    gsl_vector_set(ya,9,0.43236573);
+    gsl_vector_set(ya,10,-2*0.466203685);
+    gsl_vector_set(ya,11,-2*0.43236573);
+
+    int steps = driver(func,0,ya,2,yb,1e-4,1e-6,1e-6);
+    gsl_matrix* path_matrix = gsl_matrix_alloc(steps,13);
+    load_path(path_matrix);
+
+    gsl_vector_free(ya);
+    gsl_vector_free(yb);
+return 0;
+}
\ No newline at end of file