Thread: Debugging Physics code is by far the worst code to debug!

  1. #1
    Registered User Swarvy's Avatar
    Join Date
    Apr 2008
    Location
    United Kingdom
    Posts
    195

    Question Debugging Physics code is by far the worst code to debug!

    Hello all, I've been writing a little physics code recently to solve the radial TISE for Hydrogen-like potentials using the Numerov Algorithm, basing my code on source code found on the internet. The program seems to work almost exactly as it should, however, if the radius is increased or the eigenvalue number is large then the algorithm seems to break down. I've included some of the key functions used in the technique.

    The 'Solve' Method - For a given energy finds the nearest eigenvalue - Which in the context of Schrodinger's Equation is the energy
    Code:
    double Numerov::Solve(bool dir)
    {
    	double dE = 1.0e-02;
    	double fE;
    	double fE1;
    	double E1, E2;
    	E1 = E + dE;
    	int it = 0;
    
    	// Secant method to improve energy eigenvalue
    	while (fabs(dE) > epsilon && it < 1000)
    	{
    		u1  = u0*exp(sqrt(-2.0*m_ec2*E)/h_bar_c*h);
    		iterate(u0, u1, E, dir);
    		fE = u[0];
    		u1  = u0*exp(sqrt(-2.0*m_ec2*E1)/h_bar_c*h);
    		iterate(u0, u1, E1, dir);
    		fE1 = u[0];
    		E2  = E1 - fE1*(E1 - E)/(fE1 - fE);
    		E   = E1;
    		E1  = E2;
    		dE  = E1 - E;
    		it++;
    	}
    	// calculate u with best energy value
    	iterate(u0, u1, E, dir);
    	Norm();
    
    	cout << endl;
    	cout << "==> Energy difference: Rydberg - E = " << Ry - E << "  after step: " << it << endl;
    	return E;
    }
    This is the 'Iterate' Method - For a given energy it calculates the wavefunction using Numerov's Method. The function 'Q' is the effective potential function. I know there are no problems with the 'Q' method.
    Code:
    long int Numerov::iterate(double y0, double y1, double E, bool direction)
    {
    	double hh12 = h*h/12.0;
    	double x_i, x_im, x_ip;
    
    	// Prograde recursion: This method is not working !
    	if (direction == true)
    	{
    		u[0] = 0.0;
    		u[1] = y0;
    		u[2] = y1;
    
    		for (long int i=2; i < NGrid-2; i++)
    		{
    			x_i  = a + i*h;
    			x_im = x_i - h;
    			x_ip = x_i + h;
    			u[i+1] = ( 2.0*(1.0 - 5.0*hh12*Q(x_i, E))*u[i] -
    			           (1.0 + hh12*Q(x_im, E))*u[i-1] ) / (1.0 + hh12*Q(x_ip, E));
    			// Test condition: increasing wave function after decreasing
    			if ((fabs(u[i-1]) > fabs(u[i])) && (fabs(u[i+1]) > fabs(u[i])))	return i;
    		}
    	}
    	// Retrograde recursion
    	else if (direction == false)
    	{
    		u[NGrid-1] = y0;
    		u[NGrid-2] = y1;
    
    		for (long int i=NGrid-2; i>=1; i--)
    		{
    			x_i  = a + i*h;
    			x_im = x_i - h;
    			x_ip = x_i + h;
    			u[i-1] = ( 2*(1.0 - 5.0*hh12*Q(x_i, E))*u[i] -
    			           (1.0 + hh12*Q(x_ip, E))*u[i+1] ) / (1.0 + hh12*Q(x_im, E));
    		}
    		return 0;
    	}
    }
    The array 'u' and 'v' represent the wavefunction and probability function respectively. 'h' is the step size, 'a' and 'b' are the starting and ending points considered (a is generally considered to be zero and so by 'large radii', I am referring to large 'b'). The data which it generates which is clearly wrong is given as output: '-1.#IND00000000'.

    Thanks in advance.
    Last edited by Swarvy; 03-10-2010 at 10:00 AM.

  2. #2
    Officially An Architect brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,396
    At higher eigenvalues the energy is greater and the wavelength is shorter. Have you tried decreasing the step size? What are the stability criteria for the Numerov algorithm?
    Code:
    //try
    //{
    	if (a) do { f( b); } while(1);
    	else   do { f(!b); } while(1);
    //}

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Replies: 0
    Last Post: 02-21-2002, 06:05 PM
  2. Debugging leads to buggy code and longer hours?
    By no-one in forum A Brief History of Cprogramming.com
    Replies: 6
    Last Post: 01-28-2002, 11:14 AM
  3. anyone bored enough to help me debug some code?
    By *ClownPimp* in forum C++ Programming
    Replies: 7
    Last Post: 01-18-2002, 07:55 PM
  4. Replies: 4
    Last Post: 01-16-2002, 12:04 AM
  5. debugging code
    By DeViLs_SouL in forum C++ Programming
    Replies: 3
    Last Post: 12-05-2001, 01:34 PM