12. Integrali#
""" Import libraries """
import numpy as np
12.1. Integrazione di Newton-Cotes#
Valutazioni di integrali definiti di funzioni \(f: [a,b] \in \mathbb{R} \rightarrow \mathbb{R}\), continue sull’intervallo. Si veda l”esempio per intuire la necessità della continuità della funzione.
Formula del punto medio
Formula del trapezio
12.1.1. Metodo di integrazione del punto medio#
Il metodo di integrazione del punto medio ricorda la definizione dell’integrale di Riemann (todo aggiungere link). Data la funzione \(f:[a,b] \in \mathbb{R} \rightarrow \mathbb{R}\), e una partizione \(P = \{ a = x_0 < x_1 < \dots < x_n = b \}\) dell’intervallo \([a,b]\), l’integrale viene calcolato come la somma dell’area dei rettangoli elementari di lati \(\Delta x_k := x_k - x_{k-1}\) e \(f(\xi_k)\), con \(\xi_k \in [x_{k-1}, x_k]\),
""" Mid-point method """
def integral_rect(f, a, b, n):
"""
Inputs:
- f: function
- a, b: extreme points of the range
- n: n. of intervals of the range
"""
# Partition of the range [a,b]
# uniform partition here; refined algorithms may use non-uniform partitions
xv = a + ( b - a ) * np.arange(n+1) / n; xv[-1] = b
xc = .5 * ( xv[1:] + xv[:-1] )
dx = xv[1:] - xv[:-1]
return np.sum( dx * f(xc) )
12.1.2. Metodo di integrazione del trapezio#
Data la funzione \(f:[a,b] \in \mathbb{R} \rightarrow \mathbb{R}\), e una partizione \(P = \{ a = x_0 < x_1 < \dots < x_n = b \}\) dell’intervallo \([a,b]\), l’integrale viene calcolato come la somma dell’area dei trapezi rettangoli elementari con basi \(f(x_{k-1})\), \(f(x_k)\) e altezza \(\Delta x_k := x_k - x_{k-1}\),
""" Trapezoid method """
def integral_trapz(f, a, b, n):
"""
Inputs:
- f: function
- a, b: extreme points of the range
- n: n. of intervals of the range
"""
# Partition of the range [a,b]
# uniform partition here; refined algorithms may use non-uniform partitions
xv = a + ( b - a ) * np.arange(n+1) / n; xv[-1] = b
dx = xv[1:] - xv[:-1]
return np.sum( .5 * dx * ( f(xv[1:]) + f(xv[:-1]) ) )
12.1.3. Esempi#
12.1.3.1. Esempio 1.#
Si valuta l’integrale
con i metodi del punto medio e del trapezio. La funzione \(f(x) = x^2\) è continua ovunque e quindi è continua nell’intervallo \([0,1]\).
f = lambda x: x**2
a, b = 0., 1.
n = 10
# Evaluate integral with the rect and trapz methods
val_rect = integral_rect(f, a, b, n)
val_trapz = integral_trapz(f, a, b, n)
print(f"Value of the integral, x \in [{a}, {b}] with {n} intervals")
print(f"- mid-point method: {val_rect}")
print(f"- trapezoid method: {val_trapz}")
Value of the integral, x \in [0.0, 1.0] with 10 intervals
- mid-point method: 0.3325
- trapezoid method: 0.33499999999999996
12.1.3.2. Esempio 2. Necessità della continuità della funzione.#
12.2. Integrazione di Gauss#
L’integrazione di Gauss permette di calcolare in maniera esatta l’integrale di una funzione polinomiale \(p_n(x)\) su un intervallo \([a,b]\), come somma pesata della funzione valuatata in alcuni punti dell’intervallo,
per un numero di nodi di Gauss, \(n_{Gauss}\) che soddisfa la relazione
Per motivi di generalizzazione dell’algoritmo, nella definizione dei pesi \(w_i\) e dei nodi di Gauss \(x_i\), l’integrale viene riportato all’integrale su un intervallo di riferimento, tramite una trasformazione di coordinate. Per domini 1D, l’intervallo di riferimento per la quadratura di Gauss è l’intervallo \(\xi = [-1, 1]\) e il cambio di variabili è
così che l’integrale originale può essere scritto come
""" Gauss integration """
# Dict of Gauss weights and nodes on the reference interval [-1,1] for 1D integration
gauss_nw = { 1: {'nodes' : [ 0. ],
'weights': [ 2. ]},
2: {'nodes' : [ -1./np.sqrt(3.), 1./np.sqrt(3.) ],
'weights': [ 1., 1.]},
3: {'nodes' : [ 0., -np.sqrt(3./5.), np.sqrt(3./5.) ],
'weights': [ 8./9., 5./9., 5./9. ]},
}
def gauss_int_1(f, n, a=-1, b=1, gauss_nw=gauss_nw):
"""
Integral of the function f(x)
over a physical domain, x \in [a, b],
using n number of Gauss nodes
"""
# Cap n.nodes to the max n.nodes stored in gauss_nw dict
n = np.min([n, np.max(list(gauss_nw.keys()))])
# Gauss nodes (from reference to actual domain) and weights
xg = .5 * ( a + b + ( b - a ) * np.array(gauss_nw[n]['nodes']) )
wg = np.array(gauss_nw[n]['weights'])
return .5 * (b-a) * np.sum( wg * f(xg) )
def gauss_int_n(f, n_gauss, a, b, n_elems):
"""
Integral of the function f(x)
over a physical domain, x \in [a, b], (uniformly) splitted in n_elems
using n number of Gauss nodes per each elem
"""
# Partition of the interval [a,b]
xv = a + ( b - a ) * np.arange(n_elems+1) / n_elems; xv[-1] = b
return np.sum([ gauss_int_1(f, n_gauss, xv[i], xv[i+1]) for i in np.arange(n_elems) ])
12.2.1. Esempio 1 - integrazione di Gauss.#
Si valuta l’integrale
con il metodo di integrazione di Gauss. Si chiede di:
osservare che l’integrazione è esatta, a meno degli arrotondamenti dovuti all’aritmetica finita, poiché la funzione integranda è di grado 2, e il numero di nodi di Gauss è \(n_{Gauss} = 3\); questo è vero indipendentemente dal numero di sotto-intervalli;
cambiare le funzioni e il numero di nodi usati nelle funzioni per l’integrazione di Gauss per verificare che l’integrazione è esatta per funzioni polinomiali di grado \(2 n_{Gauss} - 1\).
""" Evaluate integrals using Gauss integration """
f = lambda x: x**2
a, b = 0., 1.
# Gauss nodes, and n.of elements
n_gauss = 3
n_elems = 10
val_1 = gauss_int_1(f, n_gauss, 0., 1.)
val_n = gauss_int_n(f, n_gauss, 0., 1., n_elems)
print(f"Value of the integral, x \in [{a}, {b}] using Gauss integration methods with {n_gauss} nodes per elem")
print(f"- using one elem: {val_1}")
print(f"- using {n_elems} elems: {val_n}")
Value of the integral, x \in [0.0, 1.0] using Gauss integration methods with 3 nodes per elem
- using one elem: 0.33333333333333337
- using 10 elems: 0.33333333333333337