Solved

I am getting many errors in the C++ code for solving schrodinger wave equation.

Posted on 2016-10-08
9
41 Views
Last Modified: 2016-11-01
code is:
/* osci.cpp: Solution of the one-dimensional Schrodinger equation for
a particle in a harmonic potential, using the shooting method.
To compile and link with gnu compiler, type: g++ -o osci osci.cpp
To run the current C++ program, simply type: osci
Plot by gnuplot: /GNUPLOT> set terminal windows
/GNUPLOT> plot "psi-osc.dat" with lines */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
int main(int argc, char*argv[])
{// Runtime constants
const static double Epsilon = 1e-10; // Defines the precision of
//... energy calculations
//APPENDIX C. C++ CODE FOR SOLVING THE SCHRODINGER EQUATION
const static int N_of_Divisions = 1000;
const static int N_max = 5; //Number of calculated Eigenstates
FILE *Wavefunction_file, *Energy_file, *Potential_file;
Wavefunction_file = fopen("psi-osc.dat","w");
Energy_file = fopen("E_n_Oszillator.dat","w");
Potential_file = fopen("HarmonicPotentialNoDim.dat", "w");
if (!(Wavefunction_file && Energy_file && Potential_file))
{ printf("Problems to create files output.\n"); exit(2); }
/* Physical parameters using dimensionless quantities.
ATTENTION: We set initially: hbar = m = omega = a = 1, and
reintroduce physical values at the end. According to Eq.(4.117),
the ground state energy then is E_n = 0.5. Since the wave function
vanishes only at -infinity and +infinity, we have to cut off the
calculation somewhere, as given by ’xRange’. If xRange is chosen
too large, the open (positive) end of the wave function can
diverge numerically in this simple shooting approach. */
const static double xRange = 12; // xRange=11.834 corresponds to a
//... physical range of -20fm < x < +20fm, see after Eq.(4.199).
const static double h_0 = xRange / N_of_Divisions;
double* E_pot = new double[N_of_Divisions+1];
double dist;
for (int i = 0; i <= N_of_Divisions; ++i)
{ // Harmonic potential, as given in Eq. (4.115), but dimensionless
dist = i*h_0 - 0.5*xRange;
E_pot[i] = 0.5*dist*dist; // E_pot[i]=0;//E_pot=0:Infinite Well!
fprintf(Potential_file, "%16.12e \t\t %16.12e\n", dist, E_pot[i]);
}
fclose(Potential_file);
/* Since the Schrodinger equation is linear, the amplitude of the
wavefunction will be fixed by normalization.
At left we set it small but nonzero. */
const static double Psi_left = 1.0e-3; // left boundary condition
const static double Psi_right = 0.0; // right boundary condition
double *Psi, *EigenEnergies;// Arrays to hold the results
Psi = new double[N_of_Divisions+1]; //N_of_Points = N_of_Divisions+1
EigenEnergies = new double[N_max+1];
Psi[0] = Psi_left;
Psi[1] = Psi_left + 1.0e-3; // Add arbitrary small value
int N_quantum;//N_quantum is Energy Quantum Number
int Nodes_plus; // Number of nodes (+1) in wavefunction
double K_square;// Square of wave vector
// Initial Eigen-energy search limits
double E_lowerLimit = 0.0;// Eigen-energy must be positive
double E_upperLimit = 10.0;
int End_sign = -1;
char Limits_are_defined ="false";
double Normalization_coefficient;
double E_trial;
// MAIN LOOP begins:-----------------------------------
for(N_quantum=1; N_quantum <= N_max; ++N_quantum)
{
// Find the eigen-values for energy. See theorems (4.1) and (4.2).
Limits_are_defined ="false";
while (Limits_are_defined == false)
{ /* First, determine an upper limit for energy, so that the wavefunction
Psi[i] has one node more than physically needed. */
Nodes_plus = 0;
E_trial = E_upperLimit;
for (int i=2; i <= N_of_Divisions; ++i)
{ K_square = 2.0*(E_trial - E_pot[i]);
// Now use the NUMEROV-equation (4.197) to calculate wavefunction
Psi[i] = 2.0*Psi[i-1]*(1.0 - (5.0*h_0*h_0*K_square / 12.0))
/(1.0 + (h_0*h_0*K_square/12.0))-Psi[i-2];
if (Psi[i]*Psi[i-1] < 0) ++Nodes_plus;
}
/* If one runs into the following condition, the modification
of the upper limit was too aggressive. */
if (E_upperLimit < E_lowerLimit)
E_upperLimit = MAX(2*E_upperLimit, -2*E_upperLimit);
if (Nodes_plus > N_quantum) E_upperLimit *= 0.7;
else if (Nodes_plus < N_quantum) E_upperLimit *= 2.0;
else Limits_are_defined =" true"; // At least one node should appear.
} // End of the loop: while (Limits_are_defined == false)
// Refine the energy by satisfying the right boundary condition.
End_sign = -End_sign;
while ((E_upperLimit - E_lowerLimit) > Epsilon)
{ E_trial = (E_upperLimit + E_lowerLimit) / 2.0;
for (int i=2; i <= N_of_Divisions; ++i)
{ // Again eq.(4.197) of the Numerov-algorithm:
K_square = 2.0*(E_trial - E_pot[i]);
Psi[i] = 2.0*Psi[i-1] * (1.0 - (5.0*h_0*h_0*K_square / 12.0))
/(1.0 + (h_0*h_0*K_square/12.0))-Psi[i-2];
}
if (End_sign*Psi[N_of_Divisions] > Psi_right) E_lowerLimit = E_trial;
else E_upperLimit = E_trial;
} // End of loop: while ((E_upperLimit - E_lowerLimit) > Epsilon)
//APPENDIX C. C++ CODE FOR SOLVING THE SCHRÖDINGER EQUATION
// Initialization for the next iteration in main loop
E_trial = (E_upperLimit+E_lowerLimit)/2;
EigenEnergies[N_quantum] = E_trial;
E_upperLimit = E_trial;
E_lowerLimit = E_trial;
// Now find the normalization coefficient
double Integral = 0.0;
for (int i=1; i <= N_of_Divisions; ++i)
{ // Simple integration
Integral += 0.5*h_0*(Psi[i-1]*Psi[i-1]+Psi[i]*Psi[i]);
}
Normalization_coefficient = sqrt(1.0/Integral);
// Output of normalized dimensionless wave function
for (int i=0; i <=N_of_Divisions; ++i)
{ fprintf(Wavefunction_file, "%16.12e \t\t %16.12e\n",
i*h_0 - 0.5*xRange, Normalization_coefficient*Psi[i]);
}
fprintf(Wavefunction_file,"\n");
} // End of MAIN LOOP. --------------------------------
fclose(Wavefunction_file);
/*Finally convert dimensionless units in real units. Note that
energy does not depend explicitly on the particle’s mass anymore:
hbar = 1.05457e-34;// Planck constant/2pi
omega = 5.34e21; // Frequency in 1/s
MeV = 1.602176487e-13; // in J
The correct normalization would be hbar*omega/MeV = 3.5148461144,
but we use the approximation 3.5 for energy-scale as in chap. 4.9 */
const static double Energyscale = 3.5;// in MeV
// Output with rescaled dimensions; assign Energy_file
printf("Quantum Harmonic Oscillator, program osci.cpp\n");
printf("Energies in MeV:\n");
printf("n \t\t E_n\n");
for (N_quantum=1; N_quantum <= N_max; ++N_quantum)
{ fprintf(Energy_file,"%d \t\t %16.12e\n", N_quantum-1,
Energyscale*EigenEnergies[N_quantum]);
printf("%d \t\t %16.12e\n", N_quantum-1,
Energyscale*EigenEnergies[N_quantum]);
}
fprintf(Energy_file,"\n");
fclose(Energy_file);
printf("Wave-Functions in File: psi_osc.dat \n");
printf("\n");
return 0;
}

Open in new window

0
Comment
Question by:Pooja Soni
  • 3
  • 3
  • 2
  • +1
9 Comments
 
LVL 32

Assisted Solution

by:phoffric
phoffric earned 150 total points (awarded by participants)
ID: 41835521
You should indent your code.
You should copy the errors and post them.
0
 
LVL 30

Accepted Solution

by:
Zoppo earned 300 total points (awarded by participants)
ID: 41835565
I tried to compile the code and found only one compiler error:
char Limits_are_defined = "false";

Open in new window

This line tries to initialize a single char with a string what's not possible. You should change this line to use a bool and remove "-s around all occurances of true and false in the rest of the code:
bool Limits_are_defined = false;

Open in new window

After changing this I still got warnings/errors, but they only are about using function fopen which is considered to be unsafe, you should use fopen_s instead.

Hope that helps,

ZOPPO

PS: I agree with phoffric, you should at least post the error message with this kind of questions, and indenting the code (and even better using a CODE-blocks) makes it much easier to read/understand such a question.
0
 
LVL 32

Expert Comment

by:sarabande
ID: 41838185
No points please.

i used my IDE to indent your code and did what Andy has suggested:

#define _CRT_SECURE_NO_WARNINGS // only needed for VC compiler

/* osci.cpp: Solution of the one-dimensional Schrodinger equation for
a particle in a harmonic potential, using the shooting method.
To compile and link with gnu compiler, type: g++ -o osci osci.cpp
To run the current C++ program, simply type: osci
Plot by gnuplot: /GNUPLOT> set terminal windows
/GNUPLOT> plot "psi-osc.dat" with lines */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX(a, b) (((a) > (b)) ? (a) : (b))

int main(int argc, char*argv[])
{// Runtime constants
    const static double Epsilon = 1e-10; // Defines the precision of
    //... energy calculations
    //APPENDIX C. C++ CODE FOR SOLVING THE SCHRODINGER EQUATION
    const static int N_of_Divisions = 1000;
    const static int N_max = 5; //Number of calculated Eigenstates
    FILE *Wavefunction_file, *Energy_file, *Potential_file;
    Wavefunction_file = fopen("psi-osc.dat","w");
    Energy_file = fopen("E_n_Oszillator.dat","w");
    Potential_file = fopen("HarmonicPotentialNoDim.dat", "w");
    if (!(Wavefunction_file && Energy_file && Potential_file))
    { printf("Problems to create files output.\n"); exit(2); }
    /* Physical parameters using dimensionless quantities.
    ATTENTION: We set initially: hbar = m = omega = a = 1, and
    reintroduce physical values at the end. According to Eq.(4.117),
    the ground state energy then is E_n = 0.5. Since the wave function
    vanishes only at -infinity and +infinity, we have to cut off the
    calculation somewhere, as given by ’xRange’. If xRange is chosen
    too large, the open (positive) end of the wave function can
    diverge numerically in this simple shooting approach. */
    const static double xRange = 12; // xRange=11.834 corresponds to a
    //... physical range of -20fm < x < +20fm, see after Eq.(4.199).
    const static double h_0 = xRange / N_of_Divisions;
    double* E_pot = new double[N_of_Divisions+1];
    double dist;
    for (int i = 0; i <= N_of_Divisions; ++i)
    { // Harmonic potential, as given in Eq. (4.115), but dimensionless
        dist = i*h_0 - 0.5*xRange;
        E_pot[i] = 0.5*dist*dist; // E_pot[i]=0;//E_pot=0:Infinite Well!
        fprintf(Potential_file, "%16.12e \t\t %16.12e\n", dist, E_pot[i]);
    }
    fclose(Potential_file);
    /* Since the Schrodinger equation is linear, the amplitude of the
    wavefunction will be fixed by normalization.
    At left we set it small but nonzero. */
    const static double Psi_left = 1.0e-3; // left boundary condition
    const static double Psi_right = 0.0; // right boundary condition
    double *Psi, *EigenEnergies;// Arrays to hold the results
    Psi = new double[N_of_Divisions+1]; //N_of_Points = N_of_Divisions+1
    EigenEnergies = new double[N_max+1];
    Psi[0] = Psi_left;
    Psi[1] = Psi_left + 1.0e-3; // Add arbitrary small value
    int N_quantum;//N_quantum is Energy Quantum Number
    int Nodes_plus; // Number of nodes (+1) in wavefunction
    double K_square;// Square of wave vector
    // Initial Eigen-energy search limits
    double E_lowerLimit = 0.0;// Eigen-energy must be positive
    double E_upperLimit = 10.0;
    int End_sign = -1;
    bool Limits_are_defined = false;
    double Normalization_coefficient;
    double E_trial;
    // MAIN LOOP begins:-----------------------------------
    for(N_quantum=1; N_quantum <= N_max; ++N_quantum)
    {
        // Find the eigen-values for energy. See theorems (4.1) and (4.2).
        Limits_are_defined = false;
        while (Limits_are_defined == false)
        { /* First, determine an upper limit for energy, so that the wavefunction
          Psi[i] has one node more than physically needed. */
            Nodes_plus = 0;
            E_trial = E_upperLimit;
            for (int i=2; i <= N_of_Divisions; ++i)
            { K_square = 2.0*(E_trial - E_pot[i]);
            // Now use the NUMEROV-equation (4.197) to calculate wavefunction
            Psi[i] = 2.0*Psi[i-1]*(1.0 - (5.0*h_0*h_0*K_square / 12.0))
                /(1.0 + (h_0*h_0*K_square/12.0))-Psi[i-2];
            if (Psi[i]*Psi[i-1] < 0) ++Nodes_plus;
            }
            /* If one runs into the following condition, the modification
            of the upper limit was too aggressive. */
            if (E_upperLimit < E_lowerLimit)
                E_upperLimit = MAX(2*E_upperLimit, -2*E_upperLimit);
            if (Nodes_plus > N_quantum) E_upperLimit *= 0.7;
            else if (Nodes_plus < N_quantum) E_upperLimit *= 2.0;
            else Limits_are_defined = true; // At least one node should appear.
        } // End of the loop: while (Limits_are_defined == false)
        // Refine the energy by satisfying the right boundary condition.
        End_sign = -End_sign;
        while ((E_upperLimit - E_lowerLimit) > Epsilon)
        { E_trial = (E_upperLimit + E_lowerLimit) / 2.0;
        for (int i=2; i <= N_of_Divisions; ++i)
        { // Again eq.(4.197) of the Numerov-algorithm:
            K_square = 2.0*(E_trial - E_pot[i]);
            Psi[i] = 2.0*Psi[i-1] * (1.0 - (5.0*h_0*h_0*K_square / 12.0))
                /(1.0 + (h_0*h_0*K_square/12.0))-Psi[i-2];
        }
        if (End_sign*Psi[N_of_Divisions] > Psi_right) E_lowerLimit = E_trial;
        else E_upperLimit = E_trial;
        } // End of loop: while ((E_upperLimit - E_lowerLimit) > Epsilon)
        //APPENDIX C. C++ CODE FOR SOLVING THE SCHRÖDINGER EQUATION
        // Initialization for the next iteration in main loop
        E_trial = (E_upperLimit+E_lowerLimit)/2;
        EigenEnergies[N_quantum] = E_trial;
        E_upperLimit = E_trial;
        E_lowerLimit = E_trial;
        // Now find the normalization coefficient
        double Integral = 0.0;
        for (int i=1; i <= N_of_Divisions; ++i)
        { // Simple integration
            Integral += 0.5*h_0*(Psi[i-1]*Psi[i-1]+Psi[i]*Psi[i]);
        }
        Normalization_coefficient = sqrt(1.0/Integral);
        // Output of normalized dimensionless wave function
        for (int i=0; i <=N_of_Divisions; ++i)
        { fprintf(Wavefunction_file, "%16.12e \t\t %16.12e\n",
        i*h_0 - 0.5*xRange, Normalization_coefficient*Psi[i]);
        }
        fprintf(Wavefunction_file,"\n");
    } // End of MAIN LOOP. --------------------------------
    fclose(Wavefunction_file);
    /*Finally convert dimensionless units in real units. Note that
    energy does not depend explicitly on the particle’s mass anymore:
    hbar = 1.05457e-34;// Planck constant/2pi
    omega = 5.34e21; // Frequency in 1/s
    MeV = 1.602176487e-13; // in J
    The correct normalization would be hbar*omega/MeV = 3.5148461144,
    but we use the approximation 3.5 for energy-scale as in chap. 4.9 */
    const static double Energyscale = 3.5;// in MeV
    // Output with rescaled dimensions; assign Energy_file
    printf("Quantum Harmonic Oscillator, program osci.cpp\n");
    printf("Energies in MeV:\n");
    printf("n \t\t E_n\n");
    for (N_quantum=1; N_quantum <= N_max; ++N_quantum)
    { fprintf(Energy_file,"%d \t\t %16.12e\n", N_quantum-1,
    Energyscale*EigenEnergies[N_quantum]);
    printf("%d \t\t %16.12e\n", N_quantum-1,
        Energyscale*EigenEnergies[N_quantum]);
    }
    fprintf(Energy_file,"\n");
    fclose(Energy_file);
    printf("Wave-Functions in File: psi_osc.dat \n");
    printf("\n");
    return 0;
}

Open in new window


the program runs and the output is:

Quantum Harmonic Oscillator, program osci.cpp
Energies in MeV:
n                E_n
0                1.749999999795e+000
1                5.249999998112e+000
2                8.749999992829e+000
3                1.224999998232e+001
4                1.574999996759e+001
Wave-Functions in File: psi_osc.dat

Open in new window


Sara
0
 

Author Comment

by:Pooja Soni
ID: 41838579
thank you but I'm still getting the following error while compiling with cygwin

DELL@DELL2 ~
$ ./osci.cpp
./osci.cpp: line 1: /bin: Is a directory
./osci.cpp: line 2: a: command not found
./osci.cpp: line 3: To: command not found
./osci.cpp: line 4: To: command not found
./osci.cpp: line 5: Plot: command not found
./osci.cpp: line 6: /GNUPLOT: No such file or directory
./osci.cpp: line 11: $'\r': command not found
./osci.cpp: line 12: syntax error near unexpected token `('
'/osci.cpp: line 12: `int main(int argc, char*argv[])
0
Highfive + Dolby Voice = No More Audio Complaints!

Poor audio quality is one of the top reasons people don’t use video conferencing. Get the crispest, clearest audio powered by Dolby Voice in every meeting. Highfive and Dolby Voice deliver the best video conferencing and audio experience for every meeting and every room.

 
LVL 32

Assisted Solution

by:phoffric
phoffric earned 150 total points (awarded by participants)
ID: 41838686
Here is one way to compile and build an executable named a.exe in Cygwin:

$ g++ osci.cpp

To run the program:
$ ./a.exe
Quantum Harmonic Oscillator, program osci.cpp
Energies in MeV:
n                E_n
0                1.749999999795e+00
1                5.249999998112e+00
2                8.749999992829e+00
3                1.224999998232e+01
4                1.574999996759e+01
Wave-Functions in File: psi_osc.dat
0
 
LVL 32

Assisted Solution

by:sarabande
sarabande earned 50 total points (awarded by participants)
ID: 41840555
/* osci.cpp: Solution of the one-dimensional Schrodinger equation for
a particle in a harmonic potential, using the shooting method.
To compile and link with gnu compiler, type: g++ -o osci osci.cpp
To run the current C++ program, simply type: osci
Plot by gnuplot: /GNUPLOT> set terminal windows
/GNUPLOT> plot "psi-osc.dat" with lines */

the above lines are between /* ... */ what means that the lines are comments. so they should be ignored both by the precompiler and the compiler.

i don't know any c or c++ compiler which would try to compile these comment lines.

did you try the build command suggested by phoffric?

note, at Windows you may consider to using visual studio rather than cygwin. especially if the code is not from unix and/or doesn't use unix specials like fork.

you can download visual studio 2015 community for free. then create a new empty c++ win32 console project (with precompiled headers disabled). add the osci.cpp (i posted) to the project and you are done.

Sara
0
 

Author Comment

by:Pooja Soni
ID: 41840559
yes I have tried the solution suggested by phoffric....
but but it is giving errors as above..........
0
 
LVL 32

Assisted Solution

by:phoffric
phoffric earned 150 total points (awarded by participants)
ID: 41840655
$ ./osci.cpp

You got those errors because you are trying to run a file that is source code.
C++ is not an interpretive language.

You have to build it as I showed you. After building it using g++, take a look in the same directory and see whether you have the executable named a.exe and if you do, then you can run it by typing
./a.exe at the dollar sign prompt. Note that there is a dot at the beginning of the command followed by a forward slash followed by the name of the executable, in this case the default name, a.exe.
0
 
LVL 32

Expert Comment

by:sarabande
ID: 41868063
The code compiles without problems after applying Zoppo's correction.

The Questioner has problems with their developer environment which were addressed by Phoffric and answered in further questions the Author has asked for.

Sara
0

Featured Post

Highfive + Dolby Voice = No More Audio Complaints!

Poor audio quality is one of the top reasons people don’t use video conferencing. Get the crispest, clearest audio powered by Dolby Voice in every meeting. Highfive and Dolby Voice deliver the best video conferencing and audio experience for every meeting and every room.

Join & Write a Comment

When writing generic code, using template meta-programming techniques, it is sometimes useful to know if a type is convertible to another type. A good example of when this might be is if you are writing diagnostic instrumentation for code to generat…
Introduction This article is the first in a series of articles about the C/C++ Visual Studio Express debugger.  It provides a quick start guide in using the debugger. Part 2 focuses on additional topics in breakpoints.  Lastly, Part 3 focuses on th…
The goal of the video will be to teach the user the difference and consequence of passing data by value vs passing data by reference in C++. An example of passing data by value as well as an example of passing data by reference will be be given. Bot…
The viewer will learn how to user default arguments when defining functions. This method of defining functions will be contrasted with the non-default-argument of defining functions.

747 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question

Need Help in Real-Time?

Connect with top rated Experts

11 Experts available now in Live!

Get 1:1 Help Now