Oh Well, Quantum Wells

083105109117108097116105110103032073110102105110105116121046046046

Our imagination is stretched to the utmost, not, as in fiction, to imagine things which are not really there, but just to comprehend those things which are there.

— Richard Feynman

It’s been a long hiatus since the last post but I hope to make it up for that here (with a source code instead of pseudo-code). So far, I have written about general algorithms. This time, I hope to explain the procedure involved in programming with an example where we will simulate quantum wells. I will assume that you have some basic knowledge of quantum mechanics and some basic programming experience with arrays and matrices. However, feel free to ask any question (no matter what) and I will do my best to clarify it.

We start by defining the problem. We want to study an electron trapped inside a one-dimensional quantum well. Within the quantum well, however, there are some potential barriers (think of other heavier ions) which influence the electron inside the well. We want quantum effects and so we use the Schrodinger equation:

i\hbar\frac{\partial\Psi}{\partial t} = -\frac{\hbar^2}{2m}\Delta\Psi + q\phi\Psi = \hat{H}\Psi

Here \Delta   is the second derivative with respect to space \left(\frac{\partial^2}{\partial x^2}\right)   and \phi   is the potential inside the quantum well. We are interested in calculating the energy of the electron inside the well at ground states and excited states along with the probability distribution (equivalent to charge distribution) at each energy level.

For the first step, we want that the electron remains trapped INSIDE the quantum well. That can be done if we assume an infinite potential barrier outside the well. In other words, \Psi   is zero outside the boundaries of the quantum well. This is good news because our “universe” is now confined to the quantum well.

We then look at what is required to be calculated — energy and probability distribution at each energy. The energy values are the “eigenvalues” and the probability distribution is given by the square of the wavefunction, the “eigenvectors”. There is one eigenvector for EACH eigenvalue. For more information, consult Eigenvalues and Eigenvectors. We now try to construct a matrix for the operator \hat{H}   (Hamiltonian matrix) and find its eigenvalues and eigenvectors. Each eigenvector would be a column vector. We also divide space (one-dimensional) into grid points. The number of matrix elements in the Hamiltonian matrix would be the square of the grid points. Before we construct our matrix, we will first simplify the double differentiation as

\frac{\partial^2\Psi_i}{\partial x^2} = \frac{\partial}{\partial x} \left( \frac{\Psi_{i+1}-\Psi_{i-1}}{2\delta} \right) = \frac{\frac{\Psi_{i+1}-\Psi_i}{\delta} - \frac{\Psi_i-\Psi_{i-1}}{\delta} }{2\delta} = \frac{\Psi_{i+1}-2\Psi_i+\Psi_{i-1}}{2\delta^2}

Here \delta   is the distance between two grid points. It defines the “fineness” of the mesh. Looking at the dependence on \Psi   , we see that the Hamiltonian will be zero everywhere except along the diagonal and two other “diagonals” forming a “tri-diagonal” matrix. Also, for the first row, we have \Psi_{i-1} = 0   and for the last row, we have \Psi_{i+1} = 0   because we have confined our “universe” to the boundaries of the quantum well. This is how our Hamiltonian matrix looks like:

\left(\begin{array}{cccccc}  \frac{\hbar^2}{2m\delta^2}+q\phi & \frac{1}{2\delta^2} & . & . & . & 0 \\  \frac{1}{2\delta^2}& \frac{\hbar^2}{2m\delta^2}+q\phi & \frac{1}{2\delta^2} & . & . & 0 \\  . & . & & & & .\\  . & & . & & & .\\  . & & & . & & .\\  0 & . & . & \frac{1}{2\delta^2}& \frac{\hbar^2}{2m\delta^2}+q\phi & \frac{1}{2\delta^2} \\  0 & . & . & . & \frac{1}{2\delta^2}& \frac{\hbar^2}{2m\delta^2}+q\phi  \end{array}\right)

The size of the matrix is N\times N   where N   is the number of grid points. Once this is done, you can use Lapack routines to diagonalize the matrix and provide you with the eigenvalues and eigenvectors.

I have written a quick program in Fortran for a particular example of a quantum well of about 500 nanometers and a potential barrier of -0.1 milli-volts between 175 and 325nm from the origin. The source code is provided here: qu_well

To compile, you will need Lapack routines. Download the required packages and you can compile this source code. If you use Linux with gfortran, you can compile the code and run the program using:


gfortran -o qu_well qu_well.f90 -llapack -lblas
./qu_well

You can change the value of the potential, add more barriers or change the sign of the barrier. Different simple systems can be simulated with this code and you don’t need a specialized laboratory with very expensive lasers to understand how an electron behaves. In case the program takes too long to run, you can reduce the number of grid points. At the end of the program, there are four files written out which contain the eigenvalues, the eigenvectors (for each eigenvalue), the probability distribution (also for each eigenvalue) and the potential barrier inside the quantum well. By careful plotting of the probability distribution near the potential barrier, you will also see that the electron with energy (eigenvalue) LOWER than the potential barrier still has a non-zero probability of being there. This is called Quantum Tunnelling and is a well-known and experimentally proven concept in scientific literature.

The aim of this post was to show how computer programming can be used to solve many real-world problems. That’s why I didn’t sit and analyze the results. The focus in this post may not have seemed to be computer programming but that’s because a process of tackling any problem is best understood with an example. The requirement to change the partial derivatives as a finite difference and the method of creating the matrix is so that the code can easily written. Computers work in bits and so it becomes necessary to modify “continuous” data into “discrete” data like we did with the derivatives and the grid. The Hamiltonian matrix is a “continuous” matrix (because of the derivatives) which makes it an infinite-sized 2-dimensional matrix. To compute its eigenvalues and eigenvectors, we create a mesh and only take the Hamiltonian at those points. The derivatives then become differences of the value between those points.

Quantum Mechanics is not the only field where one needs to compute the eigenvalues and eigenvectors. These are required in many other fields like Image Processing, Geology, Statistics and Rotation Studies. But the process of tackling a problem using computers remains similar. We first have to break our problem into something a computer can handle. I hope that the source code provided helps you to connect the dots with how a basic one-dimensional quantum well problem can be tackled with. Trying to replicate the code and then rewriting (at least modifying) similar code would enable one to understand the process on a more “hands-on” level.

068111110101046046046032087104101114101032097109032073063

Advertisements

Sort of… Sorts – I

083111114116105110103032100097116097046046046

Sorts and Big O

If you use a relatively modern computer, you are generally provided with a file browser. In case of Microsoft Windows, you would have an application called Explorer that allows you to sort your files and folders. In case of Ubuntu, you have Nautilus. Mac users generally use Finder.

You probably have already noticed that when you right-click on an empty section of your file browser, you are given an option to ‘Arrange Items’ or ‘Sort Items’ or something similar, which allows you to view the files and folders arranged by name, size, type or modification date.

The algorithms which help perform these tasks are called sorting algorithms. They sort items based on name, type, size and other criteria. A lot of studies have been done on sorting algorithms, which might seem surprising at first glance. However, if you imagine that you have a million items to sort, the time taken by non-recursive algorithms (these algorithms do not break the problem into smaller identical problems) becomes substantial.

You might ask, Where would one need to sort a million items?

Plenty of places. Imagine an electronic library that contains books, journals and articles in electronic format. Imagine a database of stars observed. Imagine a census of a country. Imagine a list of all species discovered. Now in the library, you might want to classify books by category, then by author and then by title. In the star database, you might be interested in properties of stars which lie within a given range of mass. With the census data, you might be simply interested in creating the lists of most common and least common names. In the species example, you might want to compare all organisms which use flagella.

I hope you get the idea across.

So sorting is used almost everyday in the modern world. And how do you sort? There are many ways to sort and I will try to encapsulate some of them in this blog. I will also try to provide the algorithm as well. But before we go on, we must try to understand the term time complexity given in the big O notation. The big O simply means that all the constants are removed and the highest order term is only considered. This implies that O(10n3 +5n2 +n) = O(n3 ).

Now, when you want to sort items, you would generally compare the value of different items with other items. If an algorithm does these comparisons, then it is called a comparison algorithm. Yes, non-comparison algorithms exist.

Now, it is understood that the least amount of effort for any sorting algorithm that uses comparisons is O(nlog2 n). This applies to the worst case scenario. Let’s see how we get this number. A single comparison can sort two objects. In general, we have a binary question as a comparison. If we assume that the inputs are random and if we take into account all permutations (n! for n objects), then the number of comparisons would have to distinguish between n! objects. Thus, we end up with the following inequality:

2f(n) ≥ n! => f(n) ≥ log2 n!

By using Stirling’s approximation, we get f(n) ≥ nlog2n.

In this way, we put a lower bound of the amount of work that a sorting algorithm that uses comparisons must use. Generally, the sorting algorithms which are non-recursive are of the order O(n2). Such algorithms are considered “bad”. But these are necessary to study since intuitively some of them are carried out by us. Before we begin, we must define the problem at hand. We have an unsorted list of integers and we wish to sort them in ascending order. For the sake of simplicity, the integers are stored in an array marked A(n). The index starts with 1 and ends with n, where n is the number of integers that would be sorted. In case we have integers or strings, or in case we have to sort them in descending order, the problem doesn’t change much on a conceptual level.

Selection Sort

Let’s start with the Selection Sort. In layman terms, it works like this:

First, find the smallest number in the array and swap it with that stored in the first place. Then, find the second smallest number in the array and swap it with that stored in the second place. Do the same with third, the fourth and so on.

Edit #1: You can check this algorithm at work here (A word of warning – this is a slow algorithm): AWESOME LINK!

You will notice that you only need to do it (n-1) times. Moreover, you do not need to check the initial numbers after they have been swapped or sorted. This saves slightly more time, but in the worst case scenario, the time complexity is O(n2), which means that the amount of time taken as n tends to large numbers (theoretically, infinity) is proportional to n2.

But selection sort is intuitive. Here’s the pseudo-code. I will use Fortran style since it is slightly more intuitive to understand than C. Moreover, the indexing starts from 1. In accordance with pseudo-code styles, I will skip the program and variable declarations. The variables are n, A(n), i and j – all integers. An integer type variable MIN will be used to track the smallest number. Another integer type variable TEMP will help in swapping the integers. According to Fortran style, /= is equivalent to asking “not equal”. Comments are marked by !, italics and underline.

do i = 1, n

! Assume that the ith element is the minimum.

MIN = i

! Starting from i helps us scan numbers after the sorted elements only.

do j = i, n

! If we find a number smaller, we “remember” its index.

if (A(j) < A(MIN)) then

MIN = j

end if

end do

! If the index for the smallest number is different from i, we swap the numbers.

if (MIN /= i) then

TEMP = A(MIN)

A(MIN) = A(i)

A(i) = TEMP

end if

end do

That’s it actually. One of the simplest programs to code. For small values of n, it is fast since the overhead time (time taken to initialize memory and other declarations) is low. But it is a quadratic algorithm since there are two loops, one of which runs over n, and the other over (n-i). These are nested loops so the number of comparisons is proportional to n2 plus or minus smaller order terms. The smaller order terms can be neglected when n becomes large.

Insertion Sort

Now, let’s consider the Insertion Sort. It is slightly faster than the selection sort and can be implemented very easily. It is very efficient for small data sets. However, for worst case scenarios, the time complexity is O(n2). Most people, while playing cards, use methods that are similar to the insertion sort algorithm. It works like this:

Compare the second item to the first item and swap if not in order. Compare the third item with the second item. If they are in order, go to the next item, otherwise swap. Do the same with the subsequent items.

Edit #2: The algorithm at work (Slightly better than selection sort): YET ANOTHER AWESOME LINK!

Like previously, you will only select (n-1) items. The number of comparisons can reduce heavily because you are trying to ‘insert’ an item in place. Also remember that the list above (or before) the currently selected item is already sorted. It is in this way that the insertion sort gains speed over the selection sort. In fact, in the best case scenario, i.e., when the list is already sorted, you will only do n comparisons. Selection sort will be O(n2) even in the best case. Below is the pseudocode for insertion sort in Fortran style. The variables are exactly the same as used in selection sort with the exception that we don’t need MIN anymore.

do i = 2, n

! We first compare with (i-1)th element. The TEMP variable holds the ‘selected’ item.

j = i – 1

TEMP = A(i)

! We keep moving the selected item up in the list till we either reach the first item or till we encounter a smaller item.

do while (j >= 1 .AND. A(j) > TEMP)

! We keep moving the items to proper place.

A(j+1) = A(j)

j = j – 1

end do

! We found the proper place. Now we store the selected item in its place.

A(j+1) = TEMP

end do

That’s it! If the list is already sorted, the nested DO-loop is not accessed at all, saving loads of time on sorted or partially sorted lists as compared to selection sort. However, for worst case it is still quadratic in time complexity.

There is another algorithm called the Bubble Sort. But I won’t discuss about it. It lies in the quadratic time complexity as well.

Edit #3: In case you really want to see bubble sort at work, here’s another link: AWESOME!

In the next post, I will talk about much faster sorting algorithms. These will typically be recursive in nature, thereby enabling us to divide the problem into smaller identical problems. As the division is done in two parts at each step, we will end up with the ideal limit for sorting algorithms that use comparisons – O(nlog2n).

Until next time!

068097116097032115111114116101100033033033

When (0.1 + 0.2) does not equal 0.3!

072101108108111032087111114108100033

Imagine this scenario. Person A goes to a bank and opens an account. Then he deposits 10 euros and 10 cents in the newly created account. Excellent. The next day, his friend transfers him 10 euros and 20 cents. His friend calls in and asks him to check if the money reached him. He marches on to the bank and asks for a bank statement. The bank statement shows that he has 30 euros and 31 cents.

Should he call his friend and return him one cent? Did his bank just pay him one cent as an interest of some sort?

You see, he knows everything about floating point and so he wonders if the bank made a mess up. So that’s the topic of today’s post – Floating Point. This is not an advanced topic and I believe that every person should be taught about it if he or she is interested in programming. Most courses just give an overview of the floating point representation but never talk about the consequences of the floating point.

Let’s start with our very own decimal system. It consists of 10 numbers 0 to 9 and any number above 9, would require more positions generally to the left. That is how we get numbers like 0, 1, 8, 9, 10, 23, and 8349. Then we have fractional numbers. Numbers that lie between others. These numbers use a dot (.) or a comma (,), depending on which country you are from and in which language you are writing. These numbers also require additional positions and they are filled to the right. That is how we get numbers like 0.1, 0.2, 0.002, 0.2354 and 0.43987.

The floating point representation is a way to represent small numbers, generally numbers like 0.00003432423 or 0.000001239. It is called the scientific notation and expresses the above numbers as 3.432423 x 10-5 or 1.239 x 10-6. This is a clever way to write a number with many zeros in a compact manner.

Computers sadly do not speak decimal, except through interfaces (the screen and keyboard). They operate in binary. This makes things tricky. Because not all decimal fractions can be expressed easily in binary. Why? Because binary is in base 2, which means you can use only 2 numbers, 0 and 1. This causes many fractional numbers to be recurring in binary. Take for example, 0.1 (yes, 0.1 euros is 10 cents). Let’s see the binary representation of it.

(0.1)10 = (0.0001100)2. The underline indicates that these numbers are recurring. And what about 0.2? Yes, that’s just (0.1)10 = (0.001100)2. In case you noticed, the numbers just shifted to the left. That’s what you get in binary when you multiply by 2.

Now computers do not understand recurring numbers. That’s the problem with 0.1+0.2. You see, computers work on fixed memory spaces. As a result, they have to terminate the recurring numbers. Moreover, computers do not store the binary numbers as shown above. In fact, they store them in a scientific notation. The IEEE convention for storing single precision real numbers in 32 bits indicates that one bit is stored for the sign (+ or -), 8 bits for the exponent and 23 bits for the fraction. Thus, our 0.1 will be stored in binary as:

0    01111011    10011001100110011001101

This is not good news. Because the whole recurring thing is terminated. Moreover, when you convert this back into decimal, you get

0.100000001490116119384765625

This is greater than 0.1. Let’s look at 0.2. According to the IEEE convention for single precision numbers (32 bits), 0.2 is stored in binary as:

0    01111100    10011001100110011001101

Yes, the numbers in the fractional part are the same, and the numbers in the exponent are different. Converting this into decimal gives us:

0.20000000298023223876953125

This is also greater than 0.2. To add them, we must have the same exponential term which is usually the bigger number. When we write 0.1 in that form, we get

0    01111100    11001100110011001100110

An extra one at the beginning of the fractional part because we divided by two (right shift) and moved the 1 present to the left of the decimal to the right. When we add the significand (23 bits) for 0.1 and 0.2 we get an extra bit which again must be rounded off. When we do that, we get

0    01111101   00110011001100110011010

Here, I didn’t use more than 23 bits for the addition. If we convert our answer, we get,

0.300000011920328955078125 > 0.3 (However, 0.3 is stored exactly like obtained. This means that 0.3 will always be stored in the computer as a number greater than 0.3.)

Usually for currency, we take only 2 decimal places. If we approximate it to the appropriate significant digits, we indeed get 0.30. But person A thinks, what if the bank has messed up here? Instead of rounding off nicely, what if the bank used an algorithm that while approximating to 2 significant digits, gives a value of 0.31 instead of 0.30?

Actually, thanks to the large number of bits (23 bits for the significand), this thankfully does not happen. Moreover, the error happens at the eighth significant digit in decimal notation. For most practical purposes, this can be safely ignored. However, things can get tricky when you perform operations involving large numbers (large exponents) and very small numbers (small exponents). As an example, consider currency conversion between a currency which is large (like Yen) and a currency which is very small (like the British Pound). In such cases, the errors might sneak in to the resulting value. But imagine currency conversions like the Zimbabwean dollar or the Vietnamese dong, where you have denominations of 1,000,000 or higher. This page (link) contains some of the most ‘valueless’ currencies today and compares them with the US dollar (apparently a currency having a smaller denomination has a higher ‘value’). It is easy to see with the help of those examples, that if you normalize values, you can lose precision. Resolving this issue can be done in only one way – increase the number of bits. Use double precision, triple precision or quadruple precision. There is no way around this.

Since this blog tries to consolidate different fields and domains where programming can be applied, let’s try to think of other fields where the floating point representation might cause problems. One of them is physics. You have big numbers related to the sizes of the planets or the stars or the observable universe AND you have very tiny numbers related to the sizes of the molecules or the atoms or the even tinier elementary particles. While manipulating these numbers, one has to be very careful to not forget about the floating point errors.

Now that we have understood how the floating point works, let’s get back to the question earlier. Should person A call his friend and return him one cent? Did his bank just pay him one cent as an interest of some sort?

The definite answer to the above question is that he does not owe his friend any cent. If the bank messed up with the floating point (or with a currency conversion), the money was created out of thin air. If the bank didn’t mess up, he was paid one cent as an interest of some sort.

Either way, he is one cent richer.

071111111100045098121101032087111114108100033