Metodi numerici per l'ingegneria

Risoluzione di sistemi lineari

A_nx=b con A∈ℝ^n*m

Dati un vettore b∈ℝⁿ e una matrice A∈ℝ^m×n, si cerca X∈ℝⁿ t.c. Ax=b.

CNS per avere un'unica soluzione → A non singolare, cioè det A ≠ 0 quindi A invertibile.

Soluzione: X=A^-1b COSTOSO

Non si calcola se non è necessario, A^-1 ma si cerca di risolvere il sistema. Questo è molto costoso. Nel caso di alcune matrici si hanno però dei costi contenuti -> Di seguito sono mostrate le soluzioni di tali matrici:

• Matrice A diagonale
  • A =[ d₁₁ 0 ... ][ 0 d₂₂ ... ][ ... d_mm ]
  CNS ⇒ d_ii ≠ 0 ∀i

Sistema costituito da m equazioni indipendenti tra di loro:

• x_a = b₁/d₁₁
• x_m = b_m/d_mm
• → x_k = b_k/d_kk

Costo: m ops

• A Triangolare superiore
  • A =[ u₁₁ u₁₂ u₁₃ ... u_1m ][ 0 u₂₂ u₂₃ ... u_2m ][ 0 0 u₃₃ ... ][ 0 0 0 u_mm ]
  CNS ⇒ u_ii ≠ 0 ∀i

Sistema

• x₁ = (b₁ − ∑ ^m_j=2 l_ijx_j)/ u_1m
• ...
• x_m = b_m/u_mm

• Inversione
  • A =[ ℓ₁₁ 0 0 ... ][ ℓ₂₁ ℓ₂₂ 0 ... ][ ℓ₃₁ ℓ₃₂ ℓ₃₃ ... ][ ℓ_m1 ... ℓ_mm ]
  CNS ⇒ ℓ_ii ≠ 0 ∀i

Sistema

• x₂ = b_n/l_n
• x₂ = (b₂ − l₂₁x₁)/l₂₂
• x_m = (b_m − ∑ ^m-1_j=1 l_ijx_j) / l_mm

Si parte dall'ultima equazione

Si parte dalla prima equazione.

ξ_k = b^T

ξ_k = _T

Il costo per il calcolo del generico termine k è dato da

Nel caso di una generica matrice A si potrebbe

metodo di Cramer

usare il metodo di Cramer che permette di trovare la

data da

dove Ai è ottenuta sostituendo alla colonna i il

vettore dei termini noti. Il costo per il calcolo

di un determinante è O(m!). In questo caso si

devono calcolare m+1 determinanti (A, A1, A2, ..., Am)

costo = (m+1)

Metodo molto costoso, usato solo per sistemi al massimo

4x4.

MEG

riduco la matrice A ad una matrice triangolare

Costo = O(^m3⁄₃) + M(m+1)

I metodi visti sopra servono per risolvere sistemi del tipo Ax = b evitando di calcolare esplicitamente A^-1. Se si volesse calcolare esplicitamente i coefficienti di A^-1, si potrebbe utilizzare tale metodo.

A^-1 = [c⁽¹⁾ c⁽²⁾ ... c^(m)]

per definizione AA^-1 = I

Sostituisco ad A^-1 il vettore delle colonne e ottengo il seguente sistema:

• Ac⁽ⁿ⁾ = e⁽ⁿ⁾
• Ac⁽²⁾ = e⁽²⁾
• Ac^(m) = e^(m)

con e⁽ⁿ⁾ = [1 0 ... 0]^Tcon e⁽²⁾ = [0 1 0 ... 0]^Tcon e^(m) = [0 0 ... 1]^T

Ogni equazione mi dà una colonna della matrice. Da una equazione all'altra varia solo il vettore dei termini noti, si può quindi usare MFG e calcolare un'unica volta le 2 matrici L e U, il costo di fattorizzazione si fa quindi una volta sola.

Costo = O(^M3⁄₃) ➔ M(m+1) ➔ Mper ogni equazione

≈ ^M3⁄₃ + M³ + M² ≈ ⁴⁄₃ m³ + m²

Studio della stabilità del problema

Un problema è stabile o ben condizionato se a piccoli errori sui dati corrispondono piccoli errori sul risultato.Dato il problema Ax = b

I metodi iterativi necessitano di un test di arresto che serve per decidere quando fermare la successione e considerare abbastanza buono il risultato ottenuto. Ci sono diversi test:

• test ideale: |x^(m) - x| < ε → Problema: x non è noto;
• test sulle iterate: |x^(m) - x^(m+1)| <= εQuesto test va bene per metodi con convergenza rapida. In caso contrario si avrebbe uno scostamento tra le iterate troppo piccolo e il metodo rischierebbe di interrompersi troppo presto.
• test sul residuo: ||b - Ax^(m)|| / ||b|| = ||A*e^(m)|| / ||b|| <= ε

Generalmente si usa il secondo metodo ad ogni iterazione e se il test ha successo si verifica un controllo col test sul residuo (che è più costoso).

Calcolo manuale N O(m²): Prodotto matrice vettore N O(m²) (costo per ogni iterazione = O(m²))

Conviene utilizzare un metodo iterativo invece di uno diretto se il numero di iterazioni è << M/3.

• Sistema sottodimensionatoAx = b con b ∈ Rᵐ, A ∈ Rᵐˣⁿ con m < n. In questo caso è necessario fissare a priori m-n incognite e poi risolvere normalmente il sistema.
• Sistema sovradimensionatoAx = b con b ∈ Rᵐ, A ∈ Rᵐˣⁿ con m > n. Un sistema di questo tipo può avere uno, nessuna o infinite soluzioni a seconda di A e b. Sia A ∈ Rᵐˣⁿ, m > n, se A ha rango pieno, il sistema normale ha un'unica soluzione detta soluzione del sistema nel senso dei minimi quadrati.

Consideriamo il sistema normale Aᴵ A x = Aᴵ b

Si ottiene quindi il sistema

AA_n = E

con

• Q = [a₁, ..., a_N]^T
• E = [b₁, ..., b_n] con b_i = ∫_a^b f(x) φ_i(x) dx
• A = [a_i,j] con a_i,j = ∫_a^b φ_i(x) φ_j(x) dx

Ma siccome φ_i(x) e φ_j(x) sono ortogonali si ha:

a_i,j = { 0 se i ≠ j∫_a^b φ_i²(x) dx se j = i

Sistema:

• ∫_a^b φ₁²(x) dx
• ...
• ∫_a^b φ_N²(x) dx

Siccome A_K è diagonale, il sistema è costituito da N equazioni indipendenti che si risolvono nel seguente modo:

a_K = ∫_a^b f(x) φ_K(x) dx / ∫_a^b φ_K²(x) dx = (f, φ_K)_L²(a,b) / ‖φ_K‖²_L²(a,b)

detti coefficienti da Fourier di f

Se {φ_K} è una base ortonormale (‖φ_K‖ = 1), allora a_k = (f, φ_K) e valgono le seguenti proprietà:

→ f_N è la proiezione sullo spazio S di f cioè si ha:(f - f_N, V) = 0 ∀ V ∈ S.

Se \(d^k = d^{opt}\), l'espressione dell'errore è:

E \((x^{(m+1)}) = (1 - \gamma_k) E (x^{(m)})

Metodo del gradiente

Si prende come direzione di discesa il residuo:

\(\rho^{(m)} = y_c^{(m)} - \nabla F (x^{(m)})\)

L'espressione dell'errore è:

E \((x^{(m+1)}) = \left[ \frac{\kappa(A)-1}{\kappa(A)+1} \right]^2 E (x^{(m)})

Se \(A\) è mal condizionata \((\kappa(A) \gg 1)\) la convergenza è molto lenta.

Costo per ogni iterazione: \(O(m^2)\)

Metodo del gradiente coniugato

IDEA: Se \(d^k = d^{opt}\) allora \((x^{(m)}, \rho^{(m-1)}) = 0\). Dato \(\rho^{(m-1)}\) si cerca \(\rho^{(m)}\) nel piano

\(span\{x^{(m)}, \rho^{(m-1)}\}\).

\(\rho^{(m)} = y_c^{(m)} + \beta_k \rho^{(m-1)}\)

con \(\beta_n\) da scegliere in modo che la riduzione dell'errore sia massima. Si trova che:

\(\beta_k = - \frac{(A \rho^{(m-1)}, y_c^{(m)})}{(A \rho^{(m-1)}, \rho^{(m-1)})}\)

L'espressione dell'errore è:

E \([x^{(m)}] \leq 4 \left( \frac{\sqrt{\kappa(A)} - 1}{\sqrt{\kappa(A)} + 1} \right)^2 E (x^{(m-1)})

Il numero di iterazioni necessarie per ridurre di \(\epsilon\) l'errore iniziale è proporzionale a \(\sqrt{\kappa(A)}\), mentre

nel metodo del gradiente a \(\kappa(A)\).