Esercizio analisi numerica

Analisi, algebra lineare, topologia, gruppi, anelli, campi, ...
Rispondi
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Esercizio analisi numerica

Messaggio da Wuque »

Ciao a tutti, c'è qualcuno cos'ì gentile da spiegarmi (come lo spiegherebbe a sua nonna :P ) la soluzione di questo esercizio svolto con Mathematica?
Ringrazio anticipatamente tutti coloro che mi risponderanno.

Vi riporto il testo e la soluzione dell'esercizio.

Risolvere il sistema Ax=b utizizzando il metodo di richardson.
Stabilire empiricamente il miglior valore di "alfa" per la convergenza.

A{aij}ij = 1.N aij=1 i diverso da j
aij= n+1 x= 1 ossia [1,1,1...1]

bi=2n per ogni i.


SOLUZIONE

e = 0.1;
α = 2;
n = 5;
b = Table[1, {i, 1, n}];
A = Table[1 + n*KroneckerDelta[i, j], {i, 1, n}, {j, 1, n}];
x0 = Table[1.5, {i, 1, n}];
r0 = b - A.x0;
x = x0;
r = r0;
i = 0;
While[Norm[r] ≥ ϵ*Norm,
{
i++;
x = x + α*r;
r = b - A.x;
}
]
Print["Numero di iterazioni:"];
Print
Print["Soluzione raggiunta:"]
MatrixForm[x]
Avatar utente
Marco
Site Admin
Messaggi: 1331
Iscritto il: 01 gen 1970, 01:00
Località: IMO '93

Messaggio da Marco »

Wuque ha scritto:A{aij}ij = 1.N aij=1 i diverso da j
aij= n+1 x= 1 ossia [1,1,1...1]

bi=2n per ogni i.
Ciao. Dunque, devo confessare la mia ignoranza sul metodo di Richiardson. Conosco però un po' di Mathematica, e cerco di fare l'esegesi del tuo algoritmino.

Intanto, così come lo hai enunciato, il quesito non è affatto chiaro.
Suppongo che A sia la matrice che vale 1 vuori dalla diagonale e n+1 sulla diagonale principale.
b invece è un vettore con tutte le entrate uguali a 2n.

Codice: Seleziona tutto

e = 0.1;
a = 2;
n = 5;
Qui stai inizializzando un po' di variabili. n capisco cos'è [la dimensione del problema]. e potrebbe essere un parametro di errore che serve per stabilire se la convergenza ha raggiunto un livello soddifacente. Ti faccio notare che e non viene più citato nel seguito, quindi suppongo ci sia un errore. Comunque vedi anche il commento sotteriore.

a non ho idea di cosa sia. Potrebbe essere un "magic number" tipico del metodo? Ho il forte sospetto che qui ci sia puzza d'errore. Tornerò su questo punto poi.

Codice: Seleziona tutto

b = Table[1, {i, 1, n}];
A = Table[1 + n*KroneckerDelta[i, j], {i, 1, n}, {j, 1, n}];
Questo è l'input dei parametri del problema. La cosa che qui non mi torna è che b non è inizializzato a 2n su tutti gli elementi, ma a 1. Forse avresti dovuto scrivere
b = Table[2n, {i, 1, n}];
?

Codice: Seleziona tutto

x0 = Table[1.5, {i, 1, n}];
r0 = b - A.x0;
Da come l'ho capita io, il metodo di Richardson è un metodo basato su iterazioni successive, che cerca di minimizzare il resto (ossia la differenza tra primo e secondo membro del sistema assegnato). x0 è perciò il vettore iniziale e r0 il resto corrispondente. Congetturo che il punto x0, il seme iniziale della successione, possa essere scelto a piacere. In questo caso, il vettore con tutte le entrate pari a 1.5.

Codice: Seleziona tutto

x = x0;
r = r0;
i = 0;
x e r ora sono i valori correnti della successione. i è un indice contatore e x e r quindi dovrebbero contenere sempre x_i e r_i. Giustamente, tali variabili sono inizializzate con x_0 e r_0.

Codice: Seleziona tutto

While[Norm[r] >= ?*Norm[b],
  {
    i++;
    x = x + a*r;
    r = b - A.x;
    }
  ]
Questo passo iterativo è il cuore dell'algoritmo. Nota che nella condizione della clausola While è saltato un simbolo. La cosa più ragionevole da supporre è che ci fosse un e o comunque qualcosa che dipende dal criterio di convergenza. Infatti, non appena il criterio è soddisfatto, la clausola forza l'uscita dal ciclo. A quel punto i contiene il numero di cicli utilizzati per convergere, r dovrebbe risultare sufficientemente piccolo e x essere converso alla soluzione vera.

Ho però dei forti dubbi sulla convergenza di questa roba (e quindi ho il forte sospetto che l'algoritmo giri senza mai uscire).

Ti spiego perché: l'algoritmo calcola la successione
$ $\mathbf x_{i+1} = \mathbf x_i + a \mathbf r_i $
$ $\mathbf r_{i+1} = \mathbf b - A \mathbf x_{i+1} $
L'equazione di r continua a essere quella del resto e mi torna. Non capisco però la riga che aggiorna x. x viene modificato con il resto corrente, moltiplicato per un multiplo scalare [il famoso a che mi puzza].

Ora, ci sono scelte di che possono in qualche modo portarti a qualche sorta di convergenza, ma non certo un valore arbitrario di a.

Nel tuo specifico esempio, nota che hai preso un fortunatissimo valore del punto iniziale x0 e del termine noto b, dato che sono entrambi autovettori del medesimo autospazio di A. L'autovalore è n+2, quindi il tuo problema si risolve banalmente con soluzione $ \frac{1}{n+2}\mathbf b $. Se invece supponi di far girare l'algoritmo, tu hai:
$ $\mathbf x_i = \alpha_i \mathbf b $
$ $\mathbf r_i = (1 - \alpha_i(n+2)) \mathbf b $
$ $\alpha_{i+1} = \alpha_i + 2 ( 1 - \alpha_i(n+2) $.
Affinché il tutto converga, ci attenderemmo che $ \alpha_i $ converga a $ 1/(n+2) $.

Invece i valori saltellano da una parte all'altra, allontanandosi con velocità esponenziale dal punto di equilibrio. Evidentemente, c'è qualcosa che non torna.
[i:2epswnx1]già ambasciatore ufficiale di RM in Londra[/i:2epswnx1]
- - - - -
"Well, master, we're in a fix and no mistake."
Avatar utente
Catraga
Messaggi: 302
Iscritto il: 01 gen 1970, 01:00
Località: Trieste (Univ)

Messaggio da Catraga »

Il codice l'ho scritto io, tutte le varibili inserite sono intrinseche al metodo di Richardson. In effetti il charset di mathematica non viene correttamente riprodotto nel post.
I parametri $ \alpha, \epsilon $ vanno aggiustati in modo che il metodo converga (cosa che avevo lasciato da sistemare a Wuque). 8)
Il vettore b ho sbagliato a digitarlo... :p
fph
Site Admin
Messaggi: 3956
Iscritto il: 01 gen 1970, 01:00
Località: in giro
Contatta:

Messaggio da fph »

Breve riassunto sul metodo di Richardson:
è un metodo iterativo più "ignorante" dei comuni Jacobi e Gauss-Seidel. L'idea è fare per l'appunto quella iterazione,
$ x_{i+1}=x_i+a (b-Ax_i) $
cioè
$ x_{i+1}=(I-aA)x_i+ab $
partendo da un punto x_0 a caso.

Per "cose note" sui metodi iterativi, il resto converge a zero (e quindi x converge alla soluzione) per ogni scelta di x_0 sse la matrice I-aA ha raggio spettrale minore di 1. Se A ha tutte le entrate positive, come nel nostro caso, basta prendere un a abbastanza piccolo e questa condizione è verificata.

ciao,
--federico
[tex]\frac1{\sqrt2}\bigl(\left|\text{loves me}\right\rangle+\left|\text{loves me not}\right\rangle\bigr)[/tex]
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Messaggio da Wuque »

Catraga se solo fossi in grado lo farei... :oops:
Cmq se ho ben capito c'è un errore, qual'è la "soluzione" scritta nel modo corretto?
Me la spiegate? Io sto leggendo manuali, e il mio libro di analisi numerica, ma... più lo leggo più mi sento stupido ;_; Inoltre come fate a scrivere ai nel modo giusto??
Cercate di farmi capire vi prego!

P.S
Per questioni di lavoro risponderò ai futuri post (sperando che ci siano :P ) dal venerdi alla domenica.

Grazie a tutti coloro che mi daranno una mano.
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Messaggio da Wuque »

Possibile che non ci sia nessuno disposto a spiegarmi questo esercizio? Non vado in cerca di una semplice soluzione per un compito o cos'altro, ho solo bisogno di qualcuno che mi spieghi il programma tutto qua. L'esercizio è stato preso da verifiche degli anni passati.
fph
Site Admin
Messaggi: 3956
Iscritto il: 01 gen 1970, 01:00
Località: in giro
Contatta:

Re: Esercizio analisi numerica

Messaggio da fph »

Non ho capito esattamente di cosa hai bisogno ma ci provo...
Wuque ha scritto: e = 0.1;
α = 2;
n = 5;
Inizializza qualche variabile: n è la dimensione della matrice, e credo che sia una specie di "soglia epsilon" anche se non lo ritrovo dopo, a è il parametro che compare nel metodo di Richardson
b = Table[1, {i, 1, n}];
A = Table[1 + n*KroneckerDelta[i, j], {i, 1, n}, {j, 1, n}];
x0 = Table[1.5, {i, 1, n}];
Inizializza matrici e vettori A,b e x_0 come dalle ipotesi del problema. KroneckerDelta è un nome altisonante per la funzione che ritorna 1 se i suoi due argomenti i e j sono uguali, zero altrimenti.
r0 = b - A.x0;
x = x0;
r = r0;
i = 0;
Inizializza un altro po' di roba. x e r sono la stima di x e del "resto" b-Ax che verranno updatate durante l'algoritmo. r_0 è pressoché inutile, perché la prima e la terza riga qui sopra potevano essere condensate in una. i è un contatore che conta il numero di iterazioni fatte (difatti a ogni iterazione qui sotto c'è l'istruzione i++, cioè i=i+1.
While[Norm[r] ≥ ϵ*Norm,


Esegue l'iterazione fintantoché la norma del "resto" (cioè la discrepanza tra Ax e b) non è abbastanza "piccola". Qui credo che si usi il valore "e" che hai inizializzato all'inizio, al posto del carattere sgorbio. Nota che visto che b potrebbe contenere entrate molto grandi, il valore soglia viene moltiplicato per la norma di b.


{
i++;
x = x + α*r;
r = b - A.x;
}
}


Questo è il metodo di Richardson. Se non l'hai mai visto, dovresti riguardarlo su un buon libro di analisi numerica.


Print["Numero di iterazioni:"];
Print
Print["Soluzione raggiunta:"]
MatrixForm[x]


Scrive il valore di i (numero di iterazioni) e quello di x (la stima corrente della soluzione)

spero che ora sia ok?
ciao,
--federico
[tex]\frac1{\sqrt2}\bigl(\left|\text{loves me}\right\rangle+\left|\text{loves me not}\right\rangle\bigr)[/tex]
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Messaggio da Wuque »

Grazie 1000 dell'aiuto fph, qui sotto ti riporo le domande di cosa non mi è ancora chiaro.
Citazione:

b = Table[1, {i, 1, n}];
A = Table[1 + n*KroneckerDelta[i, j], {i, 1, n}, {j, 1, n}];
x0 = Table[1.5, {i, 1, n}];



Inizializza matrici e vettori A,b e x_0 come dalle ipotesi del problema. KroneckerDelta è un nome altisonante per la funzione che ritorna 1 se i suoi due argomenti i e j sono uguali, zero altrimenti.
Le matrici sono del tipo [1,1,1..] ?
Mi potresti spiegare la sintassi di tale scrittura b = Table[1, {i, 1, n}]?
KroneckerDelta è il nome della funzione, ok, ma come fa a ritornare 1 se i e j dal testo dell'esercizio vengono dichiarati diversi?
Citazione:

r0 = b - A.x0;
x = x0;
r = r0;
i = 0;



Inizializza un altro po' di roba. x e r sono la stima di x e del "resto" b-Ax che verranno updatate durante l'algoritmo. r_0 è pressoché inutile, perché la prima e la terza riga qui sopra potevano essere condensate in una. i è un contatore che conta il numero di iterazioni fatte (difatti a ogni iterazione qui sotto c'è l'istruzione i++, cioè i=i+1.
Cosa si intende per stima e resto?
Unendo la prima e la terza cosa avrei dovuto scrivere? r = r0 = b-A.x0 ?

Citazione:

While[Norm[r] ≥ ϵ*Norm,



Esegue l'iterazione fintantoché la norma del "resto" (cioè la discrepanza tra Ax e b) non è abbastanza "piccola". Qui credo che si usi il valore "e" che hai inizializzato all'inizio, al posto del carattere sgorbio. Nota che visto che b potrebbe contenere entrate molto grandi, il valore soglia viene moltiplicato per la norma di b.


Cosa si intende per discrepanza tra Ax e b?
Cos'è la norma di b?


Citazione:

Print["Numero di iterazioni:"];
Print
Print["Soluzione raggiunta:"]
MatrixForm[x]


Scrive il valore di i (numero di iterazioni) e quello di x (la stima corrente della soluzione)


MatrixForm[x] è il comando di mathematica che fa apparire la matrice x?

Infine quando catraga scirsse:
I parametri vanno aggiustati in modo che il metodo converga.
Il vettore b ho sbagliato a digitarlo... :p

Cosa intendeva?
Grazie ancora per l'aiuto
fph
Site Admin
Messaggi: 3956
Iscritto il: 01 gen 1970, 01:00
Località: in giro
Contatta:

Messaggio da fph »

Wuque ha scritto: Le matrici sono del tipo [1,1,1..] ?
Mi potresti spiegare la sintassi di tale scrittura b = Table[1, {i, 1, n}]?
KroneckerDelta è il nome della funzione, ok, ma come fa a ritornare 1 se i e j dal testo dell'esercizio vengono dichiarati diversi?
Non conosco mathematica, ma mi sembra di capire che "Table[f(i), {i, 1, n}]" restituisce un vettore lungo n, che all'entrata i-esima (per i=1..n) contiene il valore f(i).
Quindi b=[1 ... 1 ], A=E+nI (dove I è l'identità e E è la matrice con tutti uni) e x0=[1.5 1.5 ... 1.5].
Nel caso della dichiarazione di A, i e j variano a seconda delle entrate della tabella (=matrice) nel modo in cui sei abituato, quindi a volte KroneckerDelta restituisce 1 e a volte 0.

Cosa si intende per stima e resto?
Unendo la prima e la terza cosa avrei dovuto scrivere? r = r0 = b-A.x0 ?
Basta r=b-A.x0, r0 poi non lo usi più. Stima=vuoi approssimare la soluzione, x è la tua approssimazione corrente. Il resto è b-A.x, che dovrebbe essere uguale a 0 se x è davvero una soluzione.
Cosa si intende per discrepanza tra Ax e b?
La norma del vettore resto b-Ax. Questa non è terminologia standard, è la prima parola che mi è venuta.
Cos'è la norma di b?
ehm... dovresti veramente ripassare qualcosa mi sa. ;-)
Esistono diverse "norme" per un vettore, quella che si usa di solito è la cosiddetta norma-2, cioè
$ \left\Vert v\right\Vert_2=\left(\sum_{i=1}^n v_i^2\right)^{1/2} $
MatrixForm[x] è il comando di mathematica che fa apparire la matrice x?
Boh, non conosco mathematica ma ad occhio direi di sì. Guarda sull'help di mathematica magari.
Infine quando catraga scirsse:
I parametri vanno aggiustati in modo che il metodo converga.
Il vettore b ho sbagliato a digitarlo... :p
Cosa intendeva?
(Domanda: ma il codice che stiamo commentando l'hai scritto tu o lui?)

Che nel testo dell'esercizio b=[2n 2n ... 2n], mentre nel programma si usa un altro vettore, b=[1 1 ... 1]. Probabilmente una distrazione. (tra l'altro nota che basta dividere tutto per 2n e si riconduce una soluzione all'altra).
--federico
[tex]\frac1{\sqrt2}\bigl(\left|\text{loves me}\right\rangle+\left|\text{loves me not}\right\rangle\bigr)[/tex]
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Messaggio da Wuque »

Il codice lo ha scritto catraga, la formula A=E+nI da dove l'hai presa?è quella generica? "E" Nel nostro caso sarebbe b giusto? mentre l'identità qual'è?

Se ho ben capito per rimediare alla diamenticanza di catraga dobrei dichiarare b in questo modo b= Table[2n, {i, 1, n}] giusto?


P.S
Qualcuno mi sa spiegare il perchè scrivendo questo codice su mathematica e cliccando shif+invio non succede nulla??


Fph Grazie molte del tuo aiuto mi sei moooolto Utile! :mrgreen:
fph
Site Admin
Messaggi: 3956
Iscritto il: 01 gen 1970, 01:00
Località: in giro
Contatta:

Messaggio da fph »

Wuque ha scritto:la formula A=E+nI da dove l'hai presa?è quella generica? "E" Nel nostro caso sarebbe b giusto? mentre l'identità qual'è?
E è la matrice con tutte le entrate uguali a 1, E_{ij}=1 per tutti gli i e j. I è la matrice identità (1 sulla diagonale, 0 altrove). Quindi E+nI, se ci pensi bene, è proprio uguale alla matrice A che serve a te (guarda il testo dell'esercizio...).

Il codice "corretto" per inizializzare b è quello giusto.

ciao,
--federico
[tex]\frac1{\sqrt2}\bigl(\left|\text{loves me}\right\rangle+\left|\text{loves me not}\right\rangle\bigr)[/tex]
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Messaggio da Wuque »

OK, un'ultima domanda, ma perchè su Mathematica 5.2 caricando questo esercizio e cliccando shift+invio non succede nulla?? non mi dovrebbe comparire il numero delle iterazioni e la matrice?? So che fph non fa uso di mathematica, quindi se qualcun'altro può rispondermi gli sarei grato.

Può dipendere dalla versione non ancora registrata?
Grazie
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Messaggio da Wuque »

Mi hanno detto che il motivo per cui questo progetto non converge dipende da Alpha che non e' settato correttmante per la risposta, e devo fare dei tentativi, che rientrano anche tra i "to do" dell'esercizio.

Io ho provato con 1, 0.5 e 0.1, ma non mi restituisce nulla, con 0.1 mi da errore. Dov'è che sbaglio? Cosa si intende per "rientra nei "to do" dell'esercizio? Quale ragionamento dovrei seguire per capire in che modo devo settare correttamente alfa?
Wuque
Messaggi: 13
Iscritto il: 07 ago 2006, 11:25

Messaggio da Wuque »

Son passati svariati giorni, ma io non son riuscito a far progressi, c'è qualcuno che può darmi una mano?
Rispondi