Matematica II 15.05.06

Siamo ora nella posizione di poter comprendere meglio cosa genera il comando che Octave propone per risolvere i sistemi lineari

------------------------------------------------------
Solving Linear Equations

To solve the set of linear equations ax = b, use the left division operator, `\':

octave:> a \ b

This is conceptually equivalent to a-1*b, but avoids computing the inverse of a matrix directly.
If the coefficient matrix is singular, Octave will print a warning message and compute a minimum norm solution.
------------------------------------------------------

Consideriamo un sistema lineare

a x = b,   con a matrice di tipo m per n, x in Rn, b in Rm ;

il comando

octave:> a \ b

genera un vettore s in Rn soluzione ai minimi quadrati del sistema, cioe' tale che

|| as - b || <= || at - b ||,   per ogni t in Rn,

e se ci sono piu' soluzioni ai minimi quadrati, fra di esse sceglie quella di norma minima.

Chiaramente, in senso stretto, avremo un vettore di "numeri macchina" che "approssima al meglio" una soluzione ai minimi quadrati del sistema.

------------------------------------------------------
PROBLEMA 1. Supponiamo di avere delle misurazioni dell'umidita' (u) in diversi giorni (g) di un mese:

u1   = 30   in   g1   = 1
u2   = 45   in   g2   = 2
u3   = 53   in   g3   = 13
u4   = 66   in   g4   = 15

Supponiamo che per una qualche ragione si sia interessati a  descrivere con una funzione lineare affine 

u = m g + q,   m, q, costanti reali

la variazione dell'umidita' in funzione del tempo, relativamente ai nostri dati. 
Imponendo che nei vari giorni g1, ..., g4 le umidita' risultino rispettivamente u1, ..., u4, si ha il sistema lineare
di 4 equazioni nelle due incognite m, q, che ci aspettiamo privo di soluzioni

     m + q = 30
  2 m + q = 45
13 m + q = 53
15 m + q = 66 

In realta' cosi' e'. Per ogni scelta di m, q, si avra' un errore non nullo

     m + q - 30 
  2 m + q - 45 
13 m + q - 53 
15 m + q - 66 

Possiamo allora cercare una soluzione ai minimi quadrati, cioe' una coppia m0, q0 che minimizzi la norma dell'errore, o, che e' lo stesso,
che minimizzi il quadrato della norma dell'errore

(m + q - 30)2 + (2 m + q - 45)2 + (13 m + q - 53)2 + (15 m + q - 66)2 

La funzione lineare affine che corrisponde a tali parametri m0, q0, avra' dunque la proprieta' di minimizzare la somma dei quadrati degli scarti fra i valori
che essa assume e i valori osservati, nei giorni 1, 2, 13, 15.

Poiche' le colonne della matrice dei coefficienti del sistema sono linearmente indipendenti, sappiamo che ci sara' una sola soluzione ai minimi quadrati, data
dal coefficiente di Fourier della colonna dei termini noti rispetto alla matrice dei coefficienti.

Di seguito consideriamo il problema e svolgiamo i conti con Octave.

Una volta inseriti i dati, possiamo visualizzarli utilizzando la funzione 'plot';
il comando 'plot(g,u)', dove g e u sono vettori colonna con lo stesso numero di componenti, 
rappresenta nel piano i punti con ascisse le componenti di g ed ordinate le corrispoondenti componenti di u;
il comando 'plot(g,u,'*1')' rappresenta i punti con asterischi del colore '1'.

---con Octave--------------------------------------------------------------------
>> g=[1; 2; 13; 15], u=[30; 45; 53; 66]
g =

        1
        2
       13
       15

u =

       30
       45
       53
       66

>> plot(g,u,'*1')
---------------------------------------------------------------------------------

Introduciamo la matrice dei coefficienti, assegnadola alla variabile A, e la colonna dei termini noti, assegnandola alla variabile b;
e verifichiamo con la funzione 'rref' che non ci sono soluzioni.

---con Octave--------------------------------------------------------------------
>> A=[1 1; 2 1; 13 1; 15 1], b=u
A =

        1        1
        2        1
       13        1
       15        1

b =

       30
       45
       53
       66

>> C=[A b]
C =

        1        1       30
        2        1       45
       13        1       53
       15        1       66

>> rref(C)
ans =

        1        0        0
        0        1        0
        0        0        1
        0        0        0
---------------------------------------------------------------------------------

Calcoliamo la soluzione ai minimi quadrati utilizzando la divisione a sinistra,
e calcoliamo l'errore corrispondente utilizzando la funzione norma 'norm'.

---con Octave--------------------------------------------------------------------
>> s=A\b
s =

  1.0e+01  *

  0.18614
  3.40740

>> norm(A*s - b)
ans = 11.443

---------------------------------------------------------------------------------

Otteniamo lo stesso risultato calcolando il coefficiente di Fourier di b rispetto ad A:

---con Octave--------------------------------------------------------------------
>> s1=(A'*A)^-1*A'*b
s1 =

  1.0e+01  *

  0.18614
  3.40740
---------------------------------------------------------------------------------

Dunque la funzione lineare affine che meglio approssima i nostri dati e'

u = 1.8614 g + 34.0740

Ora inseriamo questa funzione, con nome 'retta'; 
potremo allora comparare i valori osservati 'u' nei giorni 'g', 
con i valori 'retta(g)' che la funzione 'retta' assume nei giorni 'g'.

---con Octave--------------------------------------------------------------------
>> function y=retta(x)
y = 1.8614*x + 34.0740;
endfunction
>> u, retta(g)
u =

       30
       45
       53
       66

ans =

  1.0e+01  *

  3.59354
  3.77968
  5.82722
  6.19950
---------------------------------------------------------------------------------

Possiamo rappresentare sia i dati che il grafico della funzione lineare, cioe'
la retta che meglio li approssima, nel senso dei minimi quadrati;
per dare il senso della continuita', introduciamo un vettore di 'ascisse'
costituito dai numeri della progressione aritmetica di ragione 0.1 che inizia con 0 e termina con 16.

---con Octave--------------------------------------------------------------------
>> w=[0:0.1:16];
>> plot(g,u,'*1',w,retta(w),'2')
---------------------------------------------------------------------------------

------------------------------------------------------
PROBLEMA 2. Consideriamo ancora le stesse misurazioni dell'umidita' (u) in diversi giorni (g) di un mese:

u1   = 30   in   g1   = 1
u2   = 45   in   g2   = 2
u3   = 53   in   g3   = 13
u4   = 66   in   g4   = 15

Supponiamo che per una qualche ragione si sia interessati a  descrivere con una polinomio di grado al piu' due 

u = n g2  + m g + q,   n, m, q, costanti reali

la variazione dell'umidita' in funzione del tempo, relativamente ai nostri dati. 
Imponendo che nei vari giorni g1, ..., g4 le umidita' risultino rispettivamente u1, ..., u4, si ha il sistema lineare
di 4 equazioni nelle tre incognite n, m, q, che ci aspettiamo privo di soluzioni

       n  +    m + q = 30
    4 n +  2 m + q = 45
169 n +13 m + q = 53
225 n +15 m + q = 66 

In realta' cosi' e'. Per ogni scelta di n, m, q, si avra' un errore non nullo

       n  +    m + q - 30
    4 n +  2 m + q - 45
169 n +13 m + q - 53
225 n +15 m + q - 66

Possiamo allora cercare una soluzione ai minimi quadrati, cioe' una terna n0, m0, q0 che minimizzi la norma dell'errore, o, che e' lo stesso,
che minimizzi il quadrato della norma dell'errore

(n  +    m + q - 30)2 + (4 n +  2 m + q - 45)2 + (169 n +13 m + q - 53)2 + (225 n +15 m + q - 66)2 

Il polinomio di secondo grado che corrisponde a tali parametri n0, m0, q0, avra' dunque la proprieta' di minimizzare la somma dei quadrati degli scarti fra i valori
che esso assume e i valori osservati, nei giorni 1, 2, 13, 15.

Poiche' le colonne della matrice dei coefficienti del sistema sono linearmente indipendenti, sappiamo che ci sara' una sola soluzione ai minimi quadrati, data
dal coefficiente di Fourier della colonna dei termini noti rispetto alla matrice dei coefficienti.

Di seguito consideriamo il problema e svolgiamo i conti con Octave.

---con Octave--------------------------------------------------------------------
>> A2 = [1 1 1; 4 2 1; 169 13 1; 225 15 1], b=u
A2 =

        1        1        1
        4        2        1
      169       13        1
      225       15        1

b =

       30
       45
       53
       66

>> C2 = [A2 b]
C2 =

        1        1        1       30
        4        2        1       45
      169       13        1       53
      225       15        1       66

>> rref(C2)
ans =

        1        0        0        0
        0        1        0        0
        0        0        1        0
        0        0        0        1
---con Octave--------------------------------------------------------------------

Abbiamo cosi' riconosciuto che non ci sono soluzioni esatte.
Calcoliamo la soluzione ai minimi quadrati con la divisione a sinistra.
Otteniamo lo stesso risultato calcolando il coefficiente di Fourier.
Calcoliamo anche l'errore corrispondente, che, come ci si poteva aspettare,
e' minore, anche se di poco, di quello trovato al punto precedente.

---con Octave--------------------------------------------------------------------
>> s2 = A2\b
s2 =

  1.0e+01  *

  0.00852
  0.05260
  3.59266

>> (A2'*A2)^-1*A2'*b
ans =

  1.0e+01  *

  0.00852
  0.05260
  3.59266

>> norm(A2*s2 - b)
ans = 11.320
---------------------------------------------------------------------------------

Dunque il polinomio di secondo grado che meglio approssima i nostri dati e'

u = 0.0852 g2 + 0.5260 g + 35.9266

Ora inseriamo questa funzione, con nome 'parab'; 
abbiamo qui usato l'operatore '.^' che, applicato a un vettore, ne eleva a potenza tutte le componenti. 

Rappresentiamo infine i dati, il grafico della funzione lineare, 
e il grafico del polinomio di secondo grado trovati.

---con Octave--------------------------------------------------------------------
>> function y = parab(x)
y = 0.0852*x.^2 + 0.5260*x + 35.9266*x.^0;
endfunction

>> plot(g,u,'*1',w,retta(w),'2',w,parab(w),'3')
---------------------------------------------------------------------------------

=================================================================================
Il problema che abbiamo considerato e' molto frequente ... in realta' in Octave c'e' un'apposita funzione;
qui di seguito ne riportiamo la descrizione e applichiamo la funzione per riottenere le soluzioni trovate.
=================================================================================

>> help polyfit
polyfit is the user-defined function from the file
/usr/share/octave/2.1.42/m/polynomial/polyfit.m

 - Function File: [P, YF] = polyfit (X, Y, N)
     Return the coefficients of a polynomial P(X) of degree N that
     minimizes `sumsq (p(x(i)) - y(i))',  to best fit the data in the
     least squares sense.

     The polynomial coefficients are returned in a row vector if X and
     Y are both row vectors; otherwise, they are returned in a column
     vector.

     If two output arguments are requested, the second contains the
     values of the polynomial for each value of X.

---con Octave--------------------------------------------------------------------
>> polyfit(g,u,1)
ans =

  1.0e+01  *

  0.18614
  3.40740

>> [P1, P1_] = polyfit(g,u,1)
P1 =

  1.0e+01  *

  0.18614
  3.40740

P1_ =

  1.0e+01  *

  3.59354
  3.77969
  5.82724
  6.19953

>> polyfit(g,u,2)
ans =

  1.0e+01  *

  0.00852
  0.05260
  3.59266

>> [P2, P2_] = polyfit(g,u,2)
P2 =

  1.0e+01  *

  0.00852
  0.05260
  3.59266

P2_ =

  1.0e+01  *

  3.65377
  3.73193
  5.71604
  6.29826
---------------------------------------------------------------------------------