Calcolo numerico

Programmazione in Fortran 77.

Interpolazione.

 In questa rubrica presento una completa serie di algoritmi implementati in Digital Visual Fortran 90 ( totalmente compatibile con il FORTRAN 77 ).

Nel 2005 questi programmi sono serviti agli studenti dei corsi serali per studenti lavoratori per superare l'esame universitario "Calcolo numerico e programmazione" presso la facoltą di ingegneria di Padova ( corso della professoressa Michela Redivo Zaglia per allievi aereospaziali).

I programmi sono testati e funzionanti ma essi non assicurano il superamento della prova d'esame senza un approfondito studio teorico nel testo "Calcolo numerico"  dell'autrice Redivo Zaglia.

Attenzione: Nell'anno in cui questi programmi sono stati creati sono stati prodotti in sede d'esame (quindi vanno ben memorizzati nella loro essenza e non solo a memoria) usando il compilatore in linea del sistema operativo Linux installato in aula Taliercio della facoltą di ingegneria in via Belzoni  7 (PD). E' indispensabile al fine di superare la prova d'esame, sapere compilare in ambiente Linux ed effettuare la visualizzazione degli eventuali grafici con i tools di default presenti in quel sistema operativo.

I programmi  name.f  posso essere tagliati ed incollati nell'editor di programmazione per un corretto funzionamento.   Ove sia necessario bisognerą far trovare nel direttorio di compilazione i file di testo contenenti l'insieme di input "formattato" ovvero con righe e colonne al posto in cui l'algoritmo le aspetta.  Nella stessa directory troveremo il file testuale di output.

Marco Gottardo.

Interpolazione con il metodo di Lagrange
lagrange.f lagrange.in & .out
C====================================================================
C Programma per in calcolo del Polinomio interpolatore
C su piu' punti utilizzando la FORMULA DI LAGRANGE
C====================================================================
C
PROGRAM pollag
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxn, maxpt
INTEGER i, is, k
PARAMETER (maxn=50, maxpt=300)
REAL*8 x(0:maxn), y(0:maxn), xs(0:maxpt), ys(0:maxpt)
REAL*8 f, a, b, h
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, 'Formula di LAGRANGE'
print *
print *, '1 = Lettura nodi da file'
print *, '2 = Nodi equidistanti'
print *
print *, 'Inserire la scelta'
read *, is
if (is.eq.1) then
print *, ' Inserisci il numero di nodi (n+1)'
read *, n
n = n-1
print *, 'Dammi il nome del file di ingresso contenente i nodi'
read*, nomefilein
open (unit=9,file=nomefilein)
do i = 0, n
read (9,*) x(i)
end do
close (9)
else
print *, 'Dammi gli estremi a e b dell''intervallo '
print *, 'da discretizzare'
read*, a, b
print *, ' Inserisci il numero di intervalli n'
read *, n
h = (b-a)/dble(n)
x(0) = a
do i = 1, n
x(i) = x(0) + h * i
end do
end if
print *, '1 = Lettura valori da interpolare da file'
print *, '2 = Valori equidistanti'
print *
print *, 'Inserire la scelta'
read *, is
if (is.eq.1) then
print *, ' Inserisci il numero di punti (m+1) per valutare'
print *, ' l''approssimazione del polinomio interpolatore'
read *, m
m = m-1
print *, 'Dammi il nome del file contenente i valori'
read*, nomefilein
open (unit=9,file=nomefilein)
do i = 0, m
read (9,*) xs(i)
end do
close (9)
else
print *, 'Stesso intervallo: 1=SI 2=NO'
read *, is
if (is.eq.2) then
print *, 'Dammi gli estremi a e b dell''intervallo '
print *, 'da discretizzare'
read*, a, b
else
a = x(0)
b = x(n)
end if
print *, ' Inserisci il numero di intervalli m'
read *, m
h = (b-a)/dble(m)
xs(0) = a
do i = 1, m
xs(i) = xs(0) + h * i
end do
end if
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=10,file=nomefileout)
C
C Calcolo della funzione sui nodi
C
do i = 0, n
y(i) = f(x(i))
end do
C
C Ciclo esterno per i punti su cui valutare il polinomio
C interpolatore
C
do k = 0, m
C
C Esecuzione della valutazione con Formula di Lagrange
C
CALL lagrange (x, y, xs(k), ys(k), n)
end do
C
C Uscita dati e risultati
C
write(*,*)
write(*,*) ' Scrittura dei nodi e valori'
do i = 0, n
write (*,*) i, x(i), y(i)
end do
write(*,*)
write(*,*) ' Scrittura dei punti interpolati'
do k = 0, m
write (*,*) k, xs(k), ys(k)
write (10,'(2g20.15)') xs(k), ys(k)
end do
stop
end


C====================================================================
C Subroutine per la valutazione del polinomio di interpolazione
C in un punto con la FORMULA DI LAGRANGE
C====================================================================
SUBROUTINE lagrange (x, y, xs, ys, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore contenente i nodi (0:n)
C y vettore contenente i valori della funzione nei nodi (0:n)
C xs valore di cui effettuare la valutazione
C n il numero dei nodi e' (n+1)
C
C Variabili USCITA:
C ys approssimazione calcolata

C
C====================
C Parte dichiarativa
C====================
INTEGER n
INTEGER k, i, r, j
REAL*8 x(0:*), y(0:*), L(0:n)
REAL*8 xs, ys
C====================
C Parte esecutiva
C====================

do i = 0, n
L(i) = 1.0d0
do j = 0, i-1
L(i) = L(i)*(xs-x(j))/(x(i)-x(j))
end do
do j = i+1, n
L(i) = L(i)*(xs-x(j))/(x(i)-x(j))
end do
end do
ys = 0.0d0
do i = 0, n
ys = ys + y(i)*L(i)
end do
return
end

C====================================================================
C Funzione assegnata
C====================================================================
real*8 function f(x)
implicit none
real*8 x
f=dsqrt(x)
return
end
 
lagrange.in
2.0d0
2.1d0
2.2d0
2.3d0
2.4d0
 

 

 

 

 

 

 

 
lagrange.out
 Nodi e valori della funzione da interpolare
0 2.00000000000000 1.41421356237310
1 2.10000000000000 1.44913767461894
2 2.20000000000000 1.48323969741913
3 2.30000000000000 1.51657508881031
4 2.40000000000000 1.54919333848297
Scrittura del punto interpolato
nodo (xs) val.polinomio int. (ys)
2.05000000000000 1.43178207894254
 

 

Programma di interpolazione con il metodo di Neville-Aitken.

Interpolazione con il metodo di Neville-Aitken
Neville-Aitken.f  
C====================================================================
C Programma per in calcolo del Polinomio interpolatore
C su piu' punti utilizzando la FORMULA DI NEVILLE-AITKEN
C===================================================================
C
PROGRAM NevilleAitken
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxn, maxpt
INTEGER i, is, k
PARAMETER (maxn=50, maxpt=300)
REAL*8 x(0:maxn), y(0:maxn), xs(0:maxpt), ys(0:maxpt)
REAL*8 T(0:maxn,0:maxn)
REAL*8 f, a, b, h
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, 'Formula di NEVILLE-AITKEN'
print *
print *, '1 = Lettura nodi da file'
print *, '2 = Nodi equidistanti'
print *
print *, 'Inserire la scelta'
read *, is
if (is.eq.1) then
print *, ' Inserisci il numero di nodi (n+1)'
read *, n
n = n-1
print *, 'Dammi il nome del file di ingresso contenente i nodi'
read*, nomefilein
open (unit=9,file=nomefilein)
do i = 0, n
read (9,*) x(i)
end do
close (9)
else
print *, 'Dammi gli estremi a e b dell''intervallo '
print *, 'da discretizzare'
read*, a, b
print *, ' Inserisci il numero di intervalli n'
read *, n
h = (b-a)/dble(n)
x(0) = a
do i = 1, n
x(i) = x(0) + h * i
end do
end if
print *, '1 = Lettura valori da interpolare da file'
print *, '2 = Valori equidistanti'
print *
print *, 'Inserire la scelta'
read *, is
if (is.eq.1) then
print *, ' Inserisci il numero di punti (m+1) per valutare'
print *, ' l''approssimazione del polinomio interpolatore'
read *, m
m = m-1
print *, 'Dammi il nome del file contenente i valori'
read*, nomefilein
open (unit=9,file=nomefilein)
do i = 0, m
read (9,*) xs(i)
end do
close (9)
else
print *, 'Stesso intervallo: 1=SI 2=NO'
read *, is
if (is.eq.2) then
print *, 'Dammi gli estremi a e b dell''intervallo '
print *, 'da discretizzare'
read*, a, b
else
a = x(0)
b = x(n)
end if
print *, ' Inserisci il numero di intervalli m'
read *, m
h = (b-a)/dble(m)
xs(0) = a
do i = 1, m
xs(i) = xs(0) + h * i
end do
end if
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=10,file=nomefileout)
C
C Calcolo della funzione sui nodi
C
do i = 0, n
y(i) = f(x(i))
end do
C
C Ciclo esterno per i punti su cui valutare il polinomio
C interpolatore
C
do k = 0, m
C
C Esecuzione della valutazione con Formula di Neville-Aitken
C
CALL neville (x, y, xs(k), ys(k), n)
end do
C
C Uscita dati e risultati
C
write(*,*)
write(*,*) ' Scrittura dei nodi e valori'
do i = 0, n
write (*,*) i, x(i), y(i)
end do
write(*,*)
write(*,*) ' Scrittura dei punti interpolati'
do k = 0, m
write (*,*) k, xs(k), ys(k)
write (10,'(2g20.15)') xs(k), ys(k)
end do
stop
end


C====================================================================
C Subroutine per la valutazione del polinomio di interpolazione
C in un punto con la FORMULA DI NEVILLE-AITKEN
C====================================================================
SUBROUTINE neville (x, y, xs, ys, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore contenente i nodi (0:n)
C y vettore contenente i valori della funzione nei nodi (0:n)
C xs valore di cui effettuare la valutazione
C n il numero dei nodi e' (n+1)
C
C Variabili USCITA:
C ys approssimazione calcolata

C
C====================
C Parte dichiarativa
C====================
INTEGER n
INTEGER k, i
REAL*8 x(0:*), y(0:*), T(0:n,0:n)
REAL*8 xs, ys
C====================
C Parte esecutiva
C====================
do i = 0, n
T(i,0) = y(i)
end do
do k = 0, n
do i = 0, n-k
T(i,k+1)=T(i+1,k)-(((x(i+k+1)-xs)/
* (x(i+k+1)-x(i)))*(T(i+1,k)-T(i,k)))
enddo
end do
ys = 0.0d0
ys = T(0,n)
return
end

C====================================================================
C Funzione assegnata
C====================================================================
real*8 function f(x)
implicit none
real*8 x
f=dsqrt(x)
return
end
 		  

Interpolazione con il metodo di Newton con il metodo delle differenze divise.

Interpolazione con la formula di Newton alle differenze divise.
newton_diff_div.f diffdiv.out
C====================================================================
C Programma per in calcolo del Polinomio interpolatore
C su piu' punti utilizzando la FORMULA DI NEWTON DIFFERENZE DIVISE
C===================================================================
C
PROGRAM Newtondiffdiv
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxn, maxpt, maxdim
INTEGER i, is, k
PARAMETER (maxn=50, maxpt=300)
REAL*8 x(0:maxn), y(0:maxn), xs(0:maxpt), ys(0:maxpt)
c REAL*8 C(0:maxn,0:maxn)
REAL*8 f, a, b, h
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, 'Formula di NEWTON DIFF DIVISE'
print *
print *, '1 = Lettura nodi da file'
print *, '2 = Nodi equidistanti'
print *
print *, 'Inserire la scelta'
read *, is
if (is.eq.1) then
print *, ' Inserisci il numero di nodi (n+1)'
read *, n
n = n-1
print *, 'Dammi il nome del file di ingresso contenente i nodi'
read*, nomefilein
open (unit=9,file=nomefilein)
do i = 0, n
read (9,*) x(i)
end do
close (9)
else
print *, 'Dammi gli estremi a e b dell''intervallo '
print *, 'da discretizzare'
read*, a, b
print *, ' Inserisci il numero di intervalli n'
read *, n
h = (b-a)/dble(n)
x(0) = a
do i = 1, n
x(i) = x(0) + h * i
end do
end if
print *, '1 = Lettura valori da interpolare da file'
print *, '2 = Valori equidistanti'
print *
print *, 'Inserire la scelta'
read *, is
if (is.eq.1) then
print *, ' Inserisci il numero di punti (m+1) per valutare'
print *, ' l''approssimazione del polinomio interpolatore'
read *, m
m = m-1
print *, 'Dammi il nome del file contenente i valori'
read*, nomefilein
open (unit=9,file=nomefilein)
do i = 0, m
read (9,*) xs(i)
end do
close (9)
else
print *, 'Stesso intervallo: 1=SI 2=NO'
read *, is
if (is.eq.2) then
print *, 'Dammi gli estremi a e b dell''intervallo '
print *, 'da discretizzare'
read*, a, b
else
a = x(0)
b = x(n)
end if
print *, ' Inserisci il numero di intervalli m'
read *, m
h = (b-a)/dble(m)
xs(0) = a
do i = 1, m
xs(i) = xs(0) + h * i
end do
end if
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=10,file=nomefileout)
C
C Calcolo della funzione sui nodi
C
do i = 0, n
y(i) = f(x(i))
end do
C
C Ciclo esterno per i punti su cui valutare il polinomio
C interpolatore
C
do k = 0, m
C
C Esecuzione della valutazione con Formula di Newton differenze divise
C
CALL diffdiv (x, y, xs(k), ys(k), n)
end do
C
C Uscita dati e risultati
C
write(*,*)
write(*,*) ' Scrittura dei nodi e valori'
do i = 0, n
write (*,*) i, x(i), y(i)
end do
write(*,*)
write(*,*) ' Scrittura dei punti interpolati'
do k = 0, m
write (*,*) k, xs(k), ys(k)
write (10,'(2g20.15)') xs(k), ys(k)
end do
c Call scrivimat (C, n-1, n-1, maxdim, 10)
stop
end


C====================================================================
C Subroutine per la valutazione del polinomio di interpolazione
C in un punto con la FORMULA DI NEWTON DIFFERENZE FINITE
C====================================================================
SUBROUTINE diffdiv (x, y, xs, ys, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore contenente i nodi (0:n)
C y vettore contenente i valori della funzione nei nodi (0:n)
C xs valore di cui effettuare la valutazione
C n il numero dei nodi e' (n+1)
C
C Variabili USCITA:
C ys approssimazione calcolata

C
C====================
C Parte dichiarativa
C====================
INTEGER n
INTEGER k, i, j
REAL*8 x(0:*), y(0:*), C(0:n,0:n)
REAL*8 xs, ys, u
C====================
C Parte esecutiva
C====================

do i = 0, n
C(i,0) = y(i)
enddo
do k = 1, n
do i = 0, n-k
C(i,k)=(C(i+1,k-1)-C(i,k-1))/(x(i+k)-x(i))
enddo
end do
u= C(0,n)
do j = n-1, 0, -1
u= u*(xs-x(j))+C(0,j)
enddo
ys = u
return
end

C====================================================================
C Funzione assegnata
C====================================================================
real*8 function f(x)
implicit none
real*8 x
f=dsqrt(x)
return
end

C====================================================================
C Subroutine per la SCRITTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE scrivimat (C, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve scrivere la matrice
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 C(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
write (io,*) (C(i,j), j=1,m)
end do
return
end		  1.41421356237310 0.349241122458488 -4.110447228300420E-002
9.243006108147494E-003 -2.486810899326778E-003
1.44913767461894 0.341020228001887 -3.833157045055995E-002
8.248281748416783E-003 0.000000000000000E+000
1.48323969741913 0.333353913911775 -3.585708592603492E-002
0.000000000000000E+000 0.000000000000000E+000
1.51657508881031 0.326182496726568 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000
1.54919333848297 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000
9.864708832979189E-003
-4.258417860795108E-002
0.351370331388885
1.43178207894254
1.41421356237310 0.349241122458488 -4.110447228300420E-002
9.243006108147494E-003 -2.486810899326778E-003
1.44913767461894 0.341020228001887 -3.833157045055995E-002
8.248281748416783E-003 0.000000000000000E+000
1.48323969741913 0.333353913911775 -3.585708592603492E-002
0.000000000000000E+000 0.000000000000000E+000
1.51657508881031 0.326182496726568 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000
1.54919333848297 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000
9.616027743046511E-003
-4.158527367015653E-002
0.347161858774980
1.46628784118934
1.41421356237310 0.349241122458488 -4.110447228300420E-002
9.243006108147494E-003 -2.486810899326778E-003
1.44913767461894 0.341020228001887 -3.833157045055995E-002
8.248281748416783E-003 0.000000000000000E+000
1.48323969741913 0.333353913911775 -3.585708592603492E-002
0.000000000000000E+000 0.000000000000000E+000
1.51657508881031 0.326182496726568 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000
1.54919333848297 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000
8.869984473248476E-003
-3.888697616469208E-002
0.335630680800845
1.56524736873348
2.05000000000000 1.43178207894254
2.15000000000000 1.46628784118934
2.45000000000000 1.56524736873348

Programma per in calcolo del Polinomio interpolatore su un punto utilizzando lo Schema di Neville-Aitken con aggiunta di punti.

Interpolazione con lo schema di Neville_Aitken.
Neville_Aitken.F  
C=========================================================
C Programma per in calcolo del Polinomio interpolatore su un punto
C utilizzando lo Schema di Neville-Aitken con aggiunta di punti
C=========================================================
C
PROGRAM nevaitkena
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxn, maxpt
INTEGER i, is, k ,it
PARAMETER (maxn=50, maxpt=300)
REAL*8 x(0:maxn), y(0:maxn),t(0:maxn,0:maxn)
REAL*8 f, a, b, h,xs,ys
CHARACTER*20 nomefilein, nomefileout,sn
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, 'Formula di LAGRANGE'
print *
print *, '1 = Lettura nodi da file'
print *, '2 = Nodi equidistanti'
print *
print *, 'Inserire la scelta'
read *, is
if (is.eq.1) then
c n= numero di intervalli ---> n+1= numero punti
print *, ' Inserisci il numero di nodi (n+1)'
read *, n
n = n-1
print *, 'Dammi il nome del file di ingresso contenente i nodi'
read*, nomefilein
open (unit=9,file=nomefilein)
do i = 0, n
read (9,*) x(i)
end do
close (9)
else
print *, 'Dammi gli estremi a e b dell''intervallo '
print *, 'da discretizzare'
read*, a, b
print *, ' Inserisci il numero di intervalli n'
read *, n
h = (b-a)/dble(n)
x(0) = a
do i = 1, n
x(i) = x(0) + h * i
end do
end if
print *, 'dammi il valore x stellato xs in cui determinare'
print *, 'il valore del Polinomio interpolante'
read*, xs
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=10,file=nomefileout)
C
C Calcolo della funzione sui nodi
C
do i = 0, n
y(i) = f(x(i))
end do

C
C Esecuzione della valutazione con Formula di Neville-Aitken
CALL nevait (x, y, T ,maxn, xs, ys, n)
C
C Uscita dati e risultati
C
write(*,*)
write(10,*) 'Nodi e valori della funzione da interpolare'
do i = 0, n
write (10,*) i, x(i), y(i)
end do
write(10,*) ' Scrittura del punto interpolato'
write (10,'(a,a)') 'nodo (xs)',' val.polinomio int. (ys)'
write (10,'(2g20.15)') xs, ys
write(*,*) 'Nodi e valori della funzione da interpolare'
do i = 0, n
write (*,*) i, x(i), y(i)
end do
write(*,*)
write(*,*) ' Scrittura del punto interpolato'
write (*,'(a,a)') 'nodo (xs)',' val. polinomio int. (ys)'
write (*,'(2g20.15)') xs, ys
C
C EVENTUALE AGGIUNTA DI ALTRI NODI
C per meglio approssimare funzione
C
print *, 'Vuoi aggiungere un nodo (s/n) ? '
read *, sn
do while (sn.eq.'s')
print *, 'dammi il valore x aggiunto '
read*, x(n+1)

C
C Calcolo della funzione sul nodo aggiunto
C

c incremento n per poter costruire la nuova diagonale montante
n=n+1
y(n) = f(x(n))
CALL nagg (x, y, T, maxn, xs, ys, n)

write(*,*) ' Scrittura del punto interpolato'
write (*,'(a,a)') 'nodo (xs)',' val aggiornato. (ys)'
write (*,'(2g20.15)') xs, ys
print *, 'vuoi continuare s/n? '
read *, sn
enddo
close(10)
stop
end


C====================================================================
C Subroutine per la valutazione del polinomio di interpolazione
C in un punto con la Schema di Neville-Aitken
C====================================================================
SUBROUTINE nevait(x, y, T , maxn, xs, ys, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore contenente i nodi (0:n)
C y vettore contenente i valori della funzione nei nodi (0:n)
C xs valore di cui effettuare la valutazione
C n il numero dei nodi e' (n+1)
C
C Variabili USCITA:
C ys approssimazione calcolata

C
C====================
C Parte dichiarativa
C====================
INTEGER n,maxn
INTEGER k, i, r, j
REAL*8 x(0:*), y(0:*), T(0:maxn,0:*)
REAL*8 xs, ys
C====================
C Parte esecutiva
C====================

do i = 0, n
T(i,0) = y(i)
end do
do k = 0, n-1
do i=0 , n-k-1
c T(i,k+1)=T(i+1,k)-(((x(i+k+1)-xs)/(x(i+k+1)-x(i)))
c * *(T(i+1,k)-T(i,k))
T(i,k+1)=T(i+1,k)-(((x(i+k+1)-xs)/
* (x(i+k+1)-x(i)))*(T(i+1,k)-T(i,k)))
end do
end do
ys = T(0,n)
return
end

C====================================================================
C Subroutine per la valutazione del polinomio di interpolazione
C in un punto con la Schema di Neville-Aitken
C====================================================================
SUBROUTINE nagg (x, y, T, maxn, xs, ys, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore contenente i nodi (0:n)
C y vettore contenente i valori della funzione nei nodi (0:n)
C xs valore di cui effettuare la valutazione
C n il numero dei nodi e' (n+1)
C
C Variabili USCITA:
C ys approssimazione calcolata

C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxn
INTEGER k, i, r, j
REAL*8 x(0:*), y(0:*), T(0:maxn,0:*)
REAL*8 xs, ys ,w1, w2
C====================
C Parte esecutiva
C====================
C passo in colonna 0 il valore della funzione
C calcolato nel nuovo punto aggiunto
T(n,0) = y(n)
print*,T(0,0)
print*,T(1,0)
print*,T(2,0)
print*,T(3,0)
print*,T(4,0)
print*,T(5,0)
do i = 1, n
print*, i
w1=(x(n-i+1)-xs)/(x(n-i+1)-x(n-i))
w2=T(n-i+1,i-1)-T(n-i,i-1)

T(n-i,i)=T(n-i+1,i-1)-w1*w2

print*,T(n-i+1,i-1),T(n-i,i-1)
end do
ys = T(0,n)
return
end
C ====================================================================
C Funzione assegnata
C====================================================================
real*8 function f(x)
implicit none
real*8 x
f=dsqrt(x)
return
end