Risoluzione di un sistema tridiagonale simmetrico

Con queste considerazioni teoriche abbiamo spostato il problema dell'interpolazione mediante spline cubiche alla risoluzione di un sistema lineare la cui matrice dei coefficienti è tridiagonale e simmetrica.

$$T=\left(\begin{matrix} a_1 & b_2 & & & \\ b_2 & a_2 & \ddots & & \\ & & \ddots & & b_n \\ & & & & \\ & & & b_n & a_n \\ \end{matrix}\right)$$

Per risolvere in maniera efficiente questo problema facciamo le seguenti osservazioni:

• Una matrice tridiagonale $T$ può essere convenientemente rappresentata da due vettori: il primo, $a$, che contenga nelle posizioni $1..n$ gli elementi della diagonale di $T$ e un altro, $b$, che contenga nelle posizioni $2..n$ gli elementi degli elementi al di sopra e al di sotto della diagonale.
• La fattorizzazione più conveniente è quella $LDL^T$

Con la fattorizzazione $LDL^T$, la matrice $T$ viene fattorizzata in due matrici: $L$ che ha gli elementi della diagonale uguali a 1 e la sottodiagonale memorizzata nel vettore $(l_2 \ldots l_n)$. La diagonale della matrice D è il vettore $(d_1 \ldots d_n)$.

Svolgendo formalmente il prodotto $LDL^T$ si ottiene la matrice

$$\left( \begin{array}{ccccc} d_{1} & d_{1}l_{2} & & & 0 \\ & & & &... ... & & &\\ 0 & & & d_{n-1}l_{n} & (d_{n}+l_{n}^{2}d_{n-1}) \\ \end{array} \right)$$ (3.7)

Dalla (3.7) si ricava l'algoritmo per ricavare i vettori $d$ ed $l$.

$$d_1 = a_1$$

$$\left\{\begin{array}{ccl} l_i & = & \frac{b_i}{d_{i-1}} \\ d_i & = & a_i - b_{i}l_i \nonumber \\ \end{array}\right.$$

Questo algoritmo viene seguito nella procedura solveTrid per calcolare la fattorizzazione $LDL^T$ in questo caso particolare.

Successivamente la procedura risolve il sistema $LDL^T m = 6 \delta$ risolvendo tre sistemi successivi

$$\left\{ \begin{array}{ccl} Ly & = & 6 \delta \\ Dz & = & y \\ L^Tm & = & z \\ \end{array} \right.$$

La procedura per calcolare i valori assunti da una spline cubica consiste nel calcolo del vettore $h$ delle differenze tra due ascisse di interpolazione successive, quindi dei vettori che descrivono la matrice dei coefficienti $T$ e del vettore dei termini noti $\delta$. Si richiama quindi la procedura solveTrid per la risoluzione del sistema tridiagonale e si usa l'espressione (3.5) per calcolare i valori di $s_3$ nell'intervallo di tabulazione.

Nei grafici delle figure 3.11-3.13 è illustrato il comportamento dell'interpolazione mediante spline cubiche naturali della funzione di Runge $y=1/(1+x^2)$.

L'interpolazione mediante spline dà in questo caso risultati molto buoni, migliori anche dell'interpolazione polinomiale con ascisse di Cebyshev.

In questo caso la scelta delle ascisse di interpolazione a distanza costante è quella ottimale, perché con questa scelta viene sicuramente soddisfatto il seguente requisito di quasi uniformità che garantisce l'accuratezza della soluzione

$$\frac {\max h_{i}} {\min h_{i}} < k \quad \textrm{con $k$\ costante indipendente da $n$}$$

Questo algoritmo può quindi essere migliorato considerando che le ascisse siano sempre a distanza costante (cioè $h_i=h \; \forall i$) ed evitando in questo modo alcuni calcoli.

La procedura solveTridCost risolve un sistema tridiagonale del tipo

$$\left( \begin{array}{ccccc} 4h & h & & & 0 \\ h & & & & \\ & & & ... ...}{c} \delta_{1}\\ \vdots \\ \\ \\ \\ \vdots \\ \delta_{n-1} \end{array} \right)$$ (3.8)

La procedura splinecubnatCost.m interpola una funzione con una spline cubica naturale calcolata su una partizione $\Delta$ di ascisse equidistanti su $[a,b]$.

Figura 3.11: Interpolazione della funzione di Runge con una spline cubica naturale su una partizione contenente 7 punti nell'intervallo $[-5; 5]$.

Figura 3.12: Interpolazione della funzione di Runge con una spline cubica naturale su una partizione contenente 11 punti nell'intervallo $[-5; 5]$.

Figura 3.13: Interpolazione della funzione di Runge con una spline cubica naturale su una partizione contenente 19 punti nell'intervallo $[-5; 5]$.