Problem 3

Ideas

To solve this problem, I first construct the Halmiltonian Matrix to the problem.

The matrix elements [𝐀]nn′=  𝐀|n′>  

In this problem, 𝐀= 𝐇𝐇𝐚𝐫𝐦𝐨𝐧𝐢𝐜+ 𝐕𝐦𝐨𝐫𝐬𝐞− 𝐕𝐇𝐚𝐫𝐦𝐨𝐧𝐢𝐜 

Thus, the element [𝐀]nn′= < n|(𝐇𝐇𝐚𝐫𝐦𝐨𝐧𝐢𝐜+ 𝐕𝐦𝐨𝐫𝐬𝐞− 𝐕𝐇𝐚𝐫𝐦𝐨𝐧𝐢𝐜)|n′>  = <n|𝐇𝐇𝐚𝐫𝐦𝐨𝐧𝐢𝐜|n′> + <n|(𝐕𝐦𝐨𝐫𝐬𝐞− 𝐕𝐇𝐚𝐫𝐦𝐨𝐧𝐢𝐜)|n′>  

Since all the eigenstates of Harmonic oscillator are orthonormal, 𝐇𝐚𝐫𝐦𝐨𝐧𝐢𝐜|n′> = 𝐄𝐧𝛅𝐧𝐧′ 

Yet, I still need to figure out the potential part. There are two ways to work this out. One approach is to use operators a and a+. 𝐱=   (h/4πμω)^1/2( 𝐚+ 𝐚+ )/2

So, e^(−β𝐱)= e^(−β (h/4πμω)^1/2( 𝐚+ 𝐚+ )/2)   = e^(−β (h/4πμω)^1/2𝐚/2) ×e^(-β (h/4πμω)^1/2𝐚+/2)   = 1/i!(−β (h/4πμω)^1/2𝐚/2)^i *1/j!(-β (h/4πμω)^1/2𝐚+/2) j 

Then, a is a lower operator and a+ is a rising operator, if we want   𝐕|n′>  to be non-zero, the difference between the number of a and the number of a+ should be equal to  n−n′  (for the case n is less than n′, vice versa). This is the basic idea for the code.

The other method is to calculate the integral numerically. The classical approximation is to divide the whole integral space into small bins. For each bin, the exact result can be replaced by the area of trapezia. To sum over all the cells, the integral can be achieved.

I used the second one to answer this question. For all the coefficients, I used the Problem Set4 as reference. De=7.31×10^-19 J∙molecule^-1 β=1.82×10¹⁰m^-1

After I have build the Halmiltonian Matrix, the subsequent issue is to diagonalize this matrix. The eigenvalues are the trace of the diagonalized matrix. In order to transform the matrix to the desired form, I apply the classical elimination tool. The algorithm is as follows:

1) Find the first column that contains non-zero coefficient,

0  a

0    c

0    b

e.g. c, then swap the row including c with the first row and cancel all other elements (a, b) in this column by multiple some value to c.

2) Repeat 1) until it reaches the last column, step 1) once more;

Or, it reaches the last row.

3) Since the Matrix has been converted to a structure like,

a    d

0    b

0    0    c

in which all the elements at left bottom is zero and the top right contains non-zero members. We can use the right-down corner one to clear all the above in this column.

4) Repeat 3) until the we meet the first row.

Now, trace of the matrix shows the eigenvalues for this Halmitonian.

 

Codes

The platform for the execution is Microsoft Visual C++ 2008.

1) Using operators a and a+.

#include "stdafx.h"

#include "stdlib.h"

#include "iostream"

#include "math.h"

 

using namespace std;

 

const float beta = 1.0; //β. This is only for editing the code, not the real value in this problem.

 

class Integral{

 float *coefficient;

public:

Integral();

 void setCoefficient(void);

 int f(int i);

 float computeCoefficient(int longth);

 void getCoefficient();

~Integral();

};

 

Integral::Integral(){

setCoefficient();

}

 

void Integral::setCoefficient(void){

coefficient = new float [10];

 for (int i = 0; i < 10; i++){

coefficient[i] = pow(beta, i) / f(i);

}

} // Set the coefficients for the polynomial of Taylor expansion.

 

int Integral::f(int i){

 int result = 1;

 for(int j = 0; j < i; j++){

result = result * (j+1);

}

 return result;

}

 

float Integral::computeCoefficient(int longth){

 float result = 0;

 int temp = 0;

 for (int i = 0; i < (10-longth); i++){

result = result + coefficient[i] * coefficient[i+longth];

temp++;

}

 return result;

}

// Compute the result of integral. When the operators function on the states n and n′, it would

// change them to the same state. Since  n|n  is 1, we only need to consider the coefficients of

// two states n and n′.

 

void Integral::getCoefficient(void){

 for (int i = 0; i < 10; i++)

cout << coefficient[i] << " ";

cout << endl;

}

 

Integral::~Integral(){

 delete [] coefficient;

}

 

2) The numerical integral.

#include "stdafx.h"

#include "stdlib.h"

#include "iostream"

#include "math.h"

 

using namespace std;

 

/*

H_0(x) = 1

H_1(x) = 2x

H_2(x) = 4x^2-2

H_3(x) = 8x^3-12x

H_4(x) = 16x^4-48x^2+12

H_5(x) = 32x^5-160x^3+120x

H_6(x) = 64x^6-480x^4+720x^2-120

H_7(x) = 128x^7-1344x^5+3360x^3-1680x

H_8(x) = 256x^8-3584x^6+13440x^4-13440x^2+1680

H_9(x) = 512x^9-9216x^7+48384x^5-80640x^3+30240x

These are Hermite Polynomials.

*/

const double De = 0.731e-18;

const double E = 2.71; // E is the natural e.

const double Alpha = 1.644e20;

const double Beta = 1.82e10;

const double Pi = 3.14159265;

 

class NIntegral{

 int original[10];

 int use[10];

public:

NIntegral();

 void derivative(int n);

 void getPolynomial(void);

 double ex(double);

 double potentialMorse(double);

 double potentialHarmonic(double);

 double normal(int);

 double compute(int i, int j, double x);

};

 

NIntegral::NIntegral(){

original[0] = original[2] = original[4] = original[6] = original[8] = 0;

original[1] = 30240; original[3] = -80640; original[5] = 48384;

original[7] = -9216; original[9] = 512;

 for (int i = 0; i < 10; i++)

use[i] = 0;

}

// Because the integral can always have a constant, I use the derivative instead.

// By hiring the recursion relation between Hermite Polynomials, dHn/dx =2nHn−1 

// we can figure out all the lower order polynomials from H9 by taking derivative.

 

void NIntegral::derivative(int n){

 for (int i = 0; i < 10; i++)

use[i] = original[i];

 int temp[10];

 for (int i = 0; i < n; i++){

 for (int j = 1; j < 10; j++)

temp[j-1] = use[j] * j / (9-i) / 2;

temp[9] = 0;

 for (int j = 0; j < 10; j++)

use[j] = temp[j];

}

}

 

void NIntegral::getPolynomial(void){

 for (int i = 0; i < 10; i++){

derivative(i);

cout << "this is " << i << "th derivative" << endl;

 for (int j = 0; j < 10; j++){

 if (use[j]){

cout << use[j] << " x^" << j << "+";

}

}

cout << endl;

}

}

 

double NIntegral::ex(double x){

 double o = 0 - x * x * Alpha / 2;

 return pow(E, o);

} // e^−αx²

 

double NIntegral::potentialMorse(double x){

 double result;

result = 0- Beta * x;

result = pow(E, result);

result = 1 - result;

result = result * result * De;

 return result;

} // De (1−e^(−βx))²

 

double NIntegral::potentialHarmonic(double x){

 return De * Beta * Beta * x * x;

}// 1/2kx²

 

double NIntegral::normal(int n){

 double k = Alpha/Pi;

 double a = pow( double (k) , double (0.25));

 double b = 1.0;

 for (int i = 0; i < n; i++){

b = b * 2 * (i+1);

}

 double nor = 1 / pow(double(b), double (0.5)) * a;

 return nor;

} // 1/(2ⁿn!)^1/2 

 

double NIntegral::compute(int i, int j, double x){

 double h1, h2, v;

 double a = pow( double (Alpha), double (0.5));

h1 = ex(x) * normal(i);

h2 = ex(x) * normal(j);

derivative(9-i);

 double temp = 0;

 for (int t = 0; t < 10; t++){

 if(use[t]){

temp = temp + use[t] * pow(a*x , t);

}

}

h1 = temp * h1;

derivative(9-j);

temp = 0;

 for (int t = 0; t < 10; t++){

 if(use[t]){

temp = temp + use[t] * pow(a*x , t);

}

}

h2 = temp * h2;

 double result = 1.0;

result = result * h1 * h2 * (potentialMorse(x) - potentialHarmonic(x));

 return result;

} // Compute the one point result with special states i, j and sepcial coordinate x.

 

3) Construct the Halmitonian Matrix

#include "stdafx.h"

#include "stdlib.h"

#include "iostream"

#include "math.h"

 

using namespace std;

 

#include "integral.cpp"

 

const double Max = 1000; // Integral up limit.

const double Min = 0; // Integral low limit.

const int N = 1000; // The number of bins.

 

class Matrix{

 double matrice[10][10]; // Store the matrix.

public:

Matrix();

 void construct(void);

 void getMatrix(void);

};

 

Matrix::Matrix(){

 for (int i = 0; i < 10; i++){

 for (int j = 0; j < 10; j++){

matrice[i][j] = 0;

}

}

}

 

void Matrix::construct(void){

NIntegral integral;

 double top, bottom, length;

length = (Max - Min) / N / 2;

 for (int i = 0; i < 10; i++){

 for (int j = 0; j < 10; j++){

 for (int t = 0; t < 2 * N; t++){

top = Min + length * t;

bottom = top + length;

matrice[i][j] = matrice[i][j] + (integral.compute(i, j, top) + integral.compute(i, j, bottom)) * length / 2;

// Compute one point by calculating the area of trapezia.

}

 if (j == i){

matrice[i][j] = matrice[i][j] + 0.3997 * ( i + 0.5 );

// For the diagonal elements, the Harmonic part have been added.

// h/2π (k/m)^1/2  = 0.3997e-18

}

}

}

}

 

void Matrix::getMatrix(void){

 for (int i = 0; i < 10; i++){

 for (int j = 0; j < 10; j++){

cout << matrice[i][j] << "\t";

}

cout << endl;

}

}

 

int _tmain(int argc, _TCHAR* argv[])

{

Matrix mat;

mat.construct();

mat.getMatrix();

 

NIntegral integ;

cout << integ.compute(9, 9, 0.00001) << endl;

getchar();

 return 0;

}

 

4) Diagonize the Matrix.

#include "stdafx.h"

#include "stdlib.h"

#include "iostream"

#include "math.h"

 

using namespace std;

#define ROW 10 // the number of equation

#define UNKNOWN 11 // UNKNOWN = unknown + result(1)

// Ideas are as follows:

// every time, I will find the first non-zero element for each equations.

// Compare their positions and find the lowest position.

// Swap the lowest to the first index.

// Reverse the substraction process.

 

class LinearEquationSolver{

private:

 float **eq; // store equations

 float** set; // store solutions

 int index;// record the pointer lively.

 int min;// record the start of non-zero term

public:

LinearEquationSolver();

 void setCoefficient(float* e1, float* e2, float* e3, float* e4, float* e5, float* e6, float* e7, float* e8, float* e9, float* e10);

 void getEquation(void);

 void getIndividual(int);

 int firstNonZero(float*);

 void solve(void);

 void minNonZero(void);

 void substract(void);

 void ultimate(void);

 int reverse(void);

 void solution(int);

~LinearEquationSolver();

};

 

LinearEquationSolver::LinearEquationSolver(){

eq = new float* [ROW];

 for (int j = 0; j < ROW; j++){

eq[j] = new float [UNKNOWN];

}

index = 0;// Search from the top.

min = 0;// The default first non-zero is the first value in the equation

}

// Initialize all equations.

// Initially, I set all other elements of all equations to be 0.

 

void LinearEquationSolver::setCoefficient(float* e1, float* e2, float* e3, float* e4, float* e5, float* e6, float* e7, float* e8, float* e9, float* e10){

 for (int i = 0; i < UNKNOWN; i++){

eq[0][i] = e1[i];

eq[1][i] = e2[i];

eq[2][i] = e3[i];

eq[3][i] = e4[i];

eq[4][i] = e5[i];

eq[5][i] = e6[i];

eq[6][i] = e7[i];

eq[7][i] = e8[i];

eq[8][i] = e9[i];

eq[9][i] = e10[i];

}

}

// Set all equations to be the desired ones.

 

LinearEquationSolver::~LinearEquationSolver(){

 delete [] eq;

}

 

void LinearEquationSolver::getEquation(void){

 for (int i = 0; i < ROW; i++)

getIndividual(i);

}

// Print out all equations

void LinearEquationSolver::getIndividual(int i){

 for (int j = 0; j < UNKNOWN-2; j++)

cout << eq[i][j] << "x + ";

cout << eq[i][UNKNOWN-2] << "x = " << eq[i][UNKNOWN-1] << endl;

}

 

int LinearEquationSolver::firstNonZero(float* tempeq){

 for (int i = 0; i < UNKNOWN-1; i++){

 if(abs(tempeq[i]))

 return i;

}

 return UNKNOWN-1;

}

// Find the first non-zero term in one equation.

// UNKNOWN-1 is the signal that all the coefficients are 0.

 

void LinearEquationSolver::solve(void){

 while ( min < (UNKNOWN-2) && (index < ROW-1)){

minNonZero();

substract();

index++;

}

 // The procedure is as follows:

 // 1) find the first non-zero terms in equations, swapping this equation to the first place (index).

 // 2) substract all below equations by the first equations to kill all the first non-zero terms.

 // 3) move the index to the next equation, continue the loops.

 // 4) until the code find the first non-zero term is the result,

 // 5) or the index is the last equation.

 

ultimate();

}

 

void LinearEquationSolver::minNonZero(void){

min = firstNonZero(eq[index]);

 int star = index;// to record the minimum row for swapping

 for(int i = index+1; i < ROW; i++){

 if( min > firstNonZero(eq[i])){

min = firstNonZero(eq[i]);

star = i;

}

}

 if (star == index)

 return;

 else{

 float temp;

 for (int i = 0; i < UNKNOWN; i++){

temp = eq[index][i];

eq[index][i] = eq[star][i];

eq[star][i] = temp;

}

 return;

}

}

 

void LinearEquationSolver::substract(void){

 float efficient, temp;

 int start = firstNonZero(eq[index]);// the first non-zero term in the row being operated.

 if (start == UNKNOWN-1)

 return;

 // If the first non-zero term is the result, then I don't want to continue.

 // The result may be 0 because I don't check that.

 

 for (int i = index + 1; i < ROW; i++){

efficient = eq[i][start] / eq[index][start];

 if (abs(efficient)){

 for (int j = start; j < UNKNOWN; j++){

temp = eq[index][j] * efficient;

eq[i][j] = eq[i][j] - temp;

 if (abs(eq[i][j]) < 0.000001)

eq[i][j] = 0;

}

}

}

 // I have substracted the first non-zero term in all equations.

}

 

void LinearEquationSolver::ultimate(void){

 float efficient;

 int subindex = reverse();

 if ((subindex + 1) == (UNKNOWN - 1)){

cout << " This equation is full order. There is only one solution." << endl;

}else{

cout << " The order of the equations: " << ( subindex + 1 ) << "." << endl;

}

 int marker = subindex;// the marker of how many dependences

 int j;

 float temp;

 for(; subindex > -1; subindex--){

j = firstNonZero(eq[subindex]);

 for (int i = subindex-1; i > -1; i--){

efficient = eq[i][j]/eq[subindex][j];

 if (efficient){

 for (int t = j; t < UNKNOWN; t++){

temp = eq[subindex][t] * efficient;

eq[i][t] = eq[i][t] - temp;

}

}

}

}

 // transform the equations to the standard form

 

solution(marker);

}

 

int LinearEquationSolver::reverse(void){

 for (int i = 0; i < ROW; i++){

 if (firstNonZero(eq[i]) == (UNKNOWN-1)){

 return (i-1);

}

 if (firstNonZero(eq[i]) == (UNKNOWN-2)){

 return i;

}

}

 return ROW-1;

 // There are three cases:

 // 1) 0x + 0x + 0x + 4x = 2, this the standard type.

 // 2) 0X + 0x + 0x + 0x = 0, this may appear when there are more equations than unknowns.

 // 3) 0x + 2x + 3x + 5x = 3, this may appear when there are less equations than unknowns.

}

 

void LinearEquationSolver::solution(int marker){

 int depend = (UNKNOWN - 1) - (1 + marker);

set = new float* [depend + 1];

 // there are "depend" basic solutions and one special one.

 

 for (int t = 0; t < depend + 1; t++){

set[t] = new float [UNKNOWN - 1];

}

 

 int* dependence;// store the position of dependence.

dependence = new int [depend];

 int d;

 for (int t = 0; t < (depend + 1); t++)

 for (d = 0; d < UNKNOWN; d++)

set[t][d] = 0;

 // Initialize the solution set.

d = 0;

 

 bool signal = false;

 for (int t = 0; t < (UNKNOWN - 1); t++){

 for (int k = 0; k < marker + 1; k++)

 if ( t == firstNonZero(eq[k])){

signal = true;

 break;

}

 if (!signal){

dependence[d] = t;

d++;

}

signal = false;

}

 // set the position of determinant unknown to the array dependence.

 

 int element;

 for (int i = 0; i < depend; i++){

set[i][dependence[i]] = 1;

 for (int j = 0; j < marker + 1; j++){

element = firstNonZero(eq[j]);

set[i][element] = ( 0 - eq[j][dependence[i]])/eq[j][element];

}

}

 // get the basic solutions

 

 for (int j = 0; j < marker + 1; j++){

element = firstNonZero(eq[j]);

set[depend][element] = eq[j][UNKNOWN - 1]/eq[j][element];

}

 // get the special solution

 

cout << " The solution are as follows:" << endl;

cout << " The basic solution are:" << endl;

 for (int i = 0; i < depend; i++){

 for (int j = 0; j < (UNKNOWN - 1); j++){

cout << set[i][j] << ", ";

}

cout << " the determinant is the " << dependence[i]+1 << "th value.";

cout << endl;

}

cout << " The special solution is:" << endl;

 for (int j = 0; j < (UNKNOWN - 1); j++){

cout << set[depend][j] << ", ";

}

cout << endl;

 

 delete [] set;

 delete [] dependence;

}

 

int _tmain(int argc, _TCHAR* argv[])

{

LinearEquationSolver linear;

 float e1[11] = {0.1999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

 float e2[11] = {0, 0.5996, 0, 0, 0, 0, 0, 0, 0, 0, 0};

 float e3[11] = {0, 0, 0.9995, 0, 0, 0, 0, 0, 0, 0, 0};

 float e4[11] = {0, 0, 0, 0.1.3993, 0, 0, 0, 0, 0, 0, 0};

 float e5[11] = {0, 0, 0, 0, 1.7991, 0, 0, 0, 0, 0, 0};

 float e6[11] = {0, 0, 0, 0, 0, 2.1989, 0, 0, 0, 0, 0};

 float e7[11] = {0, 0, 0, 0, 0, 0, 2.5987, 0, 0, 0, 0};

 float e8[11] = {0, 0, 0, 0, 0, 0, 0, 2.9985, 0, 0, 0};

 float e9[11] = {0, 0, 0, 0, 0, 0, 0, 0, 3.3983, 0, 0};

 float e10[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 3.7981, 0};

 

linear.setCoefficient(e1, e2, e3, e4, e5, e6, e7, e8, e9, e10);

linear.getEquation();

linear.solve();

linear.getEquation();

 return 0;

}

 

Result

0.1999 0 0 0 0 0 0 0 0 0

0 0.5996 0 0 0 0 0 0 0 0

0 0 0.9995 0 0 0 0 0 0 0

0 0 0 1.3993 0 0 0 0 0 0

0 0 0 0 1.7991 0 0 0 0 0

0 0 0 0 0 2.1989 0 0 0 0

0 0 0 0 0 0 2.5987 0 0 0

0 0 0 0 0 0 0 2.9985 0 0

0 0 0 0 0 0 0 0 3.3983 0

0 0 0 0 0 0 0 0 0 3.7981

This is the matrix after diagonalization. Because the halmitonian matrix is already of the type of diagonalized, they look very similar. Then, we can know the eigenvalues are: 0.1999e-18, 0.5996e-18, 0.9995e-18, 1.3993e-18, 1.7991e-18, 2.1989e-18, 2.5987e-18, 2.9985e-18, 3.3983e-18, 3.7981e-18.

The exact eigenvalues can be earned from the below formula (from Wikipedia),

En=hv0 (n+1/2) −  (hv0( n+1/2)²/4De  

v0= β/2π ( 2De/u)^1/2

The exact values are:

0.7853e-18, 0.47e-18, 0.657e-18, 0.7288e-18, 0.6910e-18, 0.54407e-18, 0.2878e-18, -0.0778e-18, -0.5527e-18, -1.1369e-18.

From these data, we can see: for the lower quantum numbers, the eigenvalues earned from harmonic approximation are close. But, the accurate ones are very different from the harmonic ones in the higher order region.

Discussion

1) From the data, we can see the potential difference between harmonic and morse are very small compared with the harmonic eigenenergy. I am not sure if this conclusion is right. If it is not correct, the reason may be that I don’t know if the case from Problem Set4 was suitable for this calculation. It is possible that some parameters are misused. However, I have test the diagonalization code and it works well. Thus, the main concern is about the choice of parameters.

2) Since harmonic oscillator is an approximation of morse oscillator at the bottom curve, they should be close around that part. The result accords with theoretical prediction. For the higher energy, harmonic oscillator diverges from morse oscillator.