<jupyter_start><jupyter_text>Pure Python evaluation of vector normsGenerate a list of random floats of a given dimension (dim), and store its result in the variable `vec`.<jupyter_code># This is used for plots and numpy
%pylab inline 

import random as rand
dim = int(10)

vec = []
for i in range(dim):
    elem = rand.random()
    vec.append(elem)

print vec

# YOUR CODE HERE
#raise NotImplementedError()
from numpy.testing import *

assert_equal(type(vec), list)
assert_equal(len(vec), dim)
for ob in vec:
    assert_equal(type(ob), float)<jupyter_output><empty_output><jupyter_text>Write a function that evaluates the $l_p$ norm of a vector in $R^d$. We remind:$$\|v \|_{p} := \left(\sum_i (v_i)^p\right)^{1/p}$$the function should take as arguments a `list`, containing your $R^d$ vector, and a number `p` in the range $[1, \infty]$, indicating the exponent of the norm. **Note:** an infinite float number is given by `float("inf")`.Throw an assertion (look it up on google!) if the exponent is not in the range you expect.<jupyter_code>def p_norm(vector,p):
    assert(p >= 1)
    abs_vec = map(abs,vector)     # returns the abs of every elem of the list: max(abs(i) for i in ll) 
    if p != float("inf"):
        sum_tmp = sum( elem**float(p) for elem in abs_vec)
        sum_tmp = pow(sum_tmp, 1/float(p))
    else:
        sum_tmp = max(abs_vec)
    return sum_tmp

#print p_norm(vec,2)
#print p_norm(vec, "inf")

assert_equal(p_norm(range(10),1), 45.0)
assert_equal(p_norm([3,4], 2), 5.0)
assert_equal(p_norm([-1,-.5,.5], float("inf")), 1)
assert_raises(AssertionError, p_norm, [2,3], 0)
assert_raises(AssertionError, p_norm, [2,3], -1)<jupyter_output><empty_output><jupyter_text>Playing with condition numbersIn this exercise you will have to figure out what are the optimalvalues of the stepping interval when approximating derivatives usingthe finite difference method. See here_ for a short introduction onhow to run these programs on SISSA machines. 1. Finite differencesWrite a program to compute the finite difference (`FD`)approximation of the derivative of a function `f`, computed atpoint `x`, using a stepping of size `h`. Recall the definition ofapproximate derivative:$$FD(f,x,h) := \frac{f(x+h)-f(x)}{h}$$<jupyter_code>def FD(f, x, h):
    return (f(x + h) - f(x))/float(h)
assert_equal(FD(lambda x: x, 0, .125), 1.0)<jupyter_output><empty_output><jupyter_text>2. Compute FDEvaluate this function for the derivative of `sin(x)` evaluated at `x=1`, for values of `h` equal to `1e-i`, with `i=0,...,20`. Store the values of the finite differences in the list `fd1`.<jupyter_code>fd1 = []
for i in range(21):
    value_der = FD(sin, 1, pow(10,-i))
    fd1.append(value_der)

print fd1
    
    
assert_equal(len(fd1), 21)
expected = [0.067826442017785205, 0.49736375253538911, 0.53608598101186899, 0.5398814803603269, 0.54026023141862112, 0.54029809850586474, 0.54030188512133037, 0.54030226404044868, 0.54030229179602429, 0.54030235840940577, 0.54030224738710331, 0.54030113716407868, 0.54034554608506369, 0.53956838996782608, 0.53290705182007514, 0.55511151231257827, 0.0, 0.0, 0.0, 0.0, 0.0]
assert_almost_equal(fd1,expected,decimal=4)<jupyter_output><empty_output><jupyter_text>3. Error plotsPlot the error, defined as `abs(FD-cos(1.0))` where `FD` is your approximation, in `loglog` format and explain what you see. A good way to emphasize the result is to give the option `'-o'` to the plot command.<jupyter_code>error = []
for i in range(21):
    error_der = abs(FD(sin, 1, pow(10,-i)) - cos(1))
    error.append(error_der)

h = [pow(10,-i) for i in range(21)]

#figure(figsize=[10,10])
#_ = loglog(h, error)

figure(figsize=[5,5])
plt.loglog(h,error, "-o")
plt.show()<jupyter_output><empty_output><jupyter_text>We have three sources of error. FD_err : error of the numerical derivarive. It is proportional to h: when h decreases, it does the samedifference_err : error due to the difference between two numbers.If h is very small, we have less significant digits. It goes as 1/h. binary_err : h goes as like powers of 10, while it would be better that h go as power of 2At a certain point, around 10^16, the error saturates since we have reached the precision of the machine  4. Error plots base 2Repeate step 2 and 3 above, but using powers of `2` instead of powers of `10`, i.e., using `h` equal to `2**(-i)` for `i=1,...,60`. Do you see differences? How do you explain these differences? Shortly comment.  A good way to emphasize the result is to give the option `'-o'` to the plot command.<jupyter_code>error = []
for i in range(61):
    error_der = abs(FD(sin, 1, pow(2,-i)) - cos(1))
    error.append(error_der)

h = [pow(2,-i) for i in range(61)]

#figure(figsize=[10,10])
#_ = loglog(h, error)

figure(figsize=[5,5])
plt.loglog(h,error, "-o")
plt.show()<jupyter_output><empty_output><jupyter_text>The considerations are similar to the ones before, but now we do not have the binary error contribution  5. Central Finite DifferencesWrite a function that computes the central finite difference approximation (`CFD`), defined as $$CFD(f,x,h) := \frac{f(x+h)-f(x-h)}{2h}$$<jupyter_code>def CFD(f, x, h):
    return (f(x + h) - f(x - h))/(2*float(h))
assert_equal(CFD(lambda x: x**2, 0.0, .5), 0.0)
assert_equal(CFD(lambda x: x**2, 1.0, .5), 2.0)<jupyter_output><empty_output><jupyter_text>6. Error plots for CFDRepeat steps 2., 3. and 4. and explain what you see. What is the *order* of the approximation 1. and what is the order of the approximation 5.? What's the order of the cancellation errors?<jupyter_code>error = []
for i in range(21):
    error_der = abs(CFD(sin, 1, pow(10,-i)) - cos(1))
    error.append(error_der)

h = [pow(10,-i) for i in range(21)]

#figure(figsize=[10,10])
#_ = loglog(h, error)

figure(figsize=[5,5])
plt.loglog(h,error, "-o")
plt.show()
error = []
for i in range(61):
    error_der = abs(CFD(sin, 1, pow(2,-i)) - cos(1))
    error.append(error_der)

h = [pow(2,-i) for i in range(61)]

#figure(figsize=[10,10])
#_ = loglog(h, error)

figure(figsize=[5,5])
plt.loglog(h,error, "-o")
plt.show()<jupyter_output><empty_output><jupyter_text>In the CFD case we are evaluating f(x + h) and f(x - h) and subtracting them. It is as you take the FD derivative on the right and the FD derivative on the left.err_dx = c*h + O(h^2)err_sx = -c*h + O(h^2)Therefore the error goes as h^2.  NumpyNumpy provides a very powerful array container. The first line of this ipython notebook has imported all of numpy functionalities in your notebook, just as if you typed::    from numpy import *    Create a numpy array whith entries that range form 0 to 64. Use the correct numpy function to do so. Call it `x`.<jupyter_code>x = np.array([k for k in range(64)])
assert_equal(type(x), ndarray)
assert_equal(len(x), 64)
for i in xrange(64):
    assert_equal(x[i], float(i))<jupyter_output><empty_output><jupyter_text>Reshape the one dimensional array, to become a 4 rows 2 dimensional array, let numpy evaluate the correct number of culumns. Call it `y`.<jupyter_code>y = np.reshape(x,(4,16))
print y

assert_equal(shape(y), (4,16))<jupyter_output><empty_output><jupyter_text>Get the following *slices* of `y`:* All the rows and the first three colums. Name it `sl1`.* All the colums and the first three rows. Name it `sl2`.* Third to sixth (included) columns and all the rows. Name it `sl3`.* The last three columns and all the rows. Name it `sl4`.<jupyter_code>sl1 = y[0:4]
sl1 = sl1[:,0:3]

sl2 = y[:,0:]
sl2 = sl2[0:3]

sl3 = y[:,3:7]

sl4 = y[:,-3:]

print sl1
print sl2
print sl3
print sl4
assert_equal(sl1,[[0,1,2],[16,17,18],[32,33,34],[48,49,50]])
assert_equal(sl2,[[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15],[16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31],[32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47]])
assert_equal(sl3,[[3,4,5,6],[19,20,21,22],[35,36,37,38],[51,52,53,54]])
assert_equal(sl4,[[13,14,15],[29,30,31],[45,46,47],[61,62,63]])<jupyter_output><empty_output><jupyter_text>Now reshape the array, as if you wanted to feed it to a fortran routine. Call it `z`.<jupyter_code>z = np.reshape(x,(4,16), order = 'F')

print z
print y<jupyter_output>[[ 0  4  8 12 16 20 24 28 32 36 40 44 48 52 56 60]
 [ 1  5  9 13 17 21 25 29 33 37 41 45 49 53 57 61]
 [ 2  6 10 14 18 22 26 30 34 38 42 46 50 54 58 62]
 [ 3  7 11 15 19 23 27 31 35 39 43 47 51 55 59 63]]
[[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15]
 [16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31]
 [32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47]
 [48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63]]<jupyter_text>Comment on the result, what has changed with respect to `y`?  The difference is that in python elements are ordered row-wise, while in fortran they are ordered colum-wise. Set the fourth element of `x` to 666666, and print `x`, `y`, `z`. Comment on the result<jupyter_code>x[3] = 666666
print x
print y
print z<jupyter_output>[     0      1      2 666666      4      5      6      7      8      9
     10     11     12     13     14     15     16     17     18     19
     20     21     22     23     24     25     26     27     28     29
     30     31     32     33     34     35     36     37     38     39
     40     41     42     43     44     45     46     47     48     49
     50     51     52     53     54     55     56     57     58     59
     60     61     62     63]
[[     0      1      2 666666      4      5      6      7      8      9
      10     11     12     13     14     15]
 [    16     17     18     19     20     21     22     23     24     25
      26     27     28     29     30     31]
 [    32     33     34     35     36     37     38     39     40     41
      42     43     44     45     46     47]
 [    48     49     50     51     52     53     54     55     56     57
      58     59     60     61     62     63]]
[[     0      4      8     12     16     20     24     28     32     36
   [...]<jupyter_text>The element x[3] has been overwritten with a new value.  Arrays and MatricesDefine 2 arrays, `A` of dimensions (2,3) and `B` of dimension (3,4).* Perform the operation `C = A.dot(B)`. Comment the result, or the error you get.<jupyter_code>A = np.array([[i*j for i in range(1,4)] for j in range(1,3)])
B = np.array([[i*j for i in range(1,5)] for j in range(1,4)])
C = A.dot(B)
print A
print B
print C
assert_equal(A.shape,(2,3))
assert_equal(B.shape,(3,4))
assert_equal(C.shape,(2,4))
expected = sum(A[1,:]*B[:,2])
assert_equal(C[1,2],expected)<jupyter_output><empty_output><jupyter_text>YOUR ANSWER HERE * Perform the operation `C = A*(B)`. Comment the result, or the error you get.<jupyter_code>C = A*B<jupyter_output><empty_output><jupyter_text>YOUR ANSWER HERE  * Convert A and B, from arrays to matrices and perform `A*B`. Comment the result.<jupyter_code>A = np.asmatrix(A)
B = np.asmatrix(B)
C = A*B
print C
print type(C)


assert_equal(type(A),numpy.matrixlib.defmatrix.matrix)
assert_equal(type(B),numpy.matrixlib.defmatrix.matrix)
assert_equal(type(C),numpy.matrixlib.defmatrix.matrix)
assert_equal(A.shape,(2,3))
assert_equal(B.shape,(3,4))
assert_equal(C.shape,(2,4))
expected = sum(A[1,:]*B[:,2])
assert_equal(C[1,2],expected)<jupyter_output><empty_output><jupyter_text>YOUR ANSWER HERE  Playing with polynomialsThe polynomial `(1-x)^6` can be expanded to::          x^6 - 6*x^5 + 15*x^4 - 20*x^3 + 15*x^2 - 6*x + 1   The two forms above are equivalent from a mathematical point of   view, but may yield different results in a computer machine.   Compute and plot the values of this polynomial, using each of the   two forms, for 101 equally spaced points in the interval   `[0.995,1.005]`, i.e., with a spacing of 0.0001 (use linspace).      Can you explain this behavior?<jupyter_code>pol = lambda x: (1 - x)**6
pol_expanded = lambda x: x**6 - 6*x**5 + 15*x**4 - 20*x**3 + 15*x**2 - 6*x + 1

x = np.linspace(0.995,1.005,(1.005-0.995)/0.0001)

pol_value = pol(x)
pol_exp_value = pol_expanded(x)

plt.plot(x,pol_value, color = "green")
plt.plot(x,pol_exp_value, color = "blue")
plt.show()<jupyter_output><empty_output><jupyter_text>Differences between small numbers yield to big oscillations (errors), since the significant digits are few. **Playing with interpolation in python**1. Given a set of $n+1$ points $x_i$ as input (either a list of floats, or a numpy array of floats), construct a function `lagrange_basis(xi,i,x)` that returns the $i$-th Lagrangepolynomial associated to $x_i$, evaluated at $x$. The $i$-th Lagrange polynomial is defined as polynomial of degree $n$ such that $l_i(x_j) = \delta_{ij}$, where $\delta$ is one if $i == j$ and zero otherwise.   Recall the mathematical definition of the $l_i(x)$ polynomials:$$l_i(x) := \prod_{j=0, j\neq i}^{n} \frac{x-x_j}{x_i-x_j}$$<jupyter_code>def lagrange_basis(xi, i, x):   
    assert(i < len(xi) and i >= 0)
    li = np.product([(x - xi[j])/(xi[i]-xi[j]) for j in range(len(xi)) if j!=i], axis = 0)
    return li

x = linspace(0,1,5)
d = 3
xi = linspace(0,1,d)

print lagrange_basis(xi,1,x)
x = linspace(0,1,5)
d = 3
xi = linspace(0,1,d)

assert_equal(list(lagrange_basis(xi, 0, x)),[1.0, 0.375, -0.0, -0.125, 0.0])
assert_equal(list(lagrange_basis(xi, 1, x)),[0.0, 0.75, 1.0, 0.75, -0.0])
assert_equal(list(lagrange_basis(xi, 2, x)),[-0.0, -0.125, 0.0, 0.375, 1.0])

assert_raises(AssertionError, lagrange_basis, xi, -1, x)
assert_raises(AssertionError, lagrange_basis, xi, 10, x)<jupyter_output><empty_output><jupyter_text>Construct the function `lagrange_interpolation(xi,g)` that, given the set of interpolation points `xi` and a function `g`, it returns **another function** that when evaluated at **x** returns the Lagrange interpolation polynomial of `g` defined as $$ \mathcal{L} g(x) := \sum_{i=0}^n g(x_i) l_i(x)$$You could use this function as follows::    Lg = lagrange_interpolation(xi, g)    xi = linspace(0,1,101)    plot(x, g(x))    plot(x, Lg(x))    plot(xi, g(xi), 'or')<jupyter_code>def lagrange_interpolation(xi,g):
    def L_g(x,xi,g):
        summ = zeros(x.shape)
        for i in range(len(xi)):
            summ += g(xi[i])*lagrange_basis(xi, i, x)
        return summ
    return lambda x: L_g(x,xi,g)
    

# Check for polynomials. This should be **exact**
g = lambda x: x**3+x**2

xi = linspace(0,1,4)
Lg = lagrange_interpolation(xi, g)

x = linspace(0,1,1001)

assert p_norm(g(x) - Lg(x),float('inf')) < 1e-15, 'This should be zero...'<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13777003484320557, 'fraction_numerical': 0.09526132404181184, 'mean_word_length': 3.730059327620303, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 77, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '492342', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 6216, 'n_tokens_neox': 5350, 'n_words': 1741, 'path': 'assignments/release/Basics/D1-exercises.ipynb', 'repo_name': 'denocris/Introduction-to-Numerical-Analysis', 'size': '14350'}