! ! Spline for *equally* spaced data ! Derivation in Numerical Methods for Scientists and Engineers by R.W.Hamming 1973 IBSN:0486652416 ! See LaTeX documentation below ! You must provide a matrix inversion routine: MatrixInverseFunction(A) ! Charles O'Neill 16 April 2008 ! www.caselab.okstate.edu ! ! Copyright (c) 2008 Charles O'Neill ! ! Permission is hereby granted, free of charge, to any person ! obtaining a copy of this software and associated documentation ! files (the "Software"), to deal in the Software without ! restriction, including without limitation the rights to use, ! copy, modify, merge, publish, distribute, sublicense, and/or sell ! copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following ! conditions: ! ! The above copyright notice and this permission notice shall be ! included in all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES ! OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND ! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT ! HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, ! WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR ! OTHER DEALINGS IN THE SOFTWARE. module SplineRoutines implicit none public :: PlanSpline, ExecuteSpline, splinetest, WriteSpline integer, parameter :: DP = 8 ! Double Precision type SplineADT private integer :: N real(DP) :: h real(DP),allocatable :: xdata(:), fdata(:) real(DP),allocatable :: Yprimeprime(:) end type interface PlanSpline module procedure plansplineinterpolate end interface interface DestroySpline module procedure destroysplineinterpolate end interface interface WriteSpline module procedure writesplineinterpolate end interface interface ExecuteSpline module procedure executesplineinterpolate end interface contains subroutine splinetest() real(DP) :: xdata(4) real(DP) :: fdata(4) real(DP) :: EndPointCurvature(2) real(DP) :: x,f type (SplineADT) :: Spl xdata = [0, 1, 2, 3] fdata = [1, 2, 5, 10] EndPointCurvature = [2.0d0, 2.0d0] x = 1.5 Spl = PlanSpline(xdata,fdata) f = ExecuteSpline(Spl,x) call DestroySpline(Spl) Spl = PlanSpline(xdata,fdata,EndPointCurvature) f = ExecuteSpline(Spl,x) call DestroySpline(Spl) if(f==3.25) write(*,*) "Spline Passed" end subroutine subroutine destroysplineinterpolate(This) type (SplineADT) :: This deallocate(This%xdata) deallocate(This%fdata) deallocate(This%Yprimeprime) end subroutine subroutine writesplineinterpolate(This) type (SplineADT) :: This integer :: i write(*,*) "X data" do i=1,This%N write(*,*) This%xdata(i) enddo write(*,*) "Function data" do i=1,This%N write(*,*) This%fdata(i) enddo write(*,*) "Yprimeprime" do i=1,This%N write(*,*) This%Yprimeprime(i) enddo end subroutine function plansplineinterpolate(xdata,fdata,EndPointCurvatureInput) use LU implicit none type (SplineADT) :: plansplineinterpolate real(DP) :: xdata(:), fdata(:) real(DP),optional :: EndPointCurvatureInput(2) real(DP) :: EndPointCurvature(2) real(DP),allocatable :: A(:,:), B(:), D(:) integer :: i,j, N real(DP) :: h ! EndPoint Curvature Information; Otherwise assume Natural Spline if(present(EndPointCurvatureInput))then EndPointCurvature=EndPointCurvatureInput else EndPointCurvature=[0.0d0, 0.0d0] endif ! Sizing N=size(xdata) plansplineinterpolate%N = N allocate(A(N,N)) allocate(B(N)) allocate(D(N)) allocate(plansplineinterpolate%xdata(N)) allocate(plansplineinterpolate%fdata(N)) allocate(plansplineinterpolate%Yprimeprime(N)) ! Coupling Coefficients (Left Hand Side) A=0.0d0 do i=2,N-1 do j=2,N-1 if(i==j)then A(i,j) = 4.0d0 ! Diagonal Term elseif(i==(j-1))then A(i,j) = 1.0d0 ! Just Off Diagonal Term elseif(i==(j+1))then A(i,j) = 1.0d0 ! Just Off Diagonal Term else A(i,j) = 0.0d0 ! Totally Off Diagonal Term endif enddo enddo ! Spacing h = xdata(2)-xdata(1) plansplineinterpolate%h = h ! 2nd Derivative do i=2,N-1 D(i) = (fdata(i+1)-2.d0*fdata(i)+fdata(i-1))/(2.0d0*h**2) enddo ! Numerical Derivatives at Points (Right Hand Side) B(i) = 0.0d0 do i=2,N-1 B(i) = 12.0d0*D(i) enddo ! Nullify nonused endpoints (via identity matrix and zero vector) A(1,1) = 1.0d0 A(N,N) = 1.0d0 B(1) = 0.0d0 B(N) = 0.0d0 ! Correct A and B for known constraints B(2) = B(2) -EndPointCurvature(1)/h B(N-1) = B(N-1) -EndPointCurvature(2)/h ! Store Interpolation Points plansplineinterpolate%xdata = xdata plansplineinterpolate%fdata = fdata ! Solve for Spline plansplineinterpolate%Yprimeprime = matmul(MatrixInverseFunction(A),B) end function function executesplineinterpolate(This,x) real(DP) :: executesplineinterpolate real(DP) :: x type (SplineADT) :: This integer :: k, i ! x is between k and k+1 locations in xdata do i=1,(This%N-1) if( (xThis%xdata(i)) )then k=i exit endif enddo executesplineinterpolate = This%fdata(k)*(This%xdata(k+1)-x)/This%h + This%fdata(k+1)*(x-This%xdata(k))/This%h & - This%h**2/6.0d0 * This%Yprimeprime(k) * ( (This%xdata(k+1)-x)/This%h - ((This%xdata(k+1)-x)/This%h)**3 ) & - This%h**2/6.0d0 * This%Yprimeprime(k+1) * ( (x-This%xdata(k))/This%h - ((x-This%xdata(k))/This%h )**3 ) end function end module !--------------------------------------------- ! ! LaTeX Documentation ! !This spline is a polynomial of order 3 between N points.\[ !y(x)=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i}\] !The spacing is\[ !h_{i}=x_{i+1}-x_{i}\] !At each point, the following should be satisfied\[ !h_{i-1}y_{i-1}^{'}+2\left(h_{i-1}+h_{i}\right)y_{i}^{''}+h_{i}y_{i+1}^{''}=6\left(\frac{y_{i+1}-y_{i}}{h_{i}}-\frac{y_{i}-y_{i-1}}{h_{i-1}}\right)\] !Altogether this requires solving the following\[ !\left[\begin{array}{cccc} !2\left(h_{1}+h_{2}\right) & h_{2}\\ !h_{2} & 2\left(h_{2}+h_{3}\right) & \ddots\\ ! & \ddots & \ddots & h_{n-2}\\ ! & & h_{n-2} & 2\left(h_{n-2}+h_{n-1}\right)\end{array}\right]\left[\begin{array}{c} !S_{2}\\ !S_{3}\\ !\vdots\\ !S_{n-1}\end{array}\right]=\left[\begin{array}{c} !\frac{y_{3}-y_{2}}{h_{2}}-\frac{y_{2}-y_{1}}{h_{1}}\\ !\frac{y_{4}-y_{3}}{h_{3}}-\frac{y_{3}-y_{2}}{h_{2}}\\ !\vdots\\ !\frac{y_{n}-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\end{array}\right]\] !The values of $a_{i}$,$b_{i}$,$c_{i}$,$d_{i}$ are\begin{eqnarray*} !a_{i} & = & \frac{\left(S_{i+1}-S_{i}\right)}{6h_{i}}\\ !b_{i} & = & \frac{S_{i}}{2}\\ !c_{i} & = & \frac{y_{i+1}-y_{i}}{h_{i}}-\frac{2h_{i}S_{i}+h_{i}S_{i+1}}{6}\\ !d_{i} & = & y_{i}\end{eqnarray*} !For equally spaced points, Hamming uses\[ !\left[\begin{array}{cccc} !4 & 1\\ !1 & 4 & 1\\ ! & \ddots & \ddots & 1\\ ! & & 1 & 4\end{array}\right]\left[\begin{array}{c} !y_{2}^{''}\\ !y_{3}^{''}\\ !\vdots\\ !y_{n-1}^{''}\end{array}\right]=\left[\begin{array}{c} !12d_{2}-\frac{y_{1}^{''}}{h}\\ !12d_{3}\\ !\vdots\\ !12d_{n-1}-\frac{y_{n}^{''}}{h}\end{array}\right]\] !Hamming's $d_{i}$ terms are\[ !d_{i}=\frac{\left(y_{i+1}-2y_{i}+y_{i-1}\right)}{2h^{2}}\] !This gives an interpolation function\begin{eqnarray*} !f_{i}(x) & = & y_{i}\left(\frac{x_{i+1}-x}{h_{i}}\right)+y_{i+1}\left(\frac{x-x_{i}}{h_{i}}\right)\\ ! & & -\frac{h_{i}^{2}}{6}y_{i}^{''}\left(\frac{x_{i+1}-x}{h_{i}}-\left(\frac{x_{i+1}-x}{h_{i}}\right)^{3}\right)\\ ! & & -\frac{h_{i}^{2}}{6}y_{i+1}^{''}\left(\frac{x-x_{i}}{h_{i}}-\left(\frac{x-x_{i}}{h_{i}}\right)^{3}\right)\end{eqnarray*} !Now, we can interpolate betwen the given nodes and points. ! !