Użytkownik "Michoo" <m...@v...pl> napisał w wiadomości grup
dyskusyjnych:h2tcvs$kbn$...@n...onet.pl...
> mieć też równe pochodne do n-1 rzędu - jeżeli masz tylko tablicę z
> wartościami funkcji bazowej to jak zrobiłeś splajny stopnia większego niż
> 1(gdzie wychodzi łamana - jak we wzorze trapezów)?

Interpolacja spline'em. Ma się zgadzać wartość i m-1 pochodnych, ma być
wielomian stopnia m-tego pomiędzy węzłami. Z tego wychodzi jaki to spline.
Ze spline nietrudno policzyć całkę - skoro to wielomian, to całka jest
wielomianem, potem tylko granice całkowania - sumowanie po wszystkich
kawałkach. Patrz Forman S. Acton, Numerical Methods that Work, Princeton
University, (ISBN-13: 9780883854501 | ISBN-10: 0883854503) - niezła książka
na swój sposób. Nota bene, odjazd to jest jak to się robi na zespolonych
(zmiana deklaracji) - też działa, bo czemu miałoby nie działać?! :)

slawek


!***************************************************
****************************************************
****************************
!
! Spline integration
!
!***************************************************
****************************************************
****************************

double precision function intspl(y,n,h)

integer, parameter :: nmax = 16384

double precision h
double precision :: y(1)
double precision :: sa,sb,sc,sd,un
double precision, dimension(nmax) :: s,u

s(1) = 0.
s(n) = 0.
u(1) = 0.
sa = 0.
sb = 0.
sc = 0.
sd = 0.

if(n .gt. 1) then
do i=2,n-1
s(i)= -1.0 / (s(i-1) + 4.0)
u(i)=(y(i+1)-2*y(i)+y(i-1)) / h
u(i)=(12.0*u(i)/h - u(i-1))/(s(i-1)+4.0)
end do
do i=n-1,1,-1
s(i)=s(i)*s(i+1)+u(i)
end do
do i=1,n
sa = sa + (s(i+1)-s(i))
sb = sb + (s(i)/2.)
sc = sc + ((y(i+1)-y(i))/h - (2.*h*s(i)+h*s(i+1))/6.)
sd = sd + y(i)
end do
end if

intspl = h**4/4.*sa + h**3/3.*sb + h**2/2.*sc + h*sd

end function intspl