Precedente Indice Successivo



Splin


Scopo: valutazione della spline cubica naturale interpolante n punti assegnati.
Specifiche: subroutine splin (p, a, b, c, d, n, m, x, y, ifail)
integer n, m, ifail
real p(n), a(n), b(n), c(n), d(n), x(n), y(n)
Descrizione: Backgroud del problema

Dati n punti distinti xi e n valori corrispondenti yi esiste un'unica spline cubica s(x) interpolante tali punti, che soddisfi le condizioni s''(x1) = s''(x2) = 0. Tale funzione è un polinomio di terzo grado, in ogni intervallo [xi, xi+1] i=1,...,n-1,quindi, per la sua valutazione in un punto x Î [xi, xn] occorre determinare l'intervallo i-esimo contenente x e quindi, noti i coefficienti, calcolare il valore del polinomio relativo ad esso.

Descrizione dell'algoritmo

L'algoritmo, ordinate le ascisse, implementa un metodo di ricerca binaria dell'ascissa x in cui valutare la spline e, determinato l'intervallo di appartenenza e i relativi coefficienti, calcola il polinomio rappresentato nella forma a(x - x)3 + b(x - x)2 + c(x - x) + d con il metodo di Horner.

Raccomandazioni sull'uso

Il numero di punti di interpolazione deve essere maggiore di 2, e le scisse di tali punti devono essere distinte.

Bibliografia: [1], [2]
Parametri di I/O: input

p - vettore di reali, ascisse dei punti di interpolazione
a - vettore di reali, coefficienti dei termini di grado 3
b - vettore di reali, coefficienti dei termini di grado 2
c - vettore di reali, coefficienti dei termini di grado 1
d - vettore di reali, coefficienti dei termini di grado 0
n - intero, numero di punti/coefficienti
x - vettore di reali, punti in cui valutare il polinomio
m - intero, numero di punti di valutazione

output

y - vettore di reali, valori del polinomio
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.
ifail = 15 una delle ascisse di valutazione è esterna all'intervallo di definizione
Routines ausiliarie: nessuna.
Tempo d'esecuzione: La complessità asintotica dell'algoritmo di ricerca binaria è O(log2 n), mentre per la valutazione è necessario un numero di operazioni proporzionale a m.
Memoria richiesta: sono allocati due array di reali in singola precisione di dimensione n, 4 di dimensione n-1 e 2 di dimensione m ,quindi, la complessità asintotica è O(6n+2m).
Accuratezza fornita: dipende dalla precisione del sistema aritmetico.
Esempio di programma chiamante: 
      Program splind
      
c     Programma dimostrativo della subroutine splin

      external splin
      
      integer dim1,dim2
      parameter (dim1=20,dim2=50)
      
      real p(dim1)
      real a(dim1-1),b(dim1-1),c(dim1-1),d(dim1-1)
      real x(dim2),y(dim2)
      integer n,ifail
      character*12 fin1,fin2,fin3,fout
      
      print*                 
      print*,'*** PROGRAMMA DIMOSTRATIVO DELLA SUBROUTINE SPLIN ***'
      print*
      write(*,20)'Numero di punti di interpolazione (max',dim1,')= '
      read(*,*)n
      print*,'Nome del file di input (punti di interpolazione)= '
      read(*,'(A)')fin1
      open(11,FILE=fin1,STATUS='OLD')
      do 1 i=1,n
         read(11,*)p(i)
    1 continue
      close(11)
      
      print*,'Nome del file di input (coefficienti della spline)= '
      read(*,'(A)')fin2
      open(12,FILE=fin2,STATUS='OLD')
      do 2 i=1,n-1
         read(12,*)a(i)
    2 continue
      do 3 i=1,n-1
         read(12,*)b(i)
    3 continue
      do 4 i=1,n-1
         read(12,*)c(i)
    4 continue
      do 5 i=1,n-1
         read(12,*)d(i)
    5 continue
      close(12)
      
      write(*,20)'Numero di punti di tabellazione (max',dim2,')= '
      read(*,*)m
      print*,'Nome del file di input (dominio)= '
      read(*,'(A)')fin3
      open(13,FILE=fin3,STATUS='OLD')
      do 6 i=1,m
         read(13,*)x(i)
    6 continue
      close(13)
      
      call splin(a,b,c,d,n,p,x,m,y,ifail)
      
      if (ifail .eq. 11) stop 'Errore! numero di punti non valido'
      if (ifail .eq. 15) stop 'Errore! ascissa non valida'
      
      print*,'Nome del file di output= '
      read(*,'(A)')fout
      open(14,FILE=fout)
      write(14,'(E14.7)')(y(i),i=1,m)
      close(14)
      
      print*,'Tabellazione della spline :'
      write(*,30)'x=','y='
      write(*,40)(x(i),y(i),i=1,m)
            
   20 format (1X,A,I3,A)   
   30 format (1X,A2,14X,A2)
   40 format (1X,E15.7,E16.7)
          
      End


Precedente Indice Successivo