8. Sistemi lineari#
La soluzione di sistemi lineari è un problema che compare in molte altre applicazioni di calcolo numerico.
Formalismo matriciale. Con il formalismo matriciale, un sistema di equazioni lineari può essere scritto come
Classificazione. In generale, i sistemi di equazioni lineari possono essere classificati:
in base al numero di incognite \(n_u\) ed equazioni indipendenti \(n_e\): \(n_e = n_u\) sistemi determinati con un’unica soluzione; \(n_e > n_u\) sistemi sovradeterminati: con nessuna soluzione in generale; \(n_e < n_u\) sistemi indeterminati, con infinite soluzioni in generale
in base alla «struttura» del sistema:
diagonale, tridiagonale, …
in base al numero di coefficienti non-nulli della matrice \(\mathbf{A}\): sistemi con matrice \(\mathbf{A}\) piena o sparsa; questa distinzione non è netta, ma il più delle volte risulta chiara dalla particolare applicazione/metodo.
Algoritmi. Esistono due grandi classi di metodi/algoritmi per la soluzione di sistemi lineari:
i metodi diretti, che si basano su una fattorizzazione della matrice
i metodi indiretti, che si basano sul calcolo di prodotti matrice-vettore
8.1. Sistemi lineari quadrati con matrici piene#
In questa sezione si discute la soluzione di sistemi lineari quadrati con matrici piene con le funzoni disponibili nella libreria NumPy.
8.1.1. Esempio 1. Sistema quadrato determinato#
Il sistema lineare
può essere riscritto con il formalismo matriciale nella forma \(\mathbf{A} \mathbf{x} = \mathbf{b}\),
e risolto grazie alla funzione \(\texttt{solve(A, b)}\) della libreria \(\texttt{numpy.linalg}\).
"""
Linear systems with full square non-singular matrices
"""
import numpy as np
A = np.array([[1.,2.,0.], [1.,0.,1.], [1.,1.,1.]])
b = np.array([0.,-1.,1.])
x = np.linalg.solve(A, b)
print(f"Sol, x: {x}")
print(f"Proof : Ax = {A @ x}") # Check that Ax = b
print(f" b = {b}")
Sol, x: [-4. 2. 3.]
Proof : Ax = [ 0. -1. 1.]
b = [ 0. -1. 1.]
8.1.2. Esempio 2. Sistemi quadrati non determinati#
I sistemi lineari
sono due sistemi quadrati non determinati. Il primo sistema non ha soluzioni, mentre il secondo ne ha infinite della forma
Dopo aver riscritto i sistemi lineari con il formalismo matriciale, si può provare a risolverli usando la funzione \(\texttt{solve(A, b)}\) della libreria \(\texttt{numpy.linalg}\). In entrambi i casi, la funzione \(\texttt{solve(A,b)}\) resitiuisce un errore, segnalando che la matrice del sistema lineare è singolare, definizione equivalente di sistemi non determinati.
todo dare interpretazione geometrica, fare grafico?
todo spiegare motivo?
Esistono algoritmi che trovano almeno una soluzione nel caso in cui ne esistano infinite?: discutere gli algoritmi implementati nella funzione \(numpy.linalg.solve()\) e rimandare alla documentazione della libreria; discutere altri algoritmi che rendono possibile trovare una soluzione
Esistono algoritmi che trovano una soluzione approssimata nel caso in cui non ne esistano?: minimi quadrati, minimizzano l’errore, dare un’interpretazione geometrica
"""
Linear systems with full square singular matrices
"""
import numpy as np
A = np.array([[1.,2.,0.], [1.,0.,1.], [2.,2.,1.]])
b = np.array([1.,-1.,1.])
# b = np.array([0.,-1.,1.])
x = np.linalg.solve(A, b)
print(f"Sol, x: {x}")
print(f"Proof : Ax = {A @ x}") # Check that Ax = b
print(f" b = {b}")
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
Cell In[2], line 11
8 b = np.array([1.,-1.,1.])
9 # b = np.array([0.,-1.,1.])
---> 11 x = np.linalg.solve(A, b)
13 print(f"Sol, x: {x}")
14 print(f"Proof : Ax = {A @ x}") # Check that Ax = b
File <__array_function__ internals>:200, in solve(*args, **kwargs)
File ~/.local/lib/python3.8/site-packages/numpy/linalg/linalg.py:386, in solve(a, b)
384 signature = 'DD->D' if isComplexType(t) else 'dd->d'
385 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 386 r = gufunc(a, b, signature=signature, extobj=extobj)
388 return wrap(r.astype(result_t, copy=False))
File ~/.local/lib/python3.8/site-packages/numpy/linalg/linalg.py:89, in _raise_linalgerror_singular(err, flag)
88 def _raise_linalgerror_singular(err, flag):
---> 89 raise LinAlgError("Singular matrix")
LinAlgError: Singular matrix
8.2. Matrici sparse#
Una matrice sparsa ha un elevato numero di elementi nulli. Una matrice sparsa viene definita in maniera efficiente salvando in memoria solo gli elementi non nulli (limiti di memoria); gli algoritmi per le matrici sparse risultano spesso efficienti perché evitano un molte operazioni che darebbero risultati parziali nulli (velocità).
todo dire due parole sui formati
todo fare esempio di calcolo del prodotto matrice vettore per matrici sparse
8.2.1. Esempio 1 - Matrice di rigidezza di elementi finiti#
Il sistema lineare
è descritto da una matrice, \(N=5\), \(N \times N = 25\), che ha \(N+(N-1)+(N-1) = 13\) elementi non nulli. Il rapporto tra il numero di elementi non nulli e il numero di elementi totali è \(\frac{3N-2}{N^2} \sim \frac{3}{N}\). Al crescere della dimensione del problema, la matrice \(\mathbf{A}\) diventa sempre più sparsa e diventa sempre più conveniente definirla come matrice sparsa, ed usare gli algoritmi pensati per questo tipo di matrici.
"""
Linear systems with square non-singular matrices, in sparse format
"""
from scipy import sparse
# Printout level: the higher the number, the more verbose the script
printout_level = 1
n_nodes = 5
i_nodes = list(np.arange(5))
# Build sparse stiffness matrix, I: row indices, J: col indices, E: matrix elems
I = np.array(i_nodes+i_nodes[:-1]+i_nodes[ 1:])
J = np.array(i_nodes+i_nodes[ 1:]+i_nodes[:-1])
E = np.array(n_nodes*[2]+(n_nodes-1)*[-1]+(n_nodes-1)*[-1])
A = sparse.coo_array((E, (I,J))).tocsr()
if ( printout_level > 50 ): # print matrix in sparse format
print(f" I: {I}\n J: {J}\n E: {E}")
print(f" A:\n {A}")
if ( printout_level > 60 ): # convert and primt matrix in full format
print(f" A.todense(): {A.todense()}")
# RHS
b = np.array(5*[1])
# Solve linear system
x = sparse.linalg.spsolve(A, b)
print(f"Sol, x: {x}")
Sol, x: [2.5 4. 4.5 4. 2.5]