Commit 655cbf29 authored by paugier's avatar paugier
Browse files

Add subjects pratical sessions

parent 3f639f57
Pipeline #77670 passed with stage
in 1 minute and 15 seconds
%% Cell type:markdown id: tags:
# Scientific Computing : integration methods
%% Cell type:markdown id: tags:
In this practical session we are going to study several methods to approach numerically the computation of:
\begin{equation}
I = \int_a^b f(x)
\end{equation}
This situation can happen, for exemple, in the following contexts:
- One only knows the function $f$ in certain specific points. This is the case when $f$ results from a numerical simulation.
- No explicit form of the primitive is known : e.g. $f(x)=e^{-x^2}$, $f(x)=\frac{e^x}{x}$
The idea here is to get an approximation $I_{app}$ of the exact value of the integral $I$ with a precision $\epsilon$:
\begin{equation}
| I - I_{app} | \le \epsilon
\end{equation}
To do so, one can replace $f$ by a function from which it is easy to compute the integral. This is can be done through a polynomial interpolation.
# Left Rectangle Method
The following scheme is one method to approximate the integral of a function $f$ between two finite bounds $a$ and $b$:
\begin{equation}
\int_a^b f(x) dx \approx (b-a) \times f(a)
\label{lrm}
\end{equation}
It is a quadrature method of order 0, meaning that one approximates the function by a piecewise constant function (0-order polynomial function). It is called the left rectangle method (LRM). One of your task will be to illustrate why it called as such.
However, polynomial interpolation gives satisfying results only on short intervals. This implies that our domain of integration must be discretized in order to achieve decent approximations (i.e. splitted into small subdomains of integration):
\begin{equation}
\int_a^b f(x)dx = \sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}}f(x)dx
\end{equation}
where
$$a = x_0 < x_1 < \cdots < x_n = b$$
It is common to denote $h$ the distance between two successive discretization points:
$$h_i=|x_{i+1} - x_i| << 1$$
**Question** : What becomes Eq. \ref{lrm} when applied on a such a discretized space ?
\begin{equation}
\int_a^b f(x)dx = \sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}}f(x)dx \approx \sum_{i=0}^{n-1} ?
\end{equation}
%% Cell type:markdown id: tags:
**Answser**:
\begin{equation}
\int_a^b f(x)dx = \sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}}f(x)dx \approx \sum_{i=0}^{n-1} (x_{i+1} - x_{i})f(x_i)
\end{equation}
%% Cell type:markdown id: tags:
Let us define the two functions that you will use in this practical session:
\begin{equation}
f_1(x)=exp(x)
\label{f1}
\end{equation}
\begin{equation}
f_2(x) = -x^3 + 3 x^2 -6x + 7
\label{f2}
\end{equation}
You will consider the integration of $f_1$ (Eq. \ref{f1}) over $[0,2]$ and $f_2$ (Eq. \ref{f2}) over $[-5,5]$.
%% Cell type:markdown id: tags:
**Q)** Create python function (named $f_x$) for $f_1$ and $f_2$. Plot them on the interval considered in this practice session. Compute and create functions for their respective primitives ( _named prim_fx_ ). Compute the exact value of the integrals ( _named I_fx_ ) on the intervals specified before.
<div class="alert alert-block alert-info">
<b>Tip:</b> For the plots, have a look at the <a href="https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.pyplot.plot.html">Matplotlib Website</a>
</div>
<div class="alert alert-block alert-info">
<b>Recall:</b> The primitive $F$ of a function $f$ is defined as $F'(x)=f(x)$.
</div>
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
import numpy as np
# Define parameters
a1 = 0; a2 = -5
b1 = 2; b2 = 5
# For function f1
def f1(x):
return None
def prim_f1(x):
return None
# For function f2
def f2(x):
return None
# Make a function for the exact solution of f2
def prim_f2(x):
return None
# Define exact values of the integrals
I_f1 = None
I_f2 = None
```
%% Cell type:markdown id: tags:
**Q)** Complete the following function (**using a for loop**):
%% Cell type:code id: tags:
``` python
def left_rect_meth(x, func):
""" Function to approximate an integral using the left rectangle method.
Computations will be realized using a for-loop.
Parameters
----------
x : numpy array
Array representing the discretization of the integration domain
[a, a+h, a+2h, ..., b-h, b] where a and b are the bounds of
integration.
func : python function
Function to be integrated
Output
----------
I : float
Approached value of the integral
"""
return None
```
%% Cell type:markdown id: tags:
**Q**) Choose a discretization for the domain (with, for example, 99 subintervals) and apply the left rectangle method to compute an approximation of $$\int_0^2 f_1(x) dx$$
Compare it to the exact value.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
**Q)** In order to practice plotting in Python, create a function that will plot the behaviour of the left rectangle integration method. This plot will contain the following elements:
- A fine representation of the considered function $f$
- A representation of the function on the discretization used for integration
- Rectangles (or patches) that corresponds to the ones used for the integration
Use this function to schematize an integration of $f_1$ and $f_2$.
%% Cell type:markdown id: tags:
<div class="alert alert-block alert-info">
<b>Tip:</b> look for matplotlib patches on the documentation.
</div>
%% Cell type:code id: tags:
``` python
import matplotlib.patches as patches
```
%% Cell type:markdown id: tags:
# Simpson's formula
Another numerical quadrature formula is the Simpson's rule:
\begin{equation}
\int_a^b f(x) dx \approx (b-a) \left( \frac{f(a)}{6} + \frac{4}{6}f\left(\frac{a+b}{2}\right) +\frac{f(b)}{6} \right)
\end{equation}
**Q)** Complete the following python function so that it computes the approximation of an integral via the simpson method **using a for loop**.
%% Cell type:code id: tags:
``` python
def Simpson(x, func):
""" Function to approximate an integral using the Simpson's rule.
\int_a^b f(x) dx ~= (b-a)*(f(a)/6 + 4/6*f((a+b)/2) +f(b)/6)
Parameters
----------
x : numpy array
Array representing the discretization of the integration domain
func : python function
Function to be integrated
Output
----------
I : float
Approached integral computed
"""
I = 0
return None
```
%% Cell type:markdown id: tags:
# Built in methods in python
In practice, one does not usually code numerical methods by him or herself. Indeed several packages provide a lot of different numerical methods already implemented and optimized. Among this packages, scipy is one of the most common in Python. In general, methods provided by these packages will perform better that user-defined ones. However, it is always a good practice to have some knowledge about what kind of methods are being used and what are the potential problems associated with it.
**Q)** Based on the online documentation, compute the integral of either $f_1$ or $f_2$ over the domains defined previously using scipy built-in functions.
%% Cell type:code id: tags:
``` python
import scipy.integrate as integrate
```
%% Cell type:markdown id: tags:
# Method efficiency
Another very important aspect in scientific computing is the computing time required for a specific computation. Sometimes it is very convenient to take advantage of mathematical properties or programming language properties to improve efficiency. For example, in Python, loops are much less efficient than array programming. Let's illustrate that point.
%% Cell type:code id: tags:
``` python
import time
# Create some vectors populated with random values
n_elements = int(1e6) # values per vector
vec_a = np.random.randint(low=0, high=10, size=n_elements)
vec_b = np.random.randint(low=0, high=10, size=n_elements)
# Let say that we want to obtain their dot product.
# First we will compute it with a loop
dot = np.empty(vec_a.shape)
start = time.time()
for k in range(len(dot)):
dot[k] = vec_a[k] * vec_b[k]
end = time.time()
elapsed_time_loop = end - start
print(f"Elapsed time with a loop : {elapsed_time_loop}")
# And now with array programming.
start = time.time()
dot_array = vec_a * vec_b
end = time.time()
elapsed_time_array = end - start
print(f"Elapsed time with array programming : {elapsed_time_array}")
```
%% Cell type:markdown id: tags:
Using array programming is actually much more efficient ! Now, you want to make use of that knowledge to reduce the computational time of your integration routines.
**Q)** Write a function to compute the left rectangle method **without using a for loop**:
%% Cell type:code id: tags:
``` python
def left_rect_meth_vec(x, func):
""" Function to approximate an integral using the left rectangle method.
Parameters
----------
x : numpy array
Array representing the discretization of the integration domain
[a, a+h, a+2h, ..., b-h, b] where a and b are the bounds of
integration.
func : python function
Function to be integrated
Output
----------
I : float
Approached integral computed
"""
return None
```
%% Cell type:markdown id: tags:
**Q)** For a discretization with about 5e5 points compare the time required to computed the approximated value of integral of $f_1$ between $0$ and $2$. _Hint: use the function time_
%% Cell type:code id: tags:
``` python
import time
```
%% Cell type:markdown id: tags:
**Q)** We now assume that the discretisation is done on a regular grid (i.e. each point of the discretization are evenly space). Show that under this assumption the Simpson's formula can also be computed using array programming. Complete the following function **without using loops**.
%% Cell type:code id: tags:
``` python
def Simpson_reg(N, a, b, func):
""" Function to approximate an integral using the Simpson's rule on a
uniform grid.
Parameters
----------
N : int
Number of discretization points between integration boundaries.
a : float
Lower boundary for integration
b : float
Upper boundary for integration
func : python function
Function to be integrated
Output
----------
I : float
Approached integral computed
"""
return None
```
%% Cell type:markdown id: tags:
# Study of the convergence
%% Cell type:markdown id: tags:
## Convergence of the Left Rectangle Method
One wishes that the absolute error between the exact and approached solution $ \epsilon_f (h)$ tends to zero as quickly as possible when $h \rightarrow 0$. More specifically, it is common to define the convergence order $\alpha$ such as
\begin{equation}
\epsilon_f(h) = C_f \times h ^ \alpha
\end{equation}
Your objective for the next couple of questions will be to determine this convergence order.
**Q)** Compute the error of integration $\epsilon_f$ as a function of the discretization step $h$. Once computed, plot these errors on two different figures in a logarithmic scale ( _use axes.loglog() function_ ).
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
5) According to what you observe, fit a polynomial regression model ( _using numpy.polyfit_ ) of the appropriate order (justify) to the errors as a function of $h$. Add the constructed model to the previous plot to verify your results. In the legend, insert the equation of the model.
%% Cell type:markdown id: tags:
<div class="alert alert-block alert-info">
<b>Tip:</b> When fitting your polynomial model to the data, remember that your observations were drawn on a <i>loglog<i> plot.
</div>
%% Cell type:markdown id: tags:
## Convergence of the Simpson's formula
**Q)** Repeat the previous steps for the Simpson's formula
%% Cell type:code id: tags:
``` python
```
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment