Numerov as a class

CMakeLists.txt

@@ -2,6 +2,7 @@ cmake_minimum_required(VERSION 3.5)
 project(Schroedinger)
 
 # Force C++17 standard
+set(CMAKE_BUILD_TYPE             Debug)
 set(CMAKE_CXX_STANDARD              17)
 set(CMAKE_CXX_STANDARD_REQUIRED     ON)
 set(CMAKE_CXX_FLAGS -D_EMULATE_GLIBC=0)

src/Potential/Potential.cpp

@@ -1,5 +1,8 @@
 #include "Potential.h"
 
+Potential::Potential() {
+    throw std::invalid_argument("Can't initialize a potential without specifying parameters");
+}
 Potential::Potential(Base base, PotentialType type, double k, double width, double height, bool separable)
 {
     this->base      = base; 

src/Potential/Potential.h

@@ -8,7 +8,7 @@
 #include <stdexcept>
 #include <stdbool.h>
 #include <fstream>
-#include "../Basis/Base.h"
+#include <Base.h>
 
 /*! Class Potential contains the potential used in the Schroedinger equation.
  * takes the necessary input: std::vector x at definition Builder(x),
@@ -38,6 +38,7 @@ class Potential {
         FINITE_WELL_POTENTIAL = 2,   
     };
 
+    Potential();
     Potential(Base, PotentialType, double, double, double, bool);
     std::vector<double> getValues();
     std::vector<double> getCoordsFromBase();

src/Solver/Numerov.cpp

@@ -1,159 +1,160 @@
 #include "Numerov.h"
 
-/*! Integrate with the trapezoidal rule method, from a to b position in a function array*/
-double trap_array(int a, int b, double stepx, double *func) {
-    double trapez_sum = 0.;
+Numerov::Numerov(Potential potential, int nbox) {
+    this->potential      = &potential;
+    this->nbox           = nbox;
+    this->solutionEnergy = 0;
+    this->probability    = std::vector<double>(nbox+10);
+    this->wavefunction   = std::vector<double>(nbox+10);
 
-    for (int j = a + 1; j < b; j++) {
-        trapez_sum += func[j];
-    }
-    trapez_sum += (func[a] + func[b]) / 2.;
-    trapez_sum = (trapez_sum) * stepx;
-    return trapez_sum;
+    std::fill(wavefunction.begin(), wavefunction.end(), 0);
+
+    this->wavefunction.at(0) = 0;
+    this->wavefunction.at(1) = 0.1;
 }
 
-/*! Numerov Algorithm solves
-f''(x) + v(x)f(x) = 0,
-by considering
-\left( 1+ \frac{h^2}{12} v(x+h) \right) f(x+h) = 2 \left( 1 - \frac{5h^2}{12} v(x) \right) f(x) - \left( 1 + \frac{h^2}{12} v(x-h) \right) f(x-h).
-for the Shroedinger equation v(x) = V(x) - E, where V(x) is the potential and E the eigenenergy
+/*! 
+    Numerov Algorithm solves f''(x) + v(x)f(x) = 0,
+    by considering
+    \left( 1+ \frac{h^2}{12} v(x+h) \right) f(x+h) = 2 \left( 1 - \frac{5h^2}{12} v(x) \right) f(x) - \left( 1 + \frac{h^2}{12} v(x-h) \right) f(x-h).
+    for the Shroedinger equation v(x) = V(x) - E, where V(x) is the potential and E the eigenenergy
 */
-void fsol_Numerov(double Energy, int nbox, Potential V, double *wavefunction) {
+void Numerov::functionSolve(double energy) {
     double c, x;
-    std::vector<double> potential = V.getValues();
-
-    c = (2. * mass / hbar / hbar) * (dx * dx / 12.);
-    //Build Numerov f(x) solution from left.
-    for (int i = 2; i <= nbox; i++) {
-        x = (-nbox / 2 + i) * dx;
-
-        wavefunction[i] = 2 * (1. - (5 * c) * (Energy - potential[i-1])) * wavefunction[i - 1]
-                  - (1. + (c) * (Energy - potential[i-2])) * wavefunction[i - 2];
-        wavefunction[i] /= (1. + (c) * (Energy - potential[i]));
+    std::cout << "Inizio di fare fsolve";
+
+    std::vector<double> pot = potential->getValues();
+    c = (2.0 * mass / hbar / hbar) * (dx * dx / 12.0);
+
+    //Build Numerov f(x) solution from left. 
+    for (int i = 2; i <= this->nbox; i++) {
+        x = (-this->nbox / 2 + i) * dx;
+        double value = 2 * (1. - (5 * c) * (energy - pot.at(i-1) )) * wavefunction.at(i-1) - (1. + (c) * (energy - pot.at(i-2))) * wavefunction.at(i-2);
+        value /= (1.0 + (c) * (energy - pot.at(i)));    
+        this->wavefunction.insert((this->wavefunction.begin()+i), value );
+        //std::cout << i << " " << this->wavefunction.at(i) << std::endl;
     }
-
-    /* //right solution
-
-    norm = wavefunction[nbox/2];
-    wavefunction[nbox-1] = first_step;
-    for(int i=nbox-2;i>=nbox/2;i--){
-      x = (i-nbox/2)*dx;
-      wavefunction[i] = -   ( 1. + (  c)*(- Energy + (*potential)(x+2*step)) ) *  wavefunction[i+2]
-                + 2*( 1. - (5*c)*(- Energy + (*potential)(x+step))   ) * wavefunction[i+1] ;
-
-      wavefunction[i]/= ( 1. + (  c)*(- Energy + (*potential)(x)) );
-    }*/
 }
 
-/*! \brief a solver of differential equation using Numerov algorithm and selecting non-trivial solutions.
-@param (*potential) is the pointer to the potential function, takes function of 1 variable as input
-@param wavefunction, takes array of @param nbox size as input (for preconditioning)
-  and gives the output the normalizedwf solution
-
-solve_Numerov solves the pointed potential using the Numerov algorithm and
-renormalizing the output wavefunction to 1. To do this it must try the solutions
-for different energies. The natural solution to the second degree differential equation
-is the exponential. But my boundary conditions impose 0 at both beginning and end
-of the wavefunction, so you have to try until you find such solution by finding
- where the exponential solution changes sign.
+/*! 
+    \brief a solver of differential equation using Numerov algorithm and selecting non-trivial solutions.
+    @param (*potential) is the pointer to the potential function, takes function of 1 variable as input
+    @param wavefunction, takes array of @param nbox size as input (for preconditioning)
+    and gives the output the normalizedwf solution
+
+    solve_Numerov solves the pointed potential using the Numerov algorithm and
+    renormalizing the output wavefunction to 1. To do this it must try the solutions
+    for different energies. The natural solution to the second degree differential equation
+    is the exponential. But my boundary conditions impose 0 at both beginning and end
+    of the wavefunction, so you have to try until you find such solution by finding
+    where the exponential solution changes sign.
 */
-double solve_Numerov(double Emin, double Emax, double Estep,
-                   int nbox, Potential V, double *wavefunction) {
 
-    double c, x, first_step, norm, Energy, Solution_Energy = 0.;
-    int n, sign;
-
-    double *probab = new double[nbox];
+double Numerov::solve(double e_min, double e_max, double e_step) {
 
-    // scan energies to find when the Numerov solution is =0 at the right extreme of the box.
-    for (n = 0; n < (Emax - Emin) / Estep; n++) {
-        Energy = Emin + n * Estep;
-        // wavefunction[1] = first_step;
+    double c, x, first_step, norm, energy = 0.0;
+    int n, sign;
 
-        fsol_Numerov(Energy, nbox, V, wavefunction);
-        // std::coutS << "# Energy = " << Energy << "  " << wavefunction[nbox] << std::endl;
+    // scan energies to find when the Numerov solution is = 0 at the right extreme of the box.
+    for (n = 0; n < (e_max - e_min) / e_step; n++) {
+        energy = e_min + n * e_step;
 
+        this->functionSolve(energy);
 
-        if (fabs(wavefunction[nbox]) < err) {
-            std::cout << "#solution found" << wavefunction[nbox] << std::endl;
-            Solution_Energy = Energy;
+        std::cout << "\n\n Aftersolve\n\n";
+        if ( fabs(this->wavefunction.at(this->nbox)) < err ) {
+            std::cout << "Solution found" << this->wavefunction.at(this->nbox) << std::endl;
+            this->solutionEnergy = energy;
             break;
         }
 
         if (n == 0)
-            sign = (wavefunction[nbox] > 0) ? 1 : -1;
+            sign = (this->wavefunction.at(this->nbox) > 0) ? 1 : -1;
 
         // when the sign changes, means that the solution for f[nbox]=0 is in in the middle, thus calls bisection rule.
-        if (sign * wavefunction[nbox] < 0) {
-          std::cout << "#bisection " << wavefunction[nbox] << std::endl;
-          Solution_Energy = bisec_Numer(Energy - Estep, Energy + Estep, nbox, V, wavefunction);
+        if (sign * this->wavefunction.at(this->nbox) < 0) {
+            std::cout << "Bisection " << wavefunction.at(nbox) << std::endl;
+            this->solutionEnergy = this->bisection(energy - e_step, energy + e_step);
             break;
         }
     }
 
-    std::cout << "# iteration " << n << "  Energy = " << Solution_Energy << std::endl;
-
     for (int i = 0; i <= nbox; i++)
-        probab[i] = wavefunction[i] * wavefunction[i];
+        this->probability.at(i) = this->wavefunction.at(i) * this->wavefunction.at(i);
 
-    norm = trap_array(0., nbox, dx, probab);
-    std::cout << "# norm=" << norm << std::endl;
+    norm = trapezoidalRule(0.0, this->nbox, dx, this->probability);
 
-    std::ofstream myfile ("wavefunction.dat");
-    if (myfile.is_open())
-    {
-
-        for (int i = 0; i <= nbox; i++) {
-            wavefunction[i] = wavefunction[i] / sqrt(norm);
-            myfile << i << " "<< wavefunction[i] << std::endl;
-        }
-        std::cout << std::endl;
-        myfile.close();
+    for (int i = 0; i <= nbox; i++) {
+        wavefunction.at(i) = wavefunction.at(i) / sqrt(norm);
     }
 
-
-    // for (int i = 0; i <= nbox; i++)
-    //     std::cout << (-nbox / 2 + i) * dx << "  " << wavefunction[i] << " " << (*potential)((-nbox / 2 + i) * dx) << std::endl;
-
-    return Solution_Energy;
+    this->printToFile();
+    return this->solutionEnergy;
 }
 
 /*! Applies a bisection algorith to the numerov method to find
 the energy that gives the non-trivial (non-exponential) solution
 with the correct boundary conditions (@param wavefunction[0] == @param wavefunction[@param nbox] == 0)
 */
-double bisec_Numer(double Emin, double Emax, int nbox, Potential V, double *wavefunction) {
-    double Emiddle, fx1, fb, fa;
+double Numerov::bisection(double e_min, double e_max) {
+    double energy_middle, fx1, fb, fa;
     std::cout.precision(17);
 
     // The number of iterations that the bisection routine needs can be evaluated in advance
-    int itmax = ceil(log2(Emax - Emin) - log2(err)) - 1;
+    int itmax = ceil(log2(e_max - e_min) - log2(err)) - 1;
 
-    std::cout << "#itmax=" << itmax << std::endl;
     for (int i = 0; i < itmax; i++) {
-        Emiddle = (Emax + Emin) / 2.;
-        fsol_Numerov(Emiddle, nbox, V, wavefunction);
-        fx1 = wavefunction[nbox];
-        fsol_Numerov(Emax, nbox, V, wavefunction);
-        fb = wavefunction[nbox];
+        energy_middle = (e_max + e_min) / 2.0;
+
+        this->functionSolve(energy_middle);
+        fx1 = this->wavefunction.at(this->nbox);
+
+        this->functionSolve(e_max);
+        fb = this->wavefunction.at(this->nbox);
 
         if (std::abs(fx1) < err) {
-            std::cout << "#Numerov E=" << Emiddle << "f(nbox=" << nbox << ") = " << fx1 << " " << wavefunction[nbox] << std::endl;
-            return Emiddle;
+            return energy_middle;
         }
 
         if (fb * fx1 < 0.) {
-            Emin = Emiddle;
+            e_min = energy_middle;
         } else {
-            fsol_Numerov(Emin, nbox, V, wavefunction);
-            fa = wavefunction[nbox];
+            this->functionSolve(e_min);
+            fa = this->wavefunction.at(this->nbox);
 
             if (fa * fx1 < 0.)
-                Emax = Emiddle;
+                e_max = energy_middle;
         }
     }
 
-    std::cerr<< "ERROR: Solution not found in bisec_Numer " << wavefunction[nbox] << " > " << err << std::endl;
+    std::cerr << "ERROR: Solution not found using the bisection method, " << this->wavefunction.at(this->nbox) << " > " << err << std::endl;
 //    exit(9);
-    return Emiddle;
+    return energy_middle;
+}
+
+
+
+void Numerov::printToFile() {
+  std::ofstream myfile ("wavefunction.dat");
+  if (myfile.is_open())
+  {
+
+    for(int i = 0; i < this->wavefunction.size(); i ++){
+        myfile << i <<" " << this->wavefunction.at(i)<< std::endl ;
+    }
+
+    myfile.close();
+  }
+}
+
+double Numerov::getSolutionEnergy() {
+    return this->solutionEnergy;
 }
+
+std::vector<double> Numerov::getWavefunction() {
+    return this->wavefunction;
+}
+
+std::vector<double> Numerov::getProbability() {
+    return this->probability;
+}

src/Solver/Numerov.h

@@ -14,13 +14,36 @@
 #include <iostream>
 #include <vector>
 #include <string>
+#include <fstream> 
 #include <Potential.h>
-#include <fstream>
-        void wavefuctionToFile(double *wavefunction);
 
-        double trap_array(int, int, double, double *);
-        void fsol_Numerov(double, int, Potential, double *);
-        double solve_Numerov(double, double, double, int, Potential, double *);
-        double bisec_Numer(double, double, int, Potential, double *);
+class Numerov {
+        public:
+                Potential *potential;
+                int nbox;
+
+                Numerov(Potential, int);
+                std::vector<double> getWavefunction();
+                std::vector<double> getProbability();
+                double getSolutionEnergy();
+                double solve(double, double, double);
+                void printToFile();
+
+                /*! Integrate with the trapezoidal rule method, from a to b position in a function array*/
+                static double trapezoidalRule(int a, int b, double stepx, std::vector<double> function) {
+                        double sum = 0.0;
+                        for (int j = a + 1; j < b; j++) 
+                                sum += function.at(j);
+                        sum += (function.at(a) + function.at(b)) / 2.0;
+                        sum = (sum) * stepx;
+                        return sum;
+                }
+        private:
+                double solutionEnergy;
+                std::vector<double> wavefunction;
+                std::vector<double> probability;
+                void functionSolve(double);
+                double bisection(double, double);
+};
 
 #endif