```python
import numpy as np
from sympy import prime
```

### Exercise 1: Comparing the Arithmetic, Harmonic, and Geometric Means

What we usually understand as the average of a sequence of numbers $x_1,x_2,\ldots,x_n$ is also called arithmetic mean, and defined as: $$\frac1n\sum_{k=1}^nx_k$$. There exist, however, other types of means: the harmonic mean $$\frac n{\sum\limits_{k=1}^n\dfrac1{x_k}}$$ and the geometric mean $$\sqrt[n]{\prod_{k=1}^nx_k}$$.

#### Questions:
- Implement three functions that define the arithmetic, harmonic and geometric mean respectively. (Hint: Remember that you can access individual entries of a list `lst` by their index `i` by writing `lst[i]`.))

- With numerical experimentation (i.e. by plugging in a few sequences of numbers), compare the sizes of these three means. Which one seems biggest, which one smallest?


```python
def arithmetic_mean(numbers):
    total = sum(numbers)
    n = len(numbers)
    return total/n

def harmonic_mean(numbers):
    reciprocal_sum = sum([1/number for number in numbers])
    n = len(numbers)
    return n/reciprocal_sum

def geometric_mean(numbers):
    product = 1
    for number in numbers:
        product *= number
    n = len(numbers)
    return pow(product, 1/n)
```


```python
numbers = [1,5,2,10]
print("Arithmetic Mean : " + str(arithmetic_mean(numbers)))
print("Harmonic Mean : " + str(harmonic_mean(numbers)))
print("Geometric Mean : " + str(geometric_mean(numbers)))
```

    Arithmetic Mean : 4.5
    Harmonic Mean : 2.2222222222222223
    Geometric Mean : 3.1622776601683795


### Conclusion : Arithmetic mean is the biggest followed by the geometric mean and lastly the harmonic mean



### Exercise 2: Series Summation made Easier with Functions

Now that we have introduced functions, we can also simplify the code from previous weeks considerably. Consider for example your work on computing the Riemann Zeta function $$\zeta(s)=\sum_{n=1}^\infty\frac1{n^s}\;.$$ 
You already simplified your code considerably by using variable assignment, list comprehension, and the `sum()` function, writing code like the following to compute the $100$-th partial sum of the series for $\zeta(2)$.

#### Questions:
- Write a function called `zeta_sum` that gets two quantities $s$ (the value at which we want to evaluate the partial sum) and $N$ (the number of terms in the partial sum), and returns the value of the partial sum of the Riemann zeta function.
- Write a similar function for the partial product of the Riemann zeta function, using the function name `zeta_prod`.
- Compare `zeta_prod` and `zeta_sum`. 
- Compute the $100$th, $200$th and $400$th partial products for $s=2.5$, and compare against the exact values of $\zeta(s)$.


```python
def zeta_sum(s, N):
    if N < 1:
        return None
    sum = 0
    for i in range(1, N+1):
        sum = sum + 1/(i**s)
    return sum

def zeta_prod(s, N):
    if N < 1:
        return None
    product = 1
    for i in range(1, N+1):
        product = product * 1/(1-((prime(i))**(-s)))
    return product

```


```python
print("Zeta sum using s = {} and n = {} is : {}".format(str(2.5), str(100), str(zeta_sum(2.5, 100))))
print("Zeta prod using s = {} and n = {} is : {}".format(str(2.5), str(100), str(zeta_prod(2.5, 100))))
```

    Zeta sum using s = 2.5 and n = 100 is : 1.3408255697514637
    Zeta prod using s = 2.5 and n = 100 is : 1.3414770673969854


#### Conclusion : zeta_prod is larger than zeta_sum


```python
import scipy.special
zeta = scipy.special.zeta

actual_value = zeta(2.5)

print("Zeta sum using s = {} and n = {} is : {}".format(str(2.5), str(100), str(zeta_sum(2.5, 100))))
print("Zeta prod using s = {} and n = {} is : {}".format(str(2.5), str(100), str(zeta_prod(2.5, 100))))
print("Actual zeta value : " + str(actual_value))


print("Zeta sum using s = {} and n = {} is : {}".format(str(2.5), str(200), str(zeta_sum(2.5, 200))))
print("Zeta prod using s = {} and n = {} is : {}".format(str(2.5), str(200), str(zeta_prod(2.5, 200))))
print("Actual zeta value : " + str(actual_value))


print("Zeta sum using s = {} and n = {} is : {}".format(str(2.5), str(400), str(zeta_sum(2.5, 400))))
print("Zeta prod using s = {} and n = {} is : {}".format(str(2.5), str(400), str(zeta_prod(2.5, 400))))
print("Actual zeta value : " + str(actual_value))
```

    Zeta sum using s = 2.5 and n = 100 is : 1.3408255697514637
    Zeta prod using s = 2.5 and n = 100 is : 1.3414770673969854
    Actual zeta value : 1.3414872572509173
    Zeta sum using s = 2.5 and n = 200 is : 1.3412524370325856
    Zeta prod using s = 2.5 and n = 200 is : 1.3414845958264436
    Actual zeta value : 1.3414872572509173
    Zeta sum using s = 2.5 and n = 400 is : 1.3414040800048226
    Zeta prod using s = 2.5 and n = 400 is : 1.341486541451443
    Actual zeta value : 1.3414872572509173


#### Conclusion : We see that values computed by zeta_prod are closer to the actual value compared to zeta_sum and moreover, increasing the value of N for both zeta_sum and zeta_prod resulted in higher accuracy

### Exercise 3: Zipping and Unzipping.

You have encountered the built-in function `zip()` which takes two lists of the same length and returns a list of paired tuples.


```python
list(zip([1,2,3,4],['a','b','c','d']))
```




    [(1, 'a'), (2, 'b'), (3, 'c'), (4, 'd')]



#### Exercise 3.a: Write your own version of this function (without using the built-in one).  Test your function on `[1,2,3,4]` and `['a','b','c','d']`. 


```python
def manual_zip(list1, list2):
    if len(list1) != len(list2):
        return None
    zipped_list = []
    for i in range(0, len(list1)):
        zipped_list.append((list1[i], list2[i]))
    return zipped_list

print(manual_zip([1,2,3,4], ['a','b','c','d']))
        
```

    [(1, 'a'), (2, 'b'), (3, 'c'), (4, 'd')]


If the sequences are of different length, you will likely get a  `list index out of range` error.

#### Exercise 3.b: Test your function also on `[1,2,3]` and `['a','b','c','d']`, and also on `[1,2,3,4]` and `['a','b','c']`.


```python
print(manual_zip([1,2,3], ['a','b','c','d']))
print(manual_zip([1,2,3,4], ['a','b','c']))
```

    None
    None


#### Exercise 3.c: Change your function such that it does not produce an error, by restricting the range over which you zip, and test all three cases.


```python
def manual_zip_with_restricted_range(list1, list2):
    restricted_range = min(len(list1), len(list2))
    zipped_list = []
    for i in range(0, restricted_range):
        zipped_list.append((list1[i], list2[i]))
    return zipped_list

print(manual_zip_with_restricted_range([1,2,3,4], ['a','b','c','d']))
print(manual_zip_with_restricted_range([1,2,3], ['a','b','c','d']))
print(manual_zip_with_restricted_range([1,2,3,4], ['a','b','c']))
```

    [(1, 'a'), (2, 'b'), (3, 'c'), (4, 'd')]
    [(1, 'a'), (2, 'b'), (3, 'c')]
    [(1, 'a'), (2, 'b'), (3, 'c')]


### Exercise 4: Approximating the derivative

The derivative of a function $f$ at a point $x$ is defined as
$$f'(x)=\lim_{h\to0}\frac{f(x+h)-f(x)}h\;.$$
We want to numerically study the quotient on the right-hand side as $h$ tends to zero.

#### Questions:
- Implement the right hand side of the equation above as a function `g(f,x,h)`, where $f$ is another function (the one which you want to compute the numerical derivative from).
- Using lambda functions (or otherwise), evaluate the numerical derivative of $f(x)=x^2$ at $x=1$, by using a numerical step $h=0.1$. Compare this result with the actual value $f'(1)=2$.
- Explore how the numerical derivative converges when $h$ shrinks, by taking $h=10^{-n}$ with $n=1,2,\ldots,18$. Explain what you observe. 


```python
def derivative(f, x, h):
    return (f(x+h) - f(x))/h

def f(x):
    return x*x

for n in range(1, 19):
    h = pow(10, -n)
    print("Computed value of f'(x) using x = {}, h = {} is : {}".format(str(1), str(h), str(derivative(f, 1, h))))
    
print("\nActual value of f'(x) using x = {} is : {}".format(str(1), str(2)))
```

    Computed value of f'(x) using x = 1, h = 0.1 is : 2.100000000000002
    Computed value of f'(x) using x = 1, h = 0.01 is : 2.0100000000000007
    Computed value of f'(x) using x = 1, h = 0.001 is : 2.0009999999996975
    Computed value of f'(x) using x = 1, h = 0.0001 is : 2.000099999999172
    Computed value of f'(x) using x = 1, h = 1e-05 is : 2.00001000001393
    Computed value of f'(x) using x = 1, h = 1e-06 is : 2.0000009999243673
    Computed value of f'(x) using x = 1, h = 1e-07 is : 2.0000001010878066
    Computed value of f'(x) using x = 1, h = 1e-08 is : 1.999999987845058
    Computed value of f'(x) using x = 1, h = 1e-09 is : 2.000000165480742
    Computed value of f'(x) using x = 1, h = 1e-10 is : 2.000000165480742
    Computed value of f'(x) using x = 1, h = 1e-11 is : 2.000000165480742
    Computed value of f'(x) using x = 1, h = 1e-12 is : 2.000177801164682
    Computed value of f'(x) using x = 1, h = 1e-13 is : 1.9984014443252818
    Computed value of f'(x) using x = 1, h = 1e-14 is : 1.9984014443252818
    Computed value of f'(x) using x = 1, h = 1e-15 is : 2.220446049250313
    Computed value of f'(x) using x = 1, h = 1e-16 is : 0.0
    Computed value of f'(x) using x = 1, h = 1e-17 is : 0.0
    Computed value of f'(x) using x = 1, h = 1e-18 is : 0.0
    
    Actual value of f'(x) using x = 1 is : 2


#### Conclusion : We see that for values of  "n" upto n=7, the trend is consistent in the sense that the manually computed value of the derivative becomes closer and closer to the actual value of the derivative, however, for n>7, the trend breaks down and the manually computed value oscillates around 2 and for values n>15, the manually computed value is zero altogether. This is likely due to internal limitations of the python implementation of floating point arithmetic causing overflows/underflows somewhere in the bytecode.

### Exercise 5: Roots of a cubic polynomial



The discriminant of a cubic polynomial $p(x)=ax^3+bx^2+cx+d$ is
$$\Delta=b^2c^2-4ac^3-4b^3d-27a^2d^2+18abcd\;.$$
The discriminant gives us information about the roots of the polynomial $p(x)$:

- if $\Delta>0$, then $p(x)$ has 3 distinct real roots,
- if $\Delta<0$, then $p(x)$ has 2 distinct complex roots and 1 real root,
- if $\Delta=0$, then $p(x)$ has at least 2 (real or complex) roots which are the same.

Represent a cubic polynomial $p(x)=ax^3+bx^2+cx+d$ by the coefficients `a,b,c,d` of numbers. For example, the polynomial $p(x)=x^3-x+1$ is represented as `1,0,-1,1`.

#### Questions:
- Write a function called `cubic_roots()` that returns information on the roots of a given polynomial. Test this function on the polynomials $x^3$, $x(1-x)^2$, $x^3-x+1$ and $x^3-x$. 
- Write a function called `has_three_real_roots()` which takes as input the four coefficients of a cubic polynomial in decreasing order, and returns `True` if $p(x)$ has 3 real distinct roots and `False` otherwise. Test this function on the same polynomials as before.







```python
def cubic_roots(a, b, c, d):
    discriminant = (b**2)*(c**2) - 4*a*(c**3) - 4*(b**3)*d - 27*(a**2)*(d**2) + 18*a*b*c*d
    if discriminant > 0:
        return ("\np(x) has 3 distinct real roots.")
    if discriminant < 0:
        return ("\np(x) has 2 distinct complex roots and 1 real root.")
    else:
        return ("\np(x) has at least 2 (real or complex) roots which are the same.")
```


```python
print(cubic_roots(1, 0, 0, 0))
print(cubic_roots(1, -2, 1, 0))
print(cubic_roots(1, 0, -1, 1))
print(cubic_roots(1, 0, -1, 0))
```

    
    p(x) has at least 2 (real or complex) roots which are the same.
    
    p(x) has at least 2 (real or complex) roots which are the same.
    
    p(x) has 2 distinct complex roots and 1 real root.
    
    p(x) has 3 distinct real roots.






```python
def has_three_real_roots(a, b, c, d):
    discriminant = (b**2)*(c**2) - 4*a*(c**3) - 4*(b**3)*d - 27*(a**2)*(d**2) + 18*a*b*c*d
    if discriminant > 0:
        return True
    else:
        return False

print(has_three_real_roots(1, 0, 0, 0))
print(has_three_real_roots(1, -2, 1, 0))
print(has_three_real_roots(1, 0, -1, 1))
print(has_three_real_roots(1, 0, -1, 0))
```

    False
    False
    False
    True



```python

```

Once you are done, save the jupyter notebook and submit it to QMPLUS under Lab Report Week 4.


```python
print([(k**3 + 7*k - 2) for k in range(2, 11)])
```

    [20, 46, 90, 158, 256, 390, 566, 790, 1068]



```python

```

2af21680dcedba7a504b8732b118713a0129c636
ipynb
Jupyter Notebook
math-in-python/Computing Approximate Derivatives and Roots.ipynb
kapuragam/everyday-python
61603eb31fceabea0f1d749274702dc32dc5bcae
MIT
math-in-python/Computing Approximate Derivatives and Roots.ipynb
kapuragam/everyday-python
61603eb31fceabea0f1d749274702dc32dc5bcae
MIT
math-in-python/Computing Approximate Derivatives and Roots.ipynb
kapuragam/everyday-python
61603eb31fceabea0f1d749274702dc32dc5bcae
MIT
Qwen/Qwen-72B
1. YES
2. YES