Calcolo Numerico, Scuderi - Teoria con esempi

CALCOLO NUMERICO

TEORIA + MATLAB

• ARITMETICA - ERRORI
• SISTEMI LINEARI
• METODI PER RISOLUZIONE SISTEMI LINEARI
• APPROSSIMAZIONE DI DATI E DI FUNZIONI
• CRITERI DI INTERPOLAZIONE
• CONVERGENZA
• EQUAZIONI DIFFERENZIALI ORDINARIE

Lezioni Soudieri

Lezione 12 - Aritmetica Errori

Consideriamo il sitema di num. decimale, ovvero in base N=10

Es. a = 109,34 = 1x10² + 0x10¹ + 9x10⁰ + 3x10^-1 + 4x10^-2

Si def. Rappresentazione Floating-Point di un numero reale a la rappresentazione

a = pN^q p reale, q intero

Es. a = 0.015 x 10^-1 0.15 x 10^-2 = 0.0015 10⁰

La rappr. floating-point non è unica.

È univocamente determinata se si impone la condizione:

• N^-1 ≤ |p| < 1 ( => per N=10 , 0.1 ≤ |p| <1)

Nell'es. precedente l'unica che soddisfa questa condizione è 0.15x10^-2 (perchè la cifra dopo la virgola è ≠0)

Se la condizione è soddisfatta allora si dice che la rapp. floating point è NORMALIZZATA

a = pN^q → | in questo caso p=mantissa di aq= esponente o caratteristica di a

Es. Per N=100 : intero

Fissato m₀, al posto di q il calcolatore memorizzaq^* = q - m₀ ≥ 0esponente normalizzato

Non tutti i numeri reali sono rappresentabili su un calcolatore.

Definiano numeri macchina quei numeri le cui p e q sono esattamente rappresentabili, negli spazi loro riservati

Es: N=10, t = 5, m = -124, M = 128

a = 1.58291.10^-2 è num. di macchina perchè abbiamo 6 cifre e t è massima 5.

y = ^{1 - cos(x)}/_x² con x ≠ 0

Eelimino la caus. num. (dovuta a sottrazione tra 1 e cos x, evitando si nulii e autocondiz.)

1 - cos x = 2sin²(x/2)

y = ^{2 sin2(x/2)}/_x² = ¹/₂ ( ^{sin (x/2)}/_x/2 )

y = ^{x - sin x}/_{lax x} x ≠ 0

Utilizziamo Taylor → y = x ( ^x3/_3! - ^x5/_5! + ... )

= ^x3/_3! - ^x5/_5! + ...

fenomeno dei numeri num. molto min.

autoc sottrazioni

PROBLEMA NUMERICO = descrizione chiara e non ambigua di una connessione funzionale

x → y

↑ ↑

dati input dati output

∞ x ∞

f(x) f(x)

Studiare CONDIZIONAMENTO di un problema.

Se x ≈ x̃ e f(x) ≈ f(x̃) → ben condizionato

(piccole perturbazioni nei dati → piccole perturb.

nei risultati)

Se x ≈ x̃ e f(x) ≠ f(x̃) → mal condizionato

(piccole perturb. nei dati → grandi perturb. su ris.)

Per lo studio del condizionamento dei relazioni del tipo:

* ||f(x)^- β(x̃)||/f(x) ≈ K (f, x) ||x - x̃||

oppure ||f(x)^- β(x̃)||/f(x) ||^{x - x̄}/_|x|

Example: y = x₁ + x₂, x₁, x₂∈ℝ

|x₁ + x₂ - (x₁ + x₂)| = |x₁(1 + ϵ₁) - x₂(1 + ϵ₂)|

[x₁ x₂] [x₁ x₂]

x₁ = ^{3x - xⱼ}/₂

|^{x - x̄}/_|x|

Comandi MatLab per implementare questo algoritmo

A: [ . . . ; . . . ; ] b: [ . ; . ; . ] m: length(b) x: zeros(m,1) x(1) = b(1) / A(1,1) for i = 2:m     s = A(i,1) * x(1,i-1) + x(i-1)     x(i) = (b(i)-s) / A(i,i) end end

Metodo di eliminazione di Gauss

Assegniato Ax = b di ordine m, il Gauss trasforma in m-1 passi il sistema Ax = b in uno equivalente Ux = B, con U matrice Δ sup.

Ax = b ⇔ Ux = B ⇔ x    Equivalenti (cioè stessa soluzione)

Proprietà del metodo

• Soluzione invariata se scambiamo tra loro due equazioni del sistema.
• Se sostituiamo a un’equ. del sistema una c.f. dell’equ. stessa.

Consideriamo il sistema di ordine m+1 non singolare e cioè det ≠ 0.

Ax = b →

• a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + a₁₄x₄ = b₁
• a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + a₂₄x₄ = b₂
• ....

• Come lo trasformiamo in Ux = b? 1. Lascia nella I colonna solo il 1' elemento 2. Eliminiamo x₂ da III e IV equazione 3. Eliminiamo x_u da IV equazione

Passo k = 1

• Poniamo: a(¹)_ij = a_ij e b(¹)_j (notazione)
• Supponiamo c(¹)₁₁ ≠ 0, altrimenti SCAMBIO le 1° eq. con la k-esima ke c(^k)_k1 ≠ 0
• Eliminiamo x nelle eq. i = 2,3,4

In generale si considera la riga 2 come riga i-esima - quindi faccio la stessa cosa per i = 3 e 4.

Alla fine ottengo la matrice →

Alcune app della fattorizzazione PA=LU

1) Risoluzione del sistema lineare Ax=b

Ax=b ↔ PAx=Pb ↔ LUx=Pb

Ly=Pb => _n.b. m²/2 operazioni

Ux=y => _n.b. m²/2 operazioni

Costo risoluzione (O m³)

Costo fabbricazione m³ oper.

⇓ nnz (numero non zeri di L+U)

φ (m³/3) operazioni _{moltiplicazioni}

Lezione 5-6

PA=LU

P: matrice di permutazione → memorizza gli scambi richiesti dal pivoting

A: matrice di partenza

L: matrice a tc che contiene i moltiplicatori

U: matrice SυD usata x risolvere il sistema Ax=b

Altre applicazioni:

1. Calcolo del determinante di A ⇓

   Calcolo PA=LU e poi det(A)= (-1)^s∏_i=1^m u_ii

   Calcolo fatto da MATLAB quando calcoliamo det(A)

   prodotto della diagonale u₁₁ ... u_nn

   s=n. totale degli scambi effettuati

2. Calcolo dell'inversa di A (costo m³)

   A^-1=?

   Ricorda: (AB)^-1 = B^-1 A^-1 ⇐ moltiplico per P x eliminare P^-1

   PA=LU

   ↔ P^-1(PA)= P^-1(LU) ↔ A= L'U ↔ A^-1P= P^-1 P^-1 = ^-1 L^-1U^-1 P

   ↔ A^-1 = U^-1 L^-1 P

   → Costo computazionale ottimale (O m³) ma si può fare con costo inferiore

   In MATLAB inv(A) esegue questo algorítmo

3. Risoluzione di p systemi lineari aventi la stessa matrice dei coeff.

   A₁x₁ = b₁ → x₁ = A\ b₁ (da Matlab) O (m³)

   A₂x₂ = b₂ → x₂ = A\ b₂ O (m³)

   A₃x₃ = b₃ → x₃ = A b₃ O (m³)

   &O (m³) = costo totale

Esiste algorítmo che calcola questi vettori incogniti con costo minore? Si.

PA=LU ⊕ m³

• P A x_i = Pb_i →

  Ly_i = Pb_i ⊕ U x_i = y_i i=1,2,3

  Funzionallo?