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.
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}];
?
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.
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.