Nicola Stella WebSite | Metodi Numerici | Algoritmi di risoluzione per sistemi lineari

La discretizzazione spaziale e temporale delle PDE riconduce le stesse a sistemi di equazioni lineari ricunducibili ad una forma matriciale

$\mathbf{AX} = \mathbf{B}$

dove la forma della matrice dei coefficienti $\mathbf{A}$ induce ad algoritmi specifici di risoluzione, a seconda che sia densa, sparsa, diagonale, tridiagonale,pentadiagonale, etc.

Gli algoritmi risolutivi per i sistemi lineari sono classificabili essenzialmente nelle seguenti categorie

Metodi diretti;
Metodi Iterativi;
Metodi Multigrid.

Metodi Diretti

Per piccoli sistemisi può utilizzare il metodo di Newton. Un metodo diretto abbastanza efficente soprattutto per la soluzione di sistemi tridiagonali è il metodo di eliminazione di Gauss (spesso accompagnato dall'algoritmo di fattorizzazione LU), ma non molto veloce e che in alcuni casi può essere poco efficiente. Un algoritmo efficente per sistemi tridiagonali è l'algoritmo di Thomas. Sistemi che presentano $\mathbf{A}$ in forma tridiagonale a blocchi possono essere risolti con il metodo di eliminazione a blocchi che è un'estensione dell'algoitmo di Thomas.

Metodi Iterativi

Si tratta di procedure che approcciano la soluzione effettuando ripetute iterazioni. In generale si stabilisce inizialmente una soluzione al problema, per poi modificarla producendo una serie di approssimazioni successive fino a convergere alla soluzione esatta entro un margine d'errore fissato.
I metodi Iterativi popolari sono: metodo di Jacobi, metodo di Gauss-Seidel, metodi iterativi a blocchi, metodo SIP (Strongly Implicit Procedure), appartenenti in particolare alla categoria dei metodi iterativi stazionari, mentre tra i metodi iterativi non stazionari citiamo il metodo del Gradiente Coniugato.

Metodi Multigrid

Algoritmo di Gauss e Fattorizzazione LU

Un sistema lineare triangolare

$\left\lbrace \begin{array}{cccc} a_{11}x_{1} & + a_{12}x_{2} & + \cdots & a_{1n}x_{n} = b_{1}\\ &a_{22}x_{2} & + \cdots & + a_{2n}x_{n} = b_{2}\\ \vdots & & & \\ &&&a_{2n}x_{n} = b_{n} \end{array}\right$

può essere posto in forma matriciale

$\left( \begin{array}{cccc} a_{11}&a_{12}&\cdots & a_{1n}\\ 0&a_{22}&\cdots &a_{2n}\\ \vdots &\cdots &\cdots &\vdots\\ 0&0&\cdots &a_{nn} \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array} \right) = \left( \begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array} \right)$

ossia nella forma $AX = B$ .
La soluzione di tale sistema può essere effettuata con il metodo di eliminazione di Gauss:

- Si considera la soluzione n-esima $x_{n} = \frac{b_{n}}{a_{nn}}$
- Per ogni valore k che va da n-1 a 1 si calcola $x_{k} = \frac{1}{a_{kk}}\left( b_{k} - \sum_{j = k+1}^{n}a_{kj}x_{j}\right).$

procedura detta per sostituzione all'indietro. Il caso particolare appena presentato riguarda un sistema la cui matrice dei coefficienti è di tipo triangolare superiore, ma può essere risolto anche per matrici triangolari inferiori con la stessa procedura ma per sostituzione in avanti.

Se invece la matrice dei coefficienti non ha natura triangolare, ma ha la generica forma

$\left( \begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots &\cdots &\cdots &\vdots \\ a_{n1} & a_{n2} &\cdots &a_{nn} \end{array} \right)$

e quindi di un sistema lineare

$\left\lbrace \begin{array}{cccc} a_{11}x_{1} & a_{12}x_{2} & \cdots & a_{1n}x_{n} = b_{1}\\ a_{21}x_{1} & a_{22}x_{2} & \cdots & a_{2n}x_{n} = b_{2}\\ \vdots & \cdots & \cdots & \vdots \\ a_{n1}x_{1} & a_{n2}x_{2} & \cdots &a_{nn}x_{n} = b_{n} \end{array}\right$

si può risolvere il sistema sempre con l'algoritmo di Gauss. Si può infatti fattorizzare la matrice dei coefficienti
$A = LU$
in maniera tale che L sia una matrice diagonale inferiore ed U una matrice diagonale superiore. In questo modo si può risolvere il sistema lineare adottando il metodo di eliminazione di Gauss.

E' importante osservare che la matrice A non deve essere singolare. I coefficienti sulla diagonale $a_{kk}$ devono essere non nulli, ossia $a_{kk} \neq 0$ per ogni valore assunto da k.

Nei problemi computazionali, inoltre, è importante tener presente che valori molto piccoli degli elementi sulla diagonale possono rendere poco efficente la risoluzione numerica. Per ovviare a possibili inconvenienti simili si può introdurre nell'algoritmo una piccola procedura detta pivoting parziale: nel passo in cui viene eliminata la variabile k-esima viene spostata nella posizione k-esima l'equazione in cui il coefficiente della variabile che viene eliminata è massimo; se tale coefficiente massimo ha valore nullo viene riportato un messaggio di errore.

[Top]

Algoritmo di Thomas

L'algoritmo di Thomas viene applicato per sistemi la cui matrice dei coefficienti è tridiagonale. Esso può essere esteso anche per matrici pentadiagonali.

Un sistema

$\left( \begin{array}{ccccccccc} b_{1} & c_{1} & & & & & & &\\ a_{2} & b_{2} & c_{2}& & & & & &\\ & \cdot & \cdot & \cdot & & & & &\\ & & a_{i} & b_{i} & c_{i}& & & &\\ & & & \cdot & \cdot & \cdot & & &\\ & & & & & & a_{n-1} & b_{n-1} & c_{n-1}\\ & & & & & & & a_{n} & b_{n} \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{i}\\ \vdots\\ x_{n-1}\\ x_{n}\\ \end{array} \right) = \left( \begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{i}\\ \vdots\\ d_{n-1}\\ d_{n}\\ \end{array} \right)$

può essere risolto effettuando prima una sostituzione forward

$\left( \begin{array}{ccccccccc} 1 & c^{\prime}_{1} & & & & & & &\\ & 1 & c^{\prime}_{2} & & & & & &\\ & \cdot & \cdot & \cdot & & & & &\\ & & 1 & c^{\prime}_{i} & & & & &\\ & & & \cdot & \cdot & \cdot & & &\\ & & & & & & & 1 & c^{\prime}_{n-1}\\ & & & & & & & & 1 \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{i}\\ \vdots\\ x_{n-1}\\ x_{n}\\ \end{array} \right) = \left( \begin{array}{c} d^{\prime}_{1} \\ d^{\prime}_{2} \\ \vdots \\ d^{\prime}_{i} \cdot\\ \cdot\\ d^{\prime}_{n} \end{array} \right)$

partendo da

$c^{\prime}_{1} = \frac{c_{1}}{b_{1}}$ $d^{\prime}_{1} = \frac{d_{1}}{b_{1}}$

e calcolando ciclicamente

$c^{\prime}_{i} = \frac{c_{i}}{b_{i} - a_{i}c^{\prime}_{i-1}}$ $d^{\prime}_{i} = \frac{d_{i} - a_{i}d^{\prime}_{i-1}}{b_{i} - a_{i}c^{\prime}_{i-1}}$

con i che assume valori da 2 ad n.
Successivamente si effettua una sostituzione all'indietro partendo da

$x_{n} = d^{\prime}_{n}$

e calcolando ciclicamente

$x_{i} = d^{\prime}_{i} - x_{i+1}c^{\prime}_{i}$

con i che assume valori da n-1 a 1.

Affinchè l'algoritmo possa essere applicato è opportuno che la matrice dei coefficienti sia diagonale dominante; in un calcolo numerico bisogna quindi effettuare prima di tutto una verifica, accertandosi cioè che valga la condizione

$|b_{i}| > |a_{i}| + |c_{i}|$

[Top]

Algoritmo di Thomas per sistemi pentadiagonali

Nel caso in cui la matrice dei coefficienti sia pentadiagonale

$\left( \begin{array}{ccccccccc} b_{1} & c_{1} & e_{1} & & & & & &\\ a_{2} & b_{2} & c_{2}& e_{2}& & & & &\\ & \cdot & \cdot & \cdot & & & & &\\ & & f_{i} & a_{i} & b_{i} & c_{i}& e_{i} & &\\ & & & \cdot & \cdot & \cdot & & &\\ & & & & & f_{n-1} & a_{n-1} & b_{n-1} & c_{n-1}\\ & & & & & & f_{n} & a_{n} & b_{n} \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{i}\\ \vdots\\ x_{n-1}\\ x_{n}\\ \end{array} \right) = \left( \begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{i}\\ \vdots\\ d_{n-1}\\ d_{n}\\ \end{array} \right)$

l'algoritmo di Thomas può essere modificato. Si procederà in tre step:

nel primo step si eliminano i termini $f_{i}$

$\left( \begin{array}{ccccccccc} b_{1} & c_{1} & e_{1} & & & & & &\\ a_{2} & b_{2} & c_{2}& e_{2}& & & & &\\ & \cdot & \cdot & \cdot & & & & &\\ & & a^{\prime}_{i} & b^{\prime}_{i} & c^{\prime}_{i}& e^{\prime}_{i} & & &\\ & & & \cdot & \cdot & \cdot & & &\\ & & & & & & a^{\prime}_{n-1} & b^{\prime}_{n-1} & c^{\prime}_{n-1}\\ & & & & & & & a^{\prime}_{n} & b^{\prime}_{n} \end{array} \right) \left( \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{i}\\ \vdots\\ x_{n-1}\\ x_{n}\\ \end{array} \right) = \left( \begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d^{\prime}_{i}\\ \vdots\\ d^{\prime}_{n-1}\\ d^{\prime}_{n}\\ \end{array} \right)$

partendo dalla riga i=3 calcolando ciclicamente

$a^{\prime}_{i} = a_{i} - \frac{f_{i}b^{\prime}_{i-1}}{a^{\prime}_{i-1}}$

$b^{\prime}_{i} = b_{i} - \frac{f_{i}c^{\prime}_{i-1}}{a^{\prime}_{i-1}}$

$c^{\prime}_{i} = c_{i} - \frac{f_{i}e^{\prime}_{i-1}}{a^{\prime}_{i-1}}$

$d^{\prime}_{i} = d_{i} - \frac{f_{i}d^{\prime}_{i-1}}{a^{\prime}_{i-1}}$

$e^{\prime}_{i} = e_{i}$ ;

Nel secondo step si procede con la sostituzione in avanti, con

$c^{\prime\prime}_{i} = \frac{c^{\prime}_{i} - a^{\prime}_{i}e^{\prime\prime}_{i-1}}{b^{\prime}_{i} - a^{\prime}_{i}c^{\prime\prime}_{i-1}}$

$d^{\prime\prime}_{i} = \frac{d^{\prime}_{i} - a^{\prime}_{i}d^{\prime\prime}_{i-1}}{b^{\prime}_{i} - a^{\prime}_{i}c^{\prime\prime}_{i-1}}$

$e^{\prime\prime}_{i} = \frac{e^{\prime}_{i}}{b^{\prime}_{i} - a^{\prime}_{i}c^{\prime\prime}_{i-1}}$

si effettua quindi l'ultimo step, ossia la sostituzione all'indietro con

$x_{i} = d^{\prime\prime}_{i} - c^{\prime\prime}_{i}x_{i+1} - e^{\prime\prime}_{i}x_{i+2}$

[Top]

Metodo di eliminazione a blocchi

L'algoritmo di Thomas viene usato per la soluzione di una singola equazione che, discretizzata, genera un sistema tridiagonale scalare. Tuttavia in molti problemi il sistema fisico da studiare è composto da un sistema di equazioni da risolvere ed esso una volta discretizzato genera un sistema tridiagonale a blocchi. Per meglio intenderci la matrice dei coefficienti in quest'ultimo caso non è composto da elementi scalari, ma ognuno degli elementi è una sottomatrice $m\times m$ mentre le matrici colonna delle incognite e dei termini noti sono caratterizzate da elementi che corrispondono a sottovettori di dimensione $m$

$\left( \begin{array}{ccccccccc} \tilde{b}_{1} & \tilde{c}_{1} & & & & & & &\\ \tilde{a}_{2} & \tilde{b}_{2} & \tilde{c}_{2}& & & & & &\\ & \cdot & \cdot & \cdot & & & & &\\ & & \tilde{a}_{i} & \tilde{b}_{i} & \tilde{c}_{i}& & & &\\ & & & \cdot & \cdot & \cdot & & &\\ & & & & & & \tilde{a}_{n-1} & \tilde{b}_{n-1} & \tilde{c}_{n-1}\\ & & & & & & & \tilde{a}_{n} & \tilde{b}_{n} \end{array} \right) \left( \begin{array}{c} \tilde{x}_{1} \\ \tilde{x}_{2} \\ \vdots \\ \tilde{x}_{i}\\ \vdots\\ \tilde{x}_{n-1}\\ \tilde{x}_{n}\\ \end{array} \right) = \left( \begin{array}{c} \tilde{d}_{1} \\ \tilde{d}_{2} \\ \vdots \\ \tilde{d}_{i}\\ \vdots\\ \tilde{d}_{n-1}\\ \tilde{d}_{n}\\ \end{array} \right).$

L'algoritmo di Thomas può essere esteso anche a questo caso; risolviamo il sistema sue step.

Trasformiamo la matrice originaria in una matrice tridiagonale superiore eliminando le matrici $\tilde{a}_{i}$

$\tilde{c}^{\prime}_{1} = \left(\tilde{b}_{1}\right)^{-1}\tilde{c}_{1},$ $\tilde{d}^{\prime}_{1} = \left(\tilde{b}_{1}\right)^{-1}\tilde{d}_{1},$

$\tilde{b}^{\prime}_{i} = \tilde{b}_{i} - \tilde{a}_{i}\tilde{c}^{\prime}_{i-1}$ $\tilde{c}^{\prime}_{i} = \left(\tilde{b}^{\prime}_{i}\right)^{-1}\tilde{c}_{i},$ $\tilde{d}^{\prime}_{i} = \left(\tilde{b}^{\prime}_{i}\right)^{-1}\left[\tilde{d}_{i} - \tilde{a}_{i}\tilde{d}^{\prime}_{i-1}\right].$

Alla fine si otterrà un sistema in cui $\tilde{c}^{\prime}_{i}$ rimpiazza $\tilde{c}_{i},$ $\tilde{d}^{\prime}_{i}$ rimpiazza $\tilde{d}_{i},$ e la matrice unità $I$ rimpiazza $\tilde{b}_{1}.$

Nel secondo step si effettua una sostituzione all'indietro

$\tilde{x}_{n} = \tilde{d}^{\prime}_{n}$

$\tilde{x}_{i} = \tilde{d}^{\prime}_{i} - \tilde{c}^{\prime}_{i}\tilde{x}_{i+1}$

[Top]

Metodo di Gauss-Seidel

L'algoritmo di Gauss-Siedel è un algoritmo iterativo per la risoluzione di sistemi lineari per i quali i coefficienti sulla diagonale della matrice dei coefficienti hanno un valore in modulo mggiore degli elementi non fiagonali.

Dato un sistema $AX=B$ deve valere $|a_{ii}|>\sum_{i\neq j}|a_{ij}|.$

Prima di tutto riscriviamo la matrice $A$ come la somma algebrica di una matrice diagonale $D$ , una matrice triangolare inferiore $G$ e una matrice triangolare superiore $S$

$A=D - G - S;$

la soluzione del sistema sarà quindi data dall'iterazione della seguente espressione

$x^{n+1} = \left(D-G\right)^{-1}\left(Sx^{n} + B\right)$

ossia, in termini di componenti

$x^{n+1}_{i} = \frac{1}{a_{ii}}\left(b_{i}-\sum_{j\prec i}a_{ij}x^{n+1}_{j}-\sum_{j\succ i}a_{ij}x^{n}_{j}\right)$

L'implementazione si effettua fissando inizialmente una soluzione $\overline{x}$ a cui deve convergere la soluzione numerica; quindi si itera la precedente per ogni step k e al termine del ciclo si verifica se è stata raggiunta la convergenza ad $\overline{x}$ con una certa tolleranza; se la condizione non viene soddisfatta il ciclo riparte con una nuova $\overline{x}.$

[Top]

Altri Articoli...

Pagina 1 di 2

«
Inizio
Prec.
1
2
Succ.
Fine
»

News by PhysOrg.com

PhysOrg.com - latest science and technology news stories

Physorg.com internet news portal provides the latest news on science including: Physics, Nanotechnology, Life Sciences, Space Science, Earth Science, Environment, Health and Medicine.

Busted BP well no longer 'threat' to Gulf: US official
The Macondo well, which spilled an estimated 4.9 million barrels of oil into the Gulf...
China launches communications satellite
China successfully launched a satellite for radio and television broadcasts early Sunday, state media said.
Over one million join Apple's music social network Ping
More than one million users joined Apple's new music-focused social network Ping in the 48...
Kinect's Israeli partner sees a remoteless world
(AP) -- Inon Beracha envisions a world where your movements control the gadgets and...
Craigslist removes adult services section
(AP) -- Craigslist appears to have surrendered in a legal fight over erotic ads posted...
French science vessel to start second leg of climate voyage
The French yacht Tara leaves Sunday on the second leg of a three-year voyage across...
Cholera kills 41 in Central African nation of Chad
(AP) -- Health officials in Chad say an outbreak of cholera in the Central...
Texas opens inquiry into Google search rankings (Update)
(AP) -- Google Inc.'s methods for recommending websites are being reviewed by Texas' attorney...
Google to pay 8.5 million dollars to settle Buzz case
Google has agreed to pay 8.5 million dollars (US) to settle a privacy lawsuit over...
Interview: Wikipedia founder bullish on news
(AP) -- Wikipedia founder Jimmy Wales believes relief may be in sight for the...
Texting on the rise among US adults: Pew survey
More American adults are texting but they are not tapping out nearly as many messages...
German court rules against YouTube over copyright
(AP) -- A German court ruled Friday that Google Inc.'s subsidiary YouTube LLC must...
What's causing life-threatening blood clots in brain surgery patients?
One of the most severe complications of brain surgery is a life-threatening blood clot in...
Moonstruck primates: Owl monkeys need moonlight as much as a biological clock for nocturnal activity
An international collaboration led by a University of Pennsylvania anthropologist has shown that environmental factors,...
NASA satellite and International Space Station catch Earl weakening
NASA satellites and the International Space Station are keeping eyes on Hurricane Earl as it...