Course Webpage

Introduction to python for numerical computing I

Using some numerical examples, we will review some Python concepts. We will also introduce some specific numerical computation modules.

We import the modules matplotlib.pyplot, numpy and time

In [2]:
import numpy as np
import matplotlib.pyplot as plt
import time

Vectorization

We define the function $f(x)=e^x$ as a lambda function

In [3]:
f = lambda x: np.exp(x)

We create a vector with $300000$ equispaced elements in the interval $[-10,10]$

In [4]:
x = np.linspace(-10,10,300000)

Let us create see three different ways of creating a vector that conteins the values of the function f for the values of x

The syntax for the loop for is

In [5]:
for i in range(4):
    print(i)

0
1
2
3

We start with a numpy array of size zero. We will add elements with np.append. We measure the execution time t1 as the difference in the machine time before and after the loop

In [6]:
y = np.array([])                 

t = time.time()
for i in range(len(x)):           
    y = np.append(y,f(x[i]))
t1 = time.time() - t

Now we reserve memory creating a numpy array of zeros of the same shape and size than x and we fill the elements one by one with the loop

In [7]:
y = np.zeros_like(x)      

t = time.time()
for i in range(len(x)):           
    y[i] = f(x[i])
t2 = time.time() - t

And, finally, we use the numpy vectorization: its functions (and lambda is a numpy function) can work for numpy arrays.

In [8]:
t = time.time()
y = f(x)
t3 = time.time() - t

We compare the execution time

In [9]:
print('Vector growing inside the loop: ', t1,' seconds')
print('Without vectorization:          ', t2,' seconds') 
print('With vectorization:             ', t3,' seconds') 
print('t1/t2 = ', int(t1/t2))
print('t2/t3 = ', int(t2/t3))
print('t1/t3 = ', int(t1/t3))

Vector growing inside the loop:  38.43652129173279  seconds
Without vectorization:           0.5062990188598633  seconds
With vectorization:              0.002621173858642578  seconds
t1/t2 =  75
t2/t3 =  193
t1/t3 =  14663

We can see that we have executed the options from worse to better. Thus, we will avoid growing vectors inside a loop and use vectorization.

One-dimensional function graph

plot function, from module matplotlib.pyplot, really plots a polygonal line that joins function points. For instance:

We define the function to plot, $f(x)=e^x$

In [13]:
f = lambda x: np.exp(x)

We choose the interval $[a,b]$ where we will plot the function

In [14]:
a = -1.; b = 1.

We create a five equally-spaced point mesh in the $[a,b]$ interval, that is stored in a numpy array x.

In [15]:
x = np.linspace(a,b,5)
print(x)

[-1.  -0.5  0.   0.5  1. ]

The corresponding numpy array for y is computed with

In [16]:
y = f(x)
print(y)

[0.36787944 0.60653066 1.         1.64872127 2.71828183]

Let us plot the function

In [17]:
plt.figure()
plt.plot(x,y)
plt.show()

The approach is quite rough with only 5 points. By default, linspace creates a 50 point mesh. We repeat the process with 50 points.

In [18]:
x = np.linspace(a,b)
y = f(x)

plt.figure()
plt.plot(x,y)
plt.show()

It is also a polygonal line, but now there are enough points to look like a smooth curve.

Let us add some details to the graph

In [19]:
OX = 0*x                               # 50 zeros numpy array. To plot OX axis

plt.figure()
plt.plot(x,y, label = 'f')
plt.plot(x,OX,'k', label = 'OX axis')  # k, from "black" 
plt.title('Example function f plot') 
plt.legend()                           # using the "label", it plots a legend
plt.show()

Taylor series

McLaurin series (Taylor series with $x_0=0$) for function $f(x)=e^x$ is

$$ e^x = 1 + \dfrac{x^1}{1!}+\dfrac{x^2}{2!}+\dfrac{x^3}{3!}+\dfrac{x^4}{4!}+\dfrac{x^5}{5!}+\cdots $$

Theoretically, if we add enough terms we can approximate the function as much as we want. The successive polynomials are

$$ \begin{array}{rcl} P_{0}(x) & = & 1\\ P_{1}(x) & = & 1 + \dfrac{x^1}{1!}\\ P_{2}(x) & = & 1 + \dfrac{x^1}{1!}+\dfrac{x^2}{2!}\\ P_{3}(x) & = & 1 + \dfrac{x^1}{1!}+\dfrac{x^2}{2!}+\dfrac{x^3}{3!}\\ P_{4}(x) & = & 1 + \dfrac{x^1}{1!}+\dfrac{x^2}{2!}+\dfrac{x^3}{3!}+\dfrac{x^4}{4!}\\ & \ldots & \\ \end{array} $$

and $P_n$ approaches the function better than polynomial $P_{n-1}$. Let us compute the value of degree 3 polynomial in $x_0=0.5$. Taking into account that $x_0^0 = 1$ and $0!=1$

In [20]:
x0 = 0.5
polynomial = 0.
factorial  = 1.

for i in range(4):
    term = x0**i/factorial
    polynomial += term
    factorial *= i+1
    
    
print('P3(0.5)     = ', polynomial)
print('np.exp(0.5) = ', np.exp(x0))

P3(0.5)     =  1.6458333333333333
np.exp(0.5) =  1.6487212707001282

Let's transform the script into a function and let us try with a higher degree polynomial

In [21]:
def P(x0,degree):
    polynomial = 0.
    factorial = 1.

    for i in range(degree + 1):
        term = x0**i/factorial
        polynomial += term
        factorial *= i+1
        
    return polynomial
In [22]:
print('P(0.5, 10)  = ', P(0.5, 10))
print('np.exp(0.5) = ', np.exp(0.5))

P(0.5, 10)  =  1.6487212706873655
np.exp(0.5) =  1.6487212707001282

If the input argument is a numpy array, the output is also a numpy array

In [23]:
a = -1.; b = 1.
x = np.linspace(a,b,5)
print('x         = ', x)
print('P(x, 10)  = ', P(x, 10))
print('np.exp(x) = ', np.exp(x))

x         =  [-1.  -0.5  0.   0.5  1. ]
P(x, 10)  =  [0.36787946 0.60653066 1.         1.64872127 2.7182818 ]
np.exp(x) =  [0.36787944 0.60653066 1.         1.64872127 2.71828183]

If we use a 50 points numpy array, we can draw the polynomial. Let us graph the 2 degree polynomial using a red line

In [24]:
a = -1.; b = 1.
f = lambda x: np.exp(x)
x = np.linspace(a,b)
OX = 0*x    

plt.figure()
plt.plot(x,f(x), label = 'f')
plt.plot(x,OX,'k')
plt.plot(x,P(x,2),'r', label = 'P2')
plt.title('Example of function and polynomial graph')
plt.legend()                           
plt.show()

No we change the interval to $[-3,3]$ and, using a loop, we plot more polynomials

In [25]:
a = -3.; b = 3.
f = lambda x: np.exp(x)
x = np.linspace(a,b)
OX = 0*x    

plt.figure()
plt.plot(x,f(x), label = 'f')
plt.plot(x,OX,'k') 

for degree in range(1,5):
    plt.plot(x,P(x,degree), label = 'P'+str(degree))
    
plt.title('Function and approximation polynomials') 
plt.legend()                           
plt.show()

If we add plt.pause(1) we can see a small animation on screen using spyder. For this purpose:

  • Open spyder.
  • In menu, Tools $\rightarrow$ Preferences $\rightarrow$ IPython console $\rightarrow$ Graphics $\rightarrow$ Graphics backend $\rightarrow$ Backend $\rightarrow$ Automatic.
  • Close spyder and open it again.
  • Copy, paste and run the following script.
import numpy as np
import matplotlib.pyplot as plt

#%%--------------------------------------
def P(x0,degree):
    polynomial = 0.
    factorial = 1.

    for i in range(degree + 1):
        term = x0**i/factorial
        polynomial += term
        factorial *= i+1

    return polynomial 

#%%--------------------------------------
a = -3.; b = 3.
f = lambda x: np.exp(x)
x = np.linspace(a,b)
y = f(x)

OX = 0*x                               

plt.figure()
plt.plot(x,y, label = 'f')
plt.plot(x,OX,'k') 

for degree in range(1,7):
    plt.plot(x,P(x,degree), label = 'P'+str(degree))
    plt.title('Function and approximation polynomials') 
    plt.legend()   
    plt.pause(1)                        
plt.show()

An example of while loop

In [26]:
i = 0
while i < 5:
    print(i)
    i += 1

0
1
2
3
4

Exercise 1

Using a while loop, write a script, that using the McLaurin polynomials for the function $f(x)=e^x$, computes its approximate value at the point x0 = -0.4. The stopping criteria are that the last added term in absolute value is less than a tolerance value tol=1.e-8 and the maximum number of terms is 100, that is maxNumSum=100. Compare the result with the function lambda f value at the same point. Print the number of iterations.

Hints:

  • Use some logical operator and, or, not, for the while condition.
  • The absolute value of a real number a can be obtained with np.abs(a)
  • McLaurin series for function $f(x)=e^x$ was
$$ e^x = 1 + \dfrac{x^1}{1!}+\dfrac{x^2}{2!}+\dfrac{x^3}{3!}+\dfrac{x^4}{4!}+\dfrac{x^5}{5!}+\cdots $$
In [27]:
%run Exercise1.py

Function value in -0.4      = 0.6703200460356393
Approximation value in -0.4 = 0.67032004600776
Number of iterations        = 10

Exercise 2

Modify the script of Exercise 1 to compute simultaneously the function approximate value in the 50 points contained in x = np.linspace(-1,1). The stopping criteria are that the last term in absolute value (now we have 50 different $x_0$ and 50 different last terms), is less than a tolerance tol=1.e-8 and that the maximum number of terms is 100, maxNumSum=100.

Transform the script into a function funExp(x, tol, maxNumSum) whose input arguments are the numpy array x of 50 $x$ values, tol and maxNumSum and whose output argument is a numpy array y of the approximate values of the function in the $x$ points obtained using McLaurin polynomial. Use this function to plot the function. Plot also the corresponding lambda function defined from np.exp(x).

Hints:

  • np.max(np.abs(term)) gives the maximum value of all the absolute values of the values contained in the numpy array term.
  • To plot a thick yellow line plt.plot(x,y,'y',linewidth=5)
  • To plot a dashed blue line plt.plot(x,y,'b--')
In [28]:
%run Exercise2.py

Proposed exercises

Exercise 3

Using a while loop, write a script, that using the McLaurin polynomials for the function $f(x)=\sin(x)$, computes its approximate value at the point x0 = np.pi/4. The stopping criteria are that the last added term in absolute value is less than a tolerance value tol=1.e-8 and the maximum number of terms is 100, that is maxNumSum=100. Compare the result with the function lambda f value at the same point. Print the number of iterations.

Hints:

  • Use some logical operator and, or, not, for the while condition.
  • McLaurin series for function $f(x)=\sin(x)$ was
$$ \sin(x) = \dfrac{x^1}{1!}-\dfrac{x^3}{3!}+\dfrac{x^5}{5!}-\dfrac{x^7}{7!}\cdots $$
In [29]:
%run Exercise3.py

Function value       =  0.7071067811865475
Approximation value  =  0.7071067811796194
Number of iterations =  6

Exercise 4

Modify the script of the previous exercise to compute simultaneously the function approximate value in the 50 points contained in x = np.linspace(-np.pi,np.pi). The stopping criteria are that the last term in absolute value (now we have 50 different $x_0$ and 50 different last terms), is less than a tolerance tol=1.e-8 and that the maximum number of terms is 100, maxNumSum=100.

Transform the script into a function funSin(x, tol, maxNumSum) whose input arguments are the numpy array x of 50 $x$ values, tol and maxNumSum and whose output argument is a numpy array y of the approximate values of the function in the $x$ points obtained using McLaurin polynomial. Use this function to plot the function. Plot also the corresponding lambda function defined from np.sin(x).

In [30]:
%run Exercise4.py

Exercise 5

McLaurin series for function $\mathrm{sinh}$ is

$$\mathrm{sinh}\;x = \dfrac{x^1}{1!}+\dfrac{x^3}{3!}+\dfrac{x^5}{5!}+\dfrac{x^7}{7!}+\dfrac{x^9}{9!}+\cdots$$

and for $\mathrm{cosh}$

$$\mathrm{cosh}\;x = 1+\dfrac{x^2}{2!}+\dfrac{x^4}{4!}+\dfrac{x^6}{6!}+\dfrac{x^8}{8!}+\cdots$$

and as $\mathrm{tanh}$

$$ \mathrm{tanh}=\dfrac{\mathrm{sinh}\;x}{\mathrm{cosh}\;x} $$

compute the $\mathrm{tanh}$ for $x_0=0.5$ from McLaurin series of sinh and cosh, using the same number of terms in the numerator and denominator. Stop when two consecutive approximations of $\mathrm{tanh}(x_0)$ differ by less than $10^{-4}$.

In [31]:
%run Exercise5.py

Approximate value  =  0.4621171922898218
Exact value.       =  0.46211715726000974

Exercise 6

Modify the script of the previous exercise to compute simultaneously the function approximate value in the 50 points contained in x = np.linspace(-3,3). The stopping criterion is that the difference between two consecutive approximations is less than a tolerance tol=1.e-8 for all the 50 points in x.

Transform the script into a function funTanh(x, tol) whose input arguments are the numpy array x of 50 $x$ values, tol and whose output argument is a numpy array y of the approximate values of the function in the $x$ points obtained using McLaurin polynomial. Use this function to plot the function. Plot also the corresponding lambda function defined from np.tanh(x).

In [32]:
%run Exercise6.py