root-finding homework

homework/root-finding/Makefile

+CFLAGS = -Wall -O3 -std=gnu11 $(shell gsl-config --cflags)
+LDLIBS = -lm $(shell gsl-config --libs)
+
+default: simple_wave.png improved_wave.png func.png
+
+simple_wave.png: paths exact.txt Makefile
+	echo '\
+	set terminal png;\
+	set output "$@";\
+	set key top right;\
+	set tics out;\
+	set yrange [0:0.4];\
+	set grid;\
+	set xlabel "r (Bohr radius)";\
+	set ylabel "f(r)";\
+	set title "Shooting method with simple boundary condition";\
+	plot\
+	"simplepath2.txt" using 1:2 with line title "rmax = 2",\
+	"simplepath4.txt" using 1:2 with line title "rmax = 4",\
+	"simplepath6.txt" using 1:2 with line title "rmax = 6",\
+	"simplepath8.txt" using 1:2 with line title "rmax = 8",\
+	"simplepath10.txt" using 1:2 with line title "rmax = 10",\
+	"exact.txt" using 1:2 with line title "exact solution",\
+	' | gnuplot
+
+improved_wave.png: paths exact.txt Makefile
+	echo '\
+	set terminal png;\
+	set output "$@";\
+	set key top right;\
+	set tics out;\
+	set yrange [0:0.4];\
+	set grid;\
+	set xlabel "r (Bohr radius)";\
+	set ylabel "f(r)";\
+	set title "Shooting method with improved boundary condition";\
+	plot\
+	"improvedpath1.txt" using 1:2 with line title "rmax = 1",\
+	"improvedpath2.txt" using 1:2 with line title "rmax = 2",\
+	"improvedpath3.txt" using 1:2 with line title "rmax = 3",\
+	"improvedpath4.txt" using 1:2 with line title "rmax = 4",\
+	"exact.txt" using 1:2 with line title "exact solution",\
+	' | gnuplot
+
+func.png: Fe.txt Makefile
+	echo '\
+	set terminal png;\
+	set output "$@";\
+	set key top right;\
+	set tics out;\
+	set xlabel "Energy";\
+	set ylabel "Fe(rmax)";\
+	set title "Zoom on Fe to show instability of ODE";\
+	set xrange [-0.35:-0.32];\
+	set yrange [-4.6:-4.2];\
+	plot\
+	"Fe.txt" using 1:2 with line title "rmax=8",\
+	' | gnuplot
+
+paths exact.txt Fe.txt: main
+	./$< > out.txt
+	$(RM) N
+
+main: main.c functions.c
+
+
+.PHONEY: clean
+
+clean:
+	$(RM) main *.o *.txt *.png
\ No newline at end of file

homework/root-finding/functions.c

+#include "functions.h"
+
+static double E;
+void f(gsl_vector* x, gsl_vector* fx) {
+    double x0 =gsl_vector_get(x,0);
+    double x1 =gsl_vector_get(x,1);
+    gsl_vector_set(fx,0,400*pow(x0,3)-400*x0*x1+2*x0-2);
+    gsl_vector_set(fx,1,200*(x1-pow(x0,2)));
+}
+
+void radial(double r, gsl_vector* y, gsl_vector* dydr) {
+    double y0 = gsl_vector_get(y,0);
+    double y1 = gsl_vector_get(y,1);
+    gsl_vector_set(dydr,0,y1);
+    gsl_vector_set(dydr,1,2*(-1.0/r-E)*y0);
+}
+
+double Fe_simple(double energy, double rmax, char* save_path) {
+    double rmin = 1e-3;
+    gsl_vector* ya = gsl_vector_alloc(2);
+    gsl_vector* yb = gsl_vector_alloc(2);
+    gsl_vector_set(ya,0,rmin-rmin*rmin);
+    gsl_vector_set(ya,1,1.0-2*rmin);
+    E = energy;
+    //printf("%g\n",E);
+    int steps = driver(radial,rmin,ya,rmax,yb,0.1,1e-4,1e-4,save_path);
+    double result = gsl_vector_get(yb,0);
+    gsl_vector_free(ya);
+    gsl_vector_free(yb);
+    return result;
+}
+
+void print_vector(char s[], gsl_vector* v) {
+	printf("%s\n",s);
+	for(int i=0;i<v->size;i++) {
+		printf("%15g",gsl_vector_get(v,i));
+		printf("\n");
+	}
+}
+
+void print_matrix(gsl_matrix* M) {
+    for(int i = 0; i < M->size1; i++) {
+        for(int j = 0; j < M->size2; j++) {
+            printf("%10g\t", gsl_matrix_get(M, i, j));
+        }
+        printf("\n");
+    }
+}
+
+void GS_decomp(gsl_matrix* A, gsl_matrix* R) {
+    int n = A -> size1;
+    int m = A -> size2;
+    for(int i=0;i<m;i++) {
+        gsl_vector_view ai = gsl_matrix_column(A,i);
+        double Rii = gsl_blas_dnrm2(&(ai.vector));
+        gsl_matrix_set(R,i,i,Rii);
+        gsl_vector_scale(&(ai.vector),1/Rii);
+        for(int j=i+1;j<m;j++) {
+            gsl_vector_view aj = gsl_matrix_column(A,j);
+            double Rij;
+            gsl_blas_ddot(&(ai.vector),&(aj.vector),&Rij);
+            gsl_matrix_set(R,i,j,Rij);
+            gsl_vector* qiRij = gsl_vector_alloc(n);
+            gsl_vector_memcpy(qiRij,&(ai.vector));
+            gsl_vector_scale(qiRij,Rij);
+            gsl_vector_sub(&(aj.vector),qiRij);
+            gsl_vector_free(qiRij);
+        }
+    }
+}
+
+void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x) {
+    assert(Q->size1==Q->size2);
+    int m = Q -> size2;
+    gsl_blas_dgemv(CblasTrans,1,Q,b,0,x);
+    for(int i=m-1;i>=0;i--) {
+        double xi = gsl_vector_get(x,i);
+        for(int j=i+1;j<m;j++) {
+            xi += -gsl_matrix_get(R,i,j)*gsl_vector_get(x,j);
+        }
+        xi /= gsl_matrix_get(R,i,i);
+        gsl_vector_set(x,i,xi);
+    }
+}
+
+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, char* save_path) {
+    FILE* path_stream = fopen(save_path,"a");
+    if(save_path!="N") {
+        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);
+            if(save_path!="N") {
+                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;
+}
+
+int search(int n, double* x, double z){
+	assert(x[0]<=z && z<=x[n-1]);
+	int i=0, j=n-1;
+	while(j-i>1){
+		int mid=(i+j)/2;
+		if(z>x[mid]) i=mid; else j=mid;
+	}
+return i;
+}
+
+void cubic_interpol(int n, double x[], double y[], cubic_coefs* a) {
+	gsl_matrix* A = gsl_matrix_alloc(n,n);
+	gsl_vector* b = gsl_vector_alloc(n);
+	gsl_vector* B = gsl_vector_alloc(n);
+	gsl_vector_set(B,0,(3*(y[1]-y[0])/(x[1]-x[0])));
+	gsl_vector_set(B,n-1,(3*(y[n-1]-y[n-2])/(x[n-1]-x[n-2])));
+	for(int i=0;i<n-2;i++) {
+		double B_i1 = 3*((y[i+1]-y[i])/(x[i+1]-x[i])+(y[i+2]-y[i+1])/(x[i+2]-x[i+1])*(x[i+1]-x[i])/(x[i+2]-x[i+1]));
+		gsl_vector_set(B,i+1,B_i1);
+	}
+	for(int i=0;i<n;i++) {
+		for(int j=0;j<n;j++) {
+			if(abs(i-j)>1) {
+				gsl_matrix_set(A,i,j,0);
+			}
+			if((i-j)==1) {
+				gsl_matrix_set(A,i,j,1);
+			}
+			if((i-j==-1)) {
+				if(i==0) {gsl_matrix_set(A,i,j,1);}
+				else {
+					double Qi = (x[i]-x[i-1])/(x[i+1]-x[i]);
+					gsl_matrix_set(A,i,j,Qi);
+				}
+			}
+			if((i-j)==0) {
+				if(i==0 || i==n-1) {gsl_matrix_set(A,i,j,2);}
+				else {
+					double Di = 2*((x[i]-x[i-1])/(x[i+1]-x[i]))+2;
+					gsl_matrix_set(A,i,j,Di);
+				}
+			}
+		}
+	}
+	gsl_linalg_HH_solve(A,B,b);
+	for(int i=0;i<n-1;i++) {
+		double pi = (y[i+1]-y[i])/(x[i+1]-x[i]);
+		double bi = gsl_vector_get(b,i);
+		double bi1 = gsl_vector_get(b,i+1);
+		(*a).b[i] = bi;
+		(*a).c[i] = (3*pi-2*bi-bi1)/(x[i+1]-x[i]);
+		(*a).d[i] = (bi+bi1-2*pi)/pow((x[i+1]-x[i]),2);
+	}
+	gsl_matrix_free(A);
+	gsl_vector_free(b);
+	gsl_vector_free(B);
+}
+
+double cubic_eval(cubic_coefs coefs, int n, double* x, double* y, double z) {
+	int i = search(n,x,z);
+return y[i]+(coefs.b)[i]*(z-x[i])+(coefs.c)[i]*pow((z-x[i]),2)+(coefs.d)[i]*pow((z-x[i]),3);
+}
+
+void newton_root(void f(gsl_vector* x, gsl_vector* fx), gsl_vector* x, double eps) {
+    int dim = x->size;
+    gsl_matrix* J = gsl_matrix_alloc(dim,dim);
+    gsl_matrix* JR = gsl_matrix_alloc(dim,dim);
+    gsl_vector* fx = gsl_vector_alloc(dim);
+    gsl_vector* fdx = gsl_vector_alloc(dim);
+    gsl_vector* Dx = gsl_vector_alloc(dim);
+    f(x,fx);
+    //printf("Dim:%i\n",dim);
+    //printf("x: %g\n",gsl_vector_get(x,0));
+    //printf("f(x): %g\n",gsl_vector_get(fx,0));
+    double stepsize = 1;
+    while(gsl_blas_dnrm2(fx) > eps && stepsize > DX) {
+        for(int j=0;j<dim;j++) {
+            double* xk = gsl_vector_ptr(x,j);
+            for(int i=0;i<dim;i++) {
+                f(x,Dx);
+                double fi = gsl_vector_get(Dx,i);
+                //printf("%.10g\n",*xk);
+                //printf("fi: %.10g\n",fi);
+                *xk += DX; f(x,Dx);
+                double dfi = gsl_vector_get(Dx,i);
+                //printf("%.10g\n",*xk);
+                //printf("dfi: %.10g\n",dfi);
+                *xk -=DX;
+                gsl_matrix_set(J,i,j,(dfi-fi)/DX);
+            }
+        }
+        //printf("Derivative: %g\n",gsl_matrix_get(J,0,0));
+        gsl_vector_scale(fx,-1.0);
+        GS_decomp(J,JR);
+        GS_solve(J,JR,fx,Dx);
+        gsl_vector_scale(fx,-1.0);
+        double l = 1;
+        gsl_blas_daxpy(l,Dx,x);
+        f(x,fdx);
+        while(gsl_blas_dnrm2(fdx) > (1-l)*gsl_blas_dnrm2(fx) && l > (1.0/64.0)) {
+            l /= 2;
+            gsl_blas_daxpy(-l,Dx,x);
+            f(x,fdx);
+        }
+        f(x,fx);
+        //printf("x: %g\n",gsl_vector_get(x,0));
+        //printf("f(x): %g\n",gsl_vector_get(fx,0));
+        stepsize = gsl_blas_dnrm2(Dx)*l;
+    }
+    gsl_matrix_free(J);
+    gsl_matrix_free(JR);
+    gsl_vector_free(fx);
+    gsl_vector_free(fdx);
+    gsl_vector_free(Dx);
+}
\ No newline at end of file

homework/root-finding/functions.h

+#ifndef FUNCTIONS_H
+#define FUNCTIONS_H
+
+#include<math.h>
+#include<time.h>
+#include<stdio.h>
+#include<float.h>
+#include<assert.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_matrix.h>
+#include<gsl/gsl_blas.h>
+#include<gsl/gsl_linalg.h>
+#define RND (double)rand()/10000
+#define DX sqrt(DBL_EPSILON)
+typedef struct {double* b; double* c; double* d;} cubic_coefs;
+
+void f(gsl_vector* x, gsl_vector* fx);
+void radial(double r, gsl_vector* y, gsl_vector* dydr);
+double Fe_simple(double energy, double rmax, char* save_path);
+void aux(gsl_vector* x, gsl_vector* fx);
+void print_vector(char s[], gsl_vector* v);
+void print_matrix(gsl_matrix* M);
+void GS_decomp(gsl_matrix* A, gsl_matrix* R);
+void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x);
+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);
+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,char* save_path);
+int search(int n, double* x, double z);
+void cubic_interpol(int n, double x[], double y[], cubic_coefs* a);
+double cubic_eval(cubic_coefs coefs, int n, double* x, double* y, double z);
+void newton_root(void f(gsl_vector* x, gsl_vector* fx), gsl_vector* x, double eps);
+
+#endif
\ No newline at end of file

homework/root-finding/main.c

+#include "functions.h"
+#define nr 2*64
+
+static double xs[nr]; static double ys[nr];
+static double my_cubic_b[nr-1];
+static double my_cubic_c[nr-1];
+static double my_cubic_d[nr-1];
+static cubic_coefs coefs = {.b=my_cubic_b,.c=my_cubic_c,.d=my_cubic_d};
+
+void aux(gsl_vector* x, gsl_vector* fx) {
+    double energy = gsl_vector_get(x,0);
+    double M = cubic_eval(coefs,nr,xs,ys,energy);
+    gsl_vector_set(fx,0,M);
+}
+
+int main() {
+    printf("Finding root of Rosenbrock's valley function:\n");
+    gsl_vector* x0 = gsl_vector_alloc(2);
+    gsl_vector_set(x0,0,-3);
+    gsl_vector_set(x0,1,3);
+    print_vector("Starting guess:",x0);
+    newton_root(f,x0,1e-9);
+    print_vector("Found root:",x0);
+    gsl_vector_free(x0);
+
+    printf("\n\nFinding bound states of hydrogen with simple boundary:\n");
+    //Making spline of Fe(rmax)(e) since its instability makes newtons method inable to calculate proper derivates (See func.png)
+    for(int rmax=2;rmax<=10;rmax+=2) {
+        int index = 0;
+        for(double e=-2;e<0;e+=2.0/(nr)) {xs[index]=e;ys[index]=Fe_simple(e,rmax,"N");index+=1;}
+        cubic_interpol(nr,xs,ys,&coefs);
+        gsl_vector* x00 = gsl_vector_alloc(1);
+        double guess = -1.0;
+        gsl_vector_set(x00,0,guess);
+        printf("Rmax = %i\n",rmax);
+        print_vector("Starting guess:",x00);
+        newton_root(aux,x00,1e-9);
+        print_vector("Found root:",x00);
+        double energy = gsl_vector_get(x00,0);
+        char filename[30];
+        sprintf(filename,"simplepath%i.txt",rmax);
+        Fe_simple(energy,rmax,filename);
+        gsl_vector_free(x00);
+    }
+    printf("\n\nFinding bound states of hydrogen with improved boundary:\n");
+    for(int rmax=1;rmax<=4;rmax++) {
+        int index = 0;
+        double rd = (double)rmax;
+        for(double e=-2;e<0;e+=2.0/(nr)) {xs[index]=e;ys[index]=Fe_simple(e,rmax,"N")-rd*exp(-rd*sqrt(-2*e));index+=1;}
+        cubic_interpol(nr,xs,ys,&coefs);
+        gsl_vector* x00 = gsl_vector_alloc(1);
+        double guess = -1.0;
+        gsl_vector_set(x00,0,guess);
+        printf("Rmax = %i\n",rmax);
+        print_vector("Starting guess:",x00);
+        newton_root(aux,x00,1e-9);
+        print_vector("Found root:",x00);
+        double energy = gsl_vector_get(x00,0);
+        char filename[30];
+        sprintf(filename,"improvedpath%i.txt",rmax);
+        Fe_simple(energy,rmax,filename);
+        gsl_vector_free(x00);
+    }
+    FILE* exact_stream = fopen("exact.txt","w");
+    for(double r=1e-3;r<8;r+=1.0/64) {
+        fprintf(exact_stream,"%g %g\n",r,r*exp(-r));
+    }
+    fclose(exact_stream);
+    FILE* plot_stream = fopen("Fe.txt","w");
+    for(double r=-0.35;r<-0.32;r+=1.0e-4) {
+        fprintf(plot_stream,"%g %g\n",r,Fe_simple(r,8,"N"));
+    }
+    fclose(plot_stream);
+return 0;
+}
\ No newline at end of file

homework/root-finding/out.txt

+Finding root of Rosenbrock's valley function:
+Starting guess:
+             -3
+              3
+Found root:
+              1
+              1
+
+
+Finding bound states of hydrogen with simple boundary:
+Rmax = 2
+Starting guess:
+             -1
+Found root:
+      -0.127929
+Rmax = 4
+Starting guess:
+             -1
+Found root:
+      -0.483549
+Rmax = 6
+Starting guess:
+             -1
+Found root:
+      -0.499289
+Rmax = 8
+Starting guess:
+             -1
+Found root:
+      -0.499977
+Rmax = 10
+Starting guess:
+             -1
+Found root:
+           -0.5
+
+
+Finding bound states of hydrogen with improved boundary:
+Rmax = 1
+Starting guess:
+             -1
+Found root:
+      -0.500014
+Rmax = 2
+Starting guess:
+             -1
+Found root:
+      -0.500873
+Rmax = 3
+Starting guess:
+             -1
+Found root:
+      -0.500605
+Rmax = 4
+Starting guess:
+             -1
+Found root:
+      -0.500018