diff --git a/homework/root-finding/Makefile b/homework/root-finding/Makefile
new file mode 100755
index 0000000000000000000000000000000000000000..5eeaabcb69948b499c819c920dc8deb55fe0f0b5
--- /dev/null
+++ b/homework/root-finding/Makefile
@@ -0,0 +1,70 @@
+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
diff --git a/homework/root-finding/func.png b/homework/root-finding/func.png
new file mode 100644
index 0000000000000000000000000000000000000000..c68e21e16175177fd8c646b0638142dee73cabf2
Binary files /dev/null and b/homework/root-finding/func.png differ
diff --git a/homework/root-finding/functions.c b/homework/root-finding/functions.c
new file mode 100755
index 0000000000000000000000000000000000000000..e836089f4902a00b845335ab6486cab7b8638d29
--- /dev/null
+++ b/homework/root-finding/functions.c
@@ -0,0 +1,249 @@
+#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
diff --git a/homework/root-finding/functions.h b/homework/root-finding/functions.h
new file mode 100755
index 0000000000000000000000000000000000000000..80358492f4b395d3727a5c529e098aada1d76237
--- /dev/null
+++ b/homework/root-finding/functions.h
@@ -0,0 +1,33 @@
+#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
diff --git a/homework/root-finding/improved_wave.png b/homework/root-finding/improved_wave.png
new file mode 100644
index 0000000000000000000000000000000000000000..3c8fb2f4557dc16cc650cf0fb2dab8d7aec9abd3
Binary files /dev/null and b/homework/root-finding/improved_wave.png differ
diff --git a/homework/root-finding/main.c b/homework/root-finding/main.c
new file mode 100755
index 0000000000000000000000000000000000000000..456dd1e8725ad69b4b4fe481d40820c1ba855aa0
--- /dev/null
+++ b/homework/root-finding/main.c
@@ -0,0 +1,75 @@
+#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
diff --git a/homework/root-finding/out.txt b/homework/root-finding/out.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fb3a8a455ece5e03fc6a4fcf19f6bdda0e582738
--- /dev/null
+++ b/homework/root-finding/out.txt
@@ -0,0 +1,58 @@
+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
diff --git a/homework/root-finding/simple_wave.png b/homework/root-finding/simple_wave.png
new file mode 100644
index 0000000000000000000000000000000000000000..d032b1aa7191d94db1af02ef66ff14a135fdb74e
Binary files /dev/null and b/homework/root-finding/simple_wave.png differ