/* ************************************************************************ * numerov.c: Numerov method to find eigenvalues and eigenfunctions of * * a particle in a potential well. As written this program * * only finds odd eigenfunctions correctly and fails for * * even eigenfunctions and if there is no eigenvalue in the * * specified range. The final function is not normalized. * * * * taken from: "Projects in Computational Physics" by Landau and Paez * * copyrighted by John Wiley and Sons, New York * * * * written by: students in PH465/565, Computational Physics, * * at Oregon State University * * code copyrighted by RH Landau * * supported by: US National Science Foundation, Northwest Alliance * * for Computational Science and Engineering (NACSE), * * US Department of Energy * * * * UNIX (DEC OSF, IBM AIX): cc numerov.c -lm * * * ************************************************************************ */ #include #include #define steps 1000 /* half width of potential well */ #define V -0.001 /* depth of potential well */ #define eps 1E-8 /* accuracy for eigenvalues */ #define Emin -0.001 /* set to find the third */ #define Emax -0.00085 /* eigenfunction */ double k2(int i, double E); /* returns potential at x */ double diff(double E); /* difference of derivatives */ void plot(double E); /* data for final plot */ main() { double E, min, max; int i=0; /* counter for iterations */ min=Emin; max=Emax; do { i++; E=(max+min)/2.0; /* divide energy range */ if (diff(max)*diff(E)>0) max=E; /* the bisection algorithm */ else min=E; }while(fabs(diff(E))>eps); printf("eigenvalue E=%.10f\n", E); printf("after %d iterations\n", i); plot(E); } /*----------------------- end of main program ------------------------*/ /* function returns difference between left and right wavefunction */ double diff(double E) { double one, two, three, plus, minus; int i; one=0.0; two=0.00001; for (i=1; i<=1500; i++) /* left side first */ { three=(2*two*(1-5./12.*k2(i,E))- (1.+1./12.*k2(i-1,E))*one)/(1.+1./12.*k2(i+1,E)); one=two; two=three; } minus=two; /* value at matching point */ one=0.0; /* reset starting conditions */ two=0.00001; for (i=1; i<500; i++) /* now the right side */ { three=(2*two*(1.-5./12.*k2(i,E))- (1.+1./12.*k2(i+1,E))*one)/(1.+1./12.*k2(i-1,E)); one=two; two=three; } plus=two; /* value at matching point */ return((minus-plus)); } /*--------------------------------------------------------------------*/ /* function returns k-vector at depending on position i */ double k2(int i, double E) { if (i<500) return(E); /* outside the well */ if (i>=500) return (E-(V)); /* inside the well */ } /*--------------------------------------------------------------------*/ /* write data for eigenfuntion into files left.dat, right.dat */ void plot(double E) { double one, two, three; int i; FILE *right, *left; /* save data in files */ right=fopen("right.dat", "w"); left=fopen("left.dat", "w"); one=0.0; two=0.00001; for (i=1; i<=1500; i++) /* left side first */ { three=(2*two*(1-5./12.*k2(i,E))- (1+1./12.*k2(i-1,E))*one)/(1+1./12.*k2(i+1,E)); fprintf(left, "%d\t%f\n", i-1000, three/2.8); one=two; two=three; } one=0.0; /* reset starting conditions */ two=0.00001; for (i=1; i<500; i++) /* now the right side */ { three=(2*two*(1.-5./12.*k2(i,E))- (1+1./12.*k2(i+1,E))*one)/(1.+1./12.*k2(i-1,E)); fprintf(right, "%d\t%f\n", 1000-i, three/2.8); one=two; two=three; } fclose(left); fclose(right); printf("data saved in left.dat and right.dat\n"); }