Precedente Indice Successivo


Cospl

Scopo: calcolo dei coefficienti della spline cubica naturale interpolante n punti assegnati.
Specifiche: subroutine cospl(x, y, n, a, b, c, d, ifail)
integer n, ifail
real x(n), y(x), a(n-1), b(n-1), c(n-1), d(n-1)
Descrizione: Backgroud del problema

Dati n punti distinti xi e n valori corrispondenti yi, una spline cubica, interpolante tali punti è una funzione continua insieme alle sue derivate prima e seconda, e coincide con un polinomio di 3° grado in ogni intervallo [xi, xi+1], i = 1,..., n-1. Il numero di coefficienti da determinare è 4(n-1) = 4n-4, mentre le condizioni di interpolazione e quelle di continuità sono n+3(n-2) = 4n- 6, quindi mancano 2 equazioni affinché il problema sia ben posto. La spline cubica naturale è quella che sfrutta i 2 restanti gradi di libertà imponendo l'annullamento delle derivate seconde negli estremi. Si può verificare con semplici calcoli, che la spline è univocamente determinata dalle sue derivate seconde nei nodi xi, e che queste si ottengono risolvendo un sistema lineare con matrice dei coefficienti di tipo tridiagonale simmetrica a diagonale strettamente dominante e quindi definita positiva.

Descrizione dell'algoritmo

L'algoritmo, dopo l'ordinamento delle ascisse e il controllo dell'assenza di ascisse uguali, costruisce la matrice dei coefficienti e il vettore dei termini noti del sistema avente come soluzione le derivate seconde della spline nei nodi xi, i = 2,..., n-1. Essendo tale matrice simmetrica e definita positiva, essa può essere fattorizzata senza ricorrere al pivoting. La risoluzione del sistema, pertanto, avviene invocando le subroutine lut e slut. Note le derivate seconde si calcolano infine i 4 coefficienti degli n-1 polinomi rappresentati in ogni intervallo [xi, xi+1], i=1,...,n-1 nella forma seguente: a(x-xi)3+b (x-xi)2 + c(x - xi ) + d

Raccomandazioni sull'uso

Il numero di punti di interpolazione deve essere almeno 3 e inoltre le ascisse di tali punti devono essere distinte.

Bibliografia: [1], [2]
Parametri di I/O: input
x - vettore di reali, ascisse dei punti di interpolazione
y - vettore di reali, ordinate dei punti di interpolazione.
n - intero, numero di punti di interpolazione

output
a - vettore di reali, coefficienti dei termini di 3° grado
b - vettore di reali, coefficienti dei termini di 2° grado
c - vettore di reali, coefficienti dei termini di 1° grado
d - vettore di reali, coefficienti dei termini di grado zero
ifail - intero, indicatore d'errore

Indicatori d'errore: Errori gestiti dalla subroutine:
ifail = 0 nessun errore. ifail = 11 il numero dei punti è minore di tre. ifail = 12 sono state immesse almeno due ascisse uguali.
Routines ausiliarie: lut per la fattorizzazione LU di una matrice tridiagonale Slut per la risoluzione di un sistema tridiagonale
Tempo d'esecuzione: La complessità asintotica è O(n).
Memoria richiesta: sono allocati due array di reali in singola precisione di dimensione n e 4 di dimensione n-1 perciò la complessità asintotica è O(6n).
Accuratezza fornita: dipende dal condizionamento della matrice e dalla precisione del sistema aritmetico.
Esempio di programma chiamante: 
       Program cospld 
 c     Programma dimostrativo della subroutine cospl 

       external cospl 
        
       integer dim 
       parameter (dim=20) 
        
       real x(dim),y(dim),a(dim-1),b(dim-1),c(dim-1),d(dim-1) 
       integer n,ifail 
       character*12 fin,fout 
              
       print*,'*** PROGRAMMA DIMOSTRATIVO DELLA SUBROUTINE COSPL ***' 
       print* 
       write(*,10)'Numero di punti di interpolazione (max ',dim,')= ' 
       read(*,*)n 
       print*,'Nome del file di input= ' 
       read(*,'(A)')fin 
       open(11,FILE=fin,STATUS='OLD') 
       do 1 i=1,n 
          read(11,*)x(i) 
     1 continue 
       read(11,*) 
       do 2 i=1,n 
          read(11,*)y(i) 
     2 continue 
       close(11) 
        
       call cospl(x,y,n,a,b,c,d,ifail) 
        
       if (ifail .eq. 11) stop 'Errore! numero di punti non valido' 
       if (ifail .eq. 12) stop 'Errore! ascisse coincidenti' 
        
       print*,'Nome del file di output= ' 
       read(*,'(A)')fout 
       open(12,FILE=fout) 
       
       write(12,'(E14.7)')(a(i),i=1,n-1) 
       write(12,*) 
       write(12,'(E14.7)')(b(i),i=1,n-1) 
       write(12,*) 
       write(12,'(E14.7)')(c(i),i=1,n-1) 
       write(12,*) 
       write(12,'(E14.7)')(d(i),i=1,n-1) 
       close(12)  
                                              
       print*,'Coefficienti della spline= ' 
       write(*,30)'a=','b=','c=','d=' 
       write(*,20)(a(i),b(i),c(i),d(i),i=1,n-1) 
       
    10 format (1X,A,I3,A)    
    20 format (1X,E15.7,E16.7,E16.7,E16.7) 
    30 format (1X,A2,14X,A2,14X,A2,14X,A2)             
       End



Precedente Indice Successivo