Week 3#
Introduction to Object-Oriented Programming#
These weeks exercises starts you working with classes. If you want a gentler introduction, exercises for Chapter 7 in Langtangen is recommended.
Exercise 1 — Quadratic functions#
In this exercise, we will build on the example given in the lectures of implementing 2nd degree polynomials as objects of a custom defined class. A general 2nd degree polynomial, aka, quadratic function, can be written as:
where the coefficients, \(a_2\), \(a_1\), and \(a_0\) uniquely defines the polynomial.
Exercise 1a) Defining the Quadratic class#
Create a class, Quadratic, that represents a general 2nd degree polynomial. Define the following methods:
A constructor (
__init__)A call method (
__call__) The constructor should take in the three coefficients in order:a_2,a_1, anda_0, and the call method should take the free variablex.
You class should be able to handle the following test script:
f = Quadratic(1, -2, 1)
x = np.linspace(-5, 5, 101)
plt.plot(x, f(x))
plt.show()
Implement your solution here:
Use this to test your implementation:
def test_Quadratic():
f = Quadratic(1, -2, 1)
assert abs(f(-1) - 4) < 1e-8
assert abs(f(0) - 1) < 1e-8
assert abs(f(1) - 0) < 1e-8
# test_Quadratic() # uncomment this line when testing
Exercise 1b) Pretty printing#
Extend your Quadratic class with a string special method (__str__) so that you can print a Polynomial object and get the polynomial written out on a readable form. Test by creating a polynomial object and printing it out.
Exercise 1c) Adding together polynomials#
Adding together two general quadratic functions:
gives a new quadratic function:
Implement this functionality using the addition special method (__add__). This method should return a new Quadratic-object, without changing the two that are added together. Your new class should be able to handle the following test script:
f = Quadratic(1, -2, 1)
g = Quadratic(-1, 6, -3)
h = f + g
print(h)
x = np.linspace(-5, 5, 101)
plt.plot(x, h(x))
plt.show()
(Because \(a_2 + b_2 = 0\), the resulting plot should be a straight line.)
Implement your solution here:
Use this to test your implementation:
def test_Quadratic_add():
f = Quadratic(1, -2, 1)
g = Quadratic(-1, 6, -3)
h = f + g
a2, a1, a0 = h.coeffs
assert a2 == 0
assert a1 == 4
assert a0 == -2
# test_Quadratic_add() # uncomment this line when testing
Exercise 1d) Finding the roots#
The roots of a general quadratic function,
are given by the quadratic formula $\({\displaystyle x={\frac {-b\pm {\sqrt {b^{2}-4ac}}}{2a}}.}\)$
Extend your Quadratic function with a method .roots() that finds and returns the real roots of the function (ignore the imaginary ones). Return the result as a tuple with 0, 1, or 2 elements.
Test your method on the three polynomials:
\(2x^2 -2x + 2\)
\(x^2 - 2x + 1\)
\(x^2 -3x + 2\)
Implement your solution here:
Use this to test your implementation:
def test_Quadratic_root():
f1 = Quadratic(2, -2, 2)
f2 = Quadratic(1, -2, 1)
f3 = Quadratic(1, -3, 2)
assert f1.roots() == ()
assert abs(f2.roots()[0] - 1) < 1e-8
assert abs(f3.roots()[0] - 1) < 1e-8 and abs(f3.roots()[1] - 2) < 1e-8
# test_Quadratic_root() # uncomment this line when testing
Exercise 1e) Finding the intersection of two quadratic functions#
Extend your class with a method that finds and returns the intersection points (if any) between two Quadratic-objects. It should work as follows:
f = Quadratic(1, -2, 1)
g = Quadratic(2, 3, -2)
print(f.intersect(g))
Hint: The intersections are all points solving \(f(x) = g(x)\), which can be written as \((f-g)(x) = 0\).
Test your solution by plotting the two functions and their intersections:
Exercise 2 — A class for general polynomials#
We now turn to looking at general polynomials of degree \(n\). These can be written as
or more compactly as $\({\displaystyle \sum _{k=0}^{n}a_{k}x^{k}}.\)$
We want to make a class that represents such a polynomial, and can take any number of coefficients in. The constructor of such a class could for example take in a list of the coefficients: [a0, a1, ..., aN]. However, this list will always have to be of length \(N\), and say we want to specify the polynomial \(x^{1000} + 1\), it is highly inefficient to pass in such a long list, as most coefficients are actually 0.
A better approach is to use a dictionary, where we use the index as the key and the coefficient as the value. Doing this, we can then specify only the non-zero coefficients, and simply skip those that are 0. So defining \(x^{1000} + 1\) would simply be: Polynomial({0: 1, 1000: 1}).
Exercise 2a) Defining the Polynomial class#
Define the Polynomial class with the following methods
A constructor (
__init__) that takes in the coefficients of the polynomial as a dictionaryA call method (
__call__) that computes f(x) for a given xA string method (
__str__) for informative printing of the polynomial
Your class should be able to handle the following test script
coeffs = {0: 1, 5:-1, 10:1}
f = Polynomial(coeffs)
print(f)
x = linspace(-1, 1, 101)
plt.plot(x, f(x))
plt.show()
Implement your solution here:
Exercise 2b): Adding general polynomials together#
We now want to be able to add together two general polynomial objects, which should produce a new general polynomial object. Mathematically, this is just an extension of the 2nd degree polynomial case which we saw in exercise (1). If we have
the sum will be defined by
Thus, if we add together two polynomials of degree \(m\) and \(n\), then the sum will have degree \(\max(m, n)\), i.e., the largest of the two.
Extend your class to add this functionality using the addition special method (__add__).
The class should handle the following test case:
f = Polynomial({0:1, 5:-7, 10:1})
g = Polynomial({5:7, 10:1, 15:-3})
print(f+g)
Which should produce the output: \(-3x^{15} + 2x^{10} + 1\)
Implement your solution here:
Hint: You will need to create a new coefficient dictionary for the new polynomial and add in the coefficients from the two polynomials. This can be slightly tricky getting the keys right. Here collections.defaultdict can be useful, but it isn’t necessary.
Exercise 2c) Defining a AddableDictionary class#
The previous exercise would have been a lot simpler, if we could simply add two dictionary objects together as follows:
a = {0: 2, 1: 3, 2: 4}
b = {0: -1, 1:3, 2: 3, 3: 2}
c = a + b
However, if you try to do this, you get an exception:
TypeError: unsupported operand type(s) for +: 'dict' and 'dict'
This means that there is no addition special method defined for dictionaries. However, we can extend the normal dictionary class to include this by adding a special method as follows
class AddableDict(dict):
def __add__(self, other):
...
Add the necessary code, so that our new AddableDict class can add two dictionaries together as follows:
a = AddableDict({0: 2, 1: 3, 2: 4})
b = AddableDict({0: -1, 1:3, 2: 3, 3: 2})
print(a + b)
And give the ouput: {0: 1, 1: 6, 2: 7, 3: 2}.
Implement your solution here:
Use this to test your implementation:
def test_AddableDict():
a = AddableDict({0: 2, 1: 3, 2: 4})
b = AddableDict({0: -1, 1:3, 2: 3, 3: 2})
c = a + b
assert c[0] == 1
assert c[1] == 6
assert c[2] == 7
assert c[3] == 2
# test_AddableDict() # uncomment this line when testing
Having made the AddableDict class, go back and change the Polynomial constructor, so that even if the user sends in the coefficients as a normal dictionary, it is converted to an AddableDict inside the Polynomial. Having done this, rewrite Polynomial.__add__, which should be trivial.
Exercise 2d) Derivative of a polynomial#
It is also the case that the derivative of a polynomial is a polynomial, if we have
then we get
which can be written as
where \(b_{k} = (k+1)a_{k+1}\).
Implement a method, derivative, that returns the function \(f'(x)\) as a new Polynomial object. Test your function by finding the derivative of
Implement your solution here:
Use this to test your implementation:
def test_derivative():
f = Polynomial({10:1, 6:-3, 2:2, 0:1})
f_deriv = f.derivative()
assert f_deriv.coeffs == {9:10, 5:-18, 1:4}
# test_derivative() # uncomment this line when testing
Exercise 2e) Multiplying polynomials#
It is also the case that the product of two polynomials form a new polynomial. If we again define
then the product is given by
Implement this functionality using the multiplication special method (__mul__). To acomplish this, you will need two nested for-loops over the coefficient dictionaries.
Test your implementation with the code block
f = Polynomial({2: 4, 1: 1})
g = Polynomial({3: 3, 0: 1})
print(f*g)
Which should give the output:
Implement your solution here:
Use this to test your implementation:
def test_Polynomial_mul():
f = Polynomial({2: 4, 1: 1})
g = Polynomial({3: 3, 0: 1})
h = f*g
assert h.coeffs == {5:12, 4:3, 2:4, 1:1}
# test_Polynomial_mul() # uncomment this line when testing
Exercise 3 - Quantum Harmonic Oscillator in One Dimension#
The quantum harmonic oscillator wave function, \(\psi_n (x)\) is a solution to the time-independent Scrödinger equation
where \(\hat{H}\) is the quantum harmonic oscillator hamiltonian and \(E_n\) is the energy at the \(n\)-th level (\(n\) must be a positive integer, \(n = 0, 1, 2,\)…).
In this exercise you will implement a class HOWF which represents a quantum harmonic oscillator wave function in one dimension.
Exercise 3a) Hermite polynomials#
To implement the harmonic oscillator wave function, a Hermite Polynomial is needed. The Hermite Polynomials are given by the recursive formula
for \(n = 2, 3, 4, ...\), and initial conditions are
Make a class HOWF. Start by implementing the constructor and the recursive private method _compute_Hermite(self, n).
The method _compute_Hermite(self, n) should return the hermite polynomial as an object of the class Polynomial, which you implemented in the previous exercise. The method should be implemented using recursion.
*Hint: You will easily get error messages if you try to multiply numbers with your polynomial class. Remember that you can consider a constant as a polynomial of the zero-th degree!
The constructor must take the level \(n\) as a paramterer. The constructor should make a call the method _compute_Hermite(self, n) to obtain and store the hermite polynomial of the \(n\)-th order. Name the variable of the hermite polynomial H.
Implement your solution here:
def test_Hermite():
H0 = lambda x: 1
H1 = lambda x: 2*x
H2 = lambda x: 4*x*x - 2
H3 = lambda x: 8*x**3 - 12*x
H4 = lambda x: 16*x**4 - 48*x**2 + 12
H5 = lambda x: 32*x**5 - 160*x**3 + 120*x
H_table = [H0, H1, H2, H3, H4, H5]
tol = 1e-12
for n in range(0, 6):
H = HOWF(n).H
for x in [0, 0.5, 1.0/3, 1, 3/2, 2]:
expected = H_table[n](x)
computed = H(x)
msg = "The implemented Hermite Polynomial yields unexpected result for n = %d\
\n\tH(%.2f) = %.13g != %.13g" %(n, x, expected, computed)
assert abs(expected - computed) < tol, msg
# test_Hermite() # uncomment this line when testing
Exercise 3b) The quantum harmonic wave function#
In this exercise you will be asked to implement two functions. The expressions would include the mass of the particle \(m\), plancks constant \(\hbar\) and the angular frequency \(\omega\). When computing, it is an advantage to simplify the expressions as much as possible. Therefore we simplifiy the calculations by using dimensionless variables, defined by \(\chi = \sqrt{\frac{m \omega}{\hbar}} \cdot x\). This will make all the calculation dimensionless. Consequently, the dimensionless energy (in this case scale) is then given by \(\epsilon_n = \frac{2 E_n}{\hbar \omega}\).
The quantum harmonic wave function of the \(n\)-th energy level is then given by
where \(H_n\) is the Hermite polynimial of the \(n\)-th order. Let HOWF a callable class by implementing the formula above in the special method __call__(self, chi).
The energy of the \(n\)-th level is then given by
Implement a method in HOWF which calculates the energy of the wave function. Define this method as a property.
Implement your solution here:
from numpy import exp, pi
from math import factorial, sqrt
Exercise 3c) Visualising the wave functions#
You will now reproduce a quite iconic figure of the quantum harmonic oscillator wave function. Plot the wave functions in relative height to their energy, as in plot \(\psi_n(\chi) + \epsilon_n\) for \(n = 0, 1, 2, ..., N\) in the same plot.
Set \(N\) to be at least 15. A sanity check for your implementation of the wave functions is that the number of nodes of \(\psi_n\) should be \(n\). (That means that \(\psi_n\) should intersect \(\epsilon_n\) \(n\) times on the y-axis, since the wave function has been lifted by \(\epsilon_n\) on the \(y\)-axis). An appropriate range for \(\chi\) could be \(\chi \in [-8, 8]\).
Implement your solution here:
Exercise 4 — Fibonnaci with Memoization#
The Fibonnaci sequence is given by the numbers $\(1, 1, 2, 3, 5, 8, 13, 21, 34, \ldots\)\( it is defined by the recursive formula \)\(F_{i} = F_{i-1} + F_{i-2}.\)$
For this defintion to make any sense, the sequence has to start somewhere, so we define that
Exercise 4a) A Fibonacci function#
Write a recursive function fibonacci(n), that returns \(F_n\). A recursive function is a function that calls itself. To do this, you won’t need any explicit loops. You will however, need to include the base-case of \(F_0=0\) \(F_1=1\), otherwise the recursion would just continue forever.
Verify your function by writing out \(F_1, \ldots, F_{10}\) and comparing to the sequence above.
Implement your solution here:
Use this to test your implementation:
def test_fibonacci():
computed = []
for i in range(11):
computed.append(fibonacci(i))
assert computed == [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
# test_fibonacci() # uncomment this line when testing
Exercise 4b) Improving our algorithm#
We have now implemented a function that finds the \(n\)’th Fibonnaci number, which is great. However, it is very inefficient. Imagine for example we want to compute \(F_{100}\). The first call to the function, creates two new calls: F(99) and F(98). Both of these call the function twice again, so we have four function calls there. Each of those four lead to two new ones and so on all the way down to the base case of \(n=1\). This means we make something on the order of \(2^{100}\) function calls! This means the complexity of calulating \(F(n)\) grows exponentially as \(n\) grows, which is horrible. Try drawing a tree to represent the function calls for \(n=100\) and you quickly see the absurdity of the situation.
To improve our Fibonnaci algorithm, we will use a dynamic programming technique called memoization. While the name is a bit weird, the idea is fairly simple, we simply make our program remember the answers it has already computed, so that we won’t have to compute them again later.
Classes are perfect for making memoized functions, as we can use an internal dictionary to remember old solutions.
Fill in the skeleton code below so that the class computes Fibonacci numbers:
class Fibonacci:
def __init__(self):
self.memory = {0: 1, 1: 1}
def __call__(self, n):
"""
if n is in memory
- return it
if n is not in memory
- calcuate it recursively
- put it into memory
- return it
"""
pass
If you program your class correctly, you should see a huge speed up. This is because we now avoid repeating the same calculations uneccesarily. For example, if you try computing \(F(100)\) as follows:
fib = Fibonacci()
print(fib(100))
None
The memoized function will only need to compute the numbers \(F(2), \ldots F(100)\) once, meaning the number of function calls needed is now linear instead of exponential!
Implement your solution here:
Same test to make sure the code still works:
# test_fibonacci() # uncomment this line when testing
Exercise 4c) Maximum Recursion Depth#
Try finding an even larger number, like \(F_{100000}\). Does it work? It will probably fail due to a maximum recursion depth. Try to understand what this means and why it happens. Can you think of a work-around for this problem?
Hint: You need to build up the memoization dictionary from the ground up, and then it will be able to handle larger inputs. This can be done with a for-loop for example.
Implement your solution here:
Same test to make sure the code still works:
def test_fibonacci():
computed = []
for i in range(11):
computed.append(fibonacci(i))
assert computed == [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
# test_fibonacci() # uncomment this line when testing
Exercise 4d) A memoized factorial class#
The factorial is also defined recursively
with a base case of \(F(0) = 1\).
Repeat the process of the Fibonacci sequence and create a memoized class, Factorial, that computes \(N!\).
Implement your solution here:
Use this to test your implementation:
def test_Factorial():
f = Factorial()
computed = []
for i in range(6):
computed.append(f(i))
assert computed == [1, 1, 2, 6, 24, 120]
# test_Factorial() # uncomment this line when testing