# Lecture 2: Debugging, Control Statements, Functions
By Dr. Anders S. Christensen
`anders.christensen @ unibas.ch`

## Debugging Exercises:

First, we go through some common errors that you **will** encounter in this course!

I wrote the following bits of code, neither of them work!



```
# Error 1:

k = 0

# Sum all numbers from 0 to 9
for i in range(10):
  k = k + 1

 print(k)
```


```
# Error 2:

my_list = [1.0, 2.0, 3.0, 4.0, 5.0]

print("Length of list:", len(my_list))

print(my_list[1])
print(my_list[2])
print(my_list[3])
print(my_list[4])
print(my_list[5])
```


```
# Error 3:

y = [3.5, 2.1, -2.5, 1.2, 6.2, 5.2, -8.1]

# print y
print(y)

# print second entry in y
print(y(2))

print("Done!")
```


```
# Error 4:

for i in range(10)
  print(i)

print "Done!"
```

## Example 2.1: `if` and `else` statements



These statments can do things depending on a condition, for example:


```
a = 33
b = 200
if b > a:
  print("b is greater than a")
```

Or with branching:


```
a = 33
b = 20
if b > a:
  print("b is greater than a")
else:
  print("b is NOT greater than a")
```

Even more possibilites with `elif` (short for else+if):


```
a = 33
b = 20
if b > a:
  print("b is greater than a")
elif b == a:
  print("b is equal to a")
else:
  print("b is smaller than a")
```

### Other common checks:


```
a == b # Note two == signs!!
```


```
a >= b # Greater than or equal to
```


```
a > b # Greater than
```


```
a <= b # Smaller than or equal to
```


```
a < b # Smaller than
```


```
a != b # Not equal to
```

`True` and `False` are called "boolean"!

The tests above are from what is called "boolean algebra"


```
# Boolean as a variable

my_var = True
# var = False
print(my_var)

print(type(my_var))
```

What happens here?


```
if my_var:
  print("That was True")
```

Other keywords: `and`, `or`, `not`


```
a = 33
b = 20

if (a > b) and (a == 33):
  print("Complicated")
```

## Example 2.2: The `break` and `continue` statements

`Break` is used to break out of a `for` loop, for example when a condition is met:


```
# We only want to print up to i==5:

for i in range(100):

  if i > 5:
    break

  print(i)

```

The `continue` statment is used to skip only the remaining part of the current loop iteration:


```
# Print odd numbers:

for i in range(20):

  # print(i, i%2)

  # NOTE: % is the modulo operator
  if i % 2 == 0:
    continue

  print(i)
```

## Example 2.3: Defining Functions

Last time we used functions from NumPy to calcualte the exponential


```
import numpy as np

e = np.exp(1.0)

print(e)
```

We can define such functions ourself, for example
\begin{equation}
f\left(x\right) = x^{2}
\end{equation}


```
def square(x):
  f = x**2
  return f
```

Let's try our new function:


```
print(square(2))
```

### Demonstration: Define your own function!

Can you define a function such that?

\begin{equation}
f\left(x\right) = \sum_{n=0}^{10} x^{n}
\end{equation}


```
# Let's write something here!


```

#### Possible Solution:


```
def power_sum(x):

  f = 0

  for n in range(11):
    f = f + x**n
    
  return f

print(power_sum(11))
print(power_sum(12))
```

### **Note:**
A common error when defining functions is demonstrated here: What could be wrong??


```
def cubed(x):
  f = x**3
  
print(cubed(3))
```

## Code Demonstration: Heat Capacity as a function of the temperature + plotting

Same example as yesterday, but today we will look at a better code and end up with a plot.




The energy $\epsilon_n$ for a state of the harmonic oscillator is given by
\begin{equation}
\epsilon_n = n+\frac{1}{2},
\end{equation}
in units of $\hbar\omega$. We also assume that there are no degenerate states!

The probability for observing the state $n$ at a given temperature $\tau$ is given by:
\begin{equation}
p(n) = \frac{\exp{\left(-\epsilon_n / \tau \right)}}{Z},
\end{equation}

where $Z$ is the partition function (*Zustandssumme*)
\begin{equation}
Z = \sum_{n=0}^{\infty} \exp{\left(-\epsilon_n / \tau\right)}.
\end{equation}

The total *internal* energy $U$ is then given by:
\begin{equation}
U = \langle \epsilon \rangle = \sum_n \epsilon_n p(n)
\end{equation}

---

Another interesting property of this system is the *heat capactity* $C_V$ given by
\begin{align}
C_V &= \frac{\partial U}{\partial \tau} \biggl|_V \\
&= \frac{\langle \epsilon^2 \rangle - \langle \epsilon \rangle^2}{\tau^2}.
\end{align}


The expectation value $\langle \epsilon \rangle$ is given via the energy $U$:
\begin{align}
\langle \epsilon \rangle &= \sum_n \epsilon_n p(n) \\
&= \sum_n \left(n+\frac{1}{2}\right) p(n)
\end{align}

And $\langle \epsilon^2 \rangle$ is the second moment of the energy given by:
\begin{align}
\langle \epsilon^2 \rangle &= \sum_n \epsilon^2_n p(n) \\
&= \sum_n \left(n+\frac{1}{2}\right)^2 p(n)
\end{align}


#### Ugly code from last time:





```
# 7 temperatures
temperatures = [0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]
n_states = 1000

for tau in temperatures:

  # Calculate the partition function for tau
  Z = 0.0
  for n in range(n_states+1):
    Z = Z + np.exp(-(n + 0.5) / tau)

  # First moment of the energy
  U = 0.0

  # Second moment of the energy
  U2 = 0.0

  for n in range(n_states+1):
  
    epsilon = n + 0.5
    prob = np.exp(-epsilon/tau) / Z
  
    U = U + epsilon * prob
    U2 = U2 + epsilon**2 * prob
  
  Cv = (U2 - U**2) / tau

  print(tau, Cv)
```

### Example 1: Print the energy $\epsilon_n$ for each state up to $n=3$?

The energy $\epsilon_n$ for a state of the harmonic oscillator is given by:
\begin{equation}
\epsilon_n = n+\frac{1}{2}
\end{equation}

Last time:
```
for n in range(3+1):

  epsilon = n + 0.5

  print(n, epsilon)
```
Now, lets write a function for this:


```
# Energy for n'th state
def epsilon(n):
  e = n+0.5
  return e

print(epsilon(1))
```

    1.5


#### Example 2: What is the partition function $Z$ 

for $\tau$ = 2.0?

Recall that the partition function is given by:

\begin{equation}
Z(\tau) = \sum_{n=0}^{\infty} \exp{\left(-\frac{\epsilon_n}{\tau}\right)}
\end{equation}


Yesterday:
```
tau = 2.0
Z = 0.0

for n in range(10000):
  
  Z = Z + np.exp(-(n + 0.5) / tau)
```




```
# Import exponential function
import numpy as np

# Partition function at temperature Tau
def partition(tau):
  Z = 0

  for n in range(100):

      # Code from last time
      # Z = Z + np.exp(-(n + 0.5) / tau)
      
      Z = Z + np.exp(-epsilon(n) / tau)

  return Z


```


```
Z_2 = partition(2.0)

print(Z_2)
```

**Question:** When should we stop the loop??
*    Can we use `break`?

#### Example 3: What is the internal energy $U$ for $\tau$ = 2.0?

\begin{align}
\langle \epsilon \rangle &= \sum_n \epsilon_n p(n) \\
\langle \epsilon^2 \rangle &= \sum_n \epsilon_n^2 p(n)
\end{align}

Also recall that:

\begin{equation}
p(n, \tau) = \frac{\exp{\left(-\epsilon_n / \tau \right)}}{Z(\tau)},
\end{equation}



```
# Probability for n'th state at temperature Tau
def prob(n, tau):

  Z = partition(tau)

  p = np.exp(-epsilon(n) / tau) / Z

  return p


print( prob(0, 2.0) )
print( prob(1, 2.0) )
print( prob(4, 2.0) )
print( prob(100, 2.0) )
```


```
# Internal energy at temperature Tau
def internal_energy(tau):

  u = 0

  for n in range(100):
    u = u + epsilon(n) * prob(n, tau)
  
  return u


# Second moment of the energy at temperature Tau
def second_moment(tau):

  u2 = 0

  for n in range(100):
    u2 = u2 + epsilon(n)**2 * prob(n, tau)
  
  return u2


print( internal_energy(2.0))
print( second_moment(2.0))

```

### Example 4: The Heat capacity as a function of $\tau$

Definition of heat capacity $C_V$:
\begin{equation}
C_V = \frac{\langle \epsilon^2 \rangle - \langle \epsilon \rangle^2}{\tau^2}.
\end{equation}



```
def heat_capacity(tau):

  u = internal_energy(tau)
  u2 = second_moment(tau)

  cv = (u2 - u**2) / tau**2

  return cv

print(heat_capacity(2.0))
```

Finally print $C_V$ for every $\tau$:


```
temperatures = [0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0]

for tau in temperatures:

  cv = heat_capacity(tau)
  print(tau, cv)

```

## Plotting everything with Pyplot/Matplotlib

First, how do we save the results in a list?


```
temperatures = [0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2.0, 3.0]

# Initialize an empty list
heat_capacities = []

for tau in temperatures:

  # calculate Cv
  cv = heat_capacity(tau)

  # Add Cv to the list!
  heat_capacities.append(cv)


print(heat_capacities)

```

To plot a graph, we can use "pyplot" from matplotlib:


```
import matplotlib.pyplot as plt

# Construct the plot
plt.plot(temperatures, heat_capacities)

# Draw the plot
plt.show()
```

Matplotlib has tons of options for customization, for example labels, titles, etc.

Take a look at the matplotlib gallery: https://matplotlib.org/gallery.html


```
# Figure size in inches
plt.figure(figsize=(8, 8))

# Add labels to axes
plt.xlabel("Tau")
plt.ylabel("Heat Capacity")

# Set a descriptive title
plt.title("Quantum Harmonic Oscillator")

# Display a grid
plt.grid(True)

plt.plot(temperatures, heat_capacities)
plt.show()
```

404c8e1f891a23a4a5454f46b0d38244041850f4
ipynb
Jupyter Notebook
intro_lecture2_2020.ipynb
andersx/python-intro
8409c89da7dd9cea21e3702a0f0f47aae816eb58
CC0-1.0
2020-05-03T11:59:01.000Z
2021-11-15T12:33:39.000Z
intro_lecture2_2020.ipynb
andersx/python-intro
8409c89da7dd9cea21e3702a0f0f47aae816eb58
CC0-1.0
intro_lecture2_2020.ipynb
andersx/python-intro
8409c89da7dd9cea21e3702a0f0f47aae816eb58
CC0-1.0
2020-05-10T21:15:15.000Z
2021-12-05T15:13:54.000Z
Qwen/Qwen-72B
1. YES
2. YES