Whats wrong with my Hermite Interpolation in Fortran?

Question

Hermite Interpolation woes

I am trying to find the Newton Dividing Differences for the function and derivative values of a given set of x's. I'm running into serious problems with my code working for tiny examples, but failing on bigger one's. As is clearly visible, my answers are very much larger than they original function values.

Does anybody have any idea what I'm doing wrong?

program inter
  implicit none
  integer ::n,m
  integer ::i
  real(kind=8),allocatable ::xVals(:),fxVals(:),newtonDivDiff(:),dxVals(:),zxVals(:),zdxVals(:),zfxVals(:)
  real(kind=8) ::Px
  real(kind=8) ::x


  Open(Unit=8,File="data/xVals")
  Open(Unit=9,File="data/fxVals")
  Open(Unit=10,File="data/dxVals")

  n = 4 ! literal number of data pts
  m = n*2+1

  !after we get the data points allocate the space
  allocate(xVals(0:n))
  allocate(fxVals(0:n))
  allocate(dxVals(0:n))
  allocate(newtonDivDiff(0:n))

  !allocate the zvalue arrays
  allocate(zxVals(0:m))
  allocate(zdxVals(0:m))
  allocate(zfxVals(0:m))

  !since the size is the same we can read in one loop
  do i=0,n
    Read(8,*) xVals(i)
    Read(9,*) fxVals(i)
    Read(10,*) dxVals(i)
  end do  
 ! contstruct the z illusion 
  do i=0,m,2
    zxVals(i) = xVals(i/2)
    zxVals(i+1) = xVals(i/2)

    zdxVals(i) = dxVals(i/2)
    zdxVals(i+1) = dxVals(i/2)

    zfxVals(i) = fxVals(i/2)
    zfxVals(i+1) = fxVals(i/2)

  end do
  !slightly modified business as usual
  call getNewtonDivDiff(zxVals,zdxVals,zfxVals,newtonDivDiff,m)
  do i=0,n 
    call evaluatePolynomial(m,newtonDivDiff,xVals(i),Px,zxVals)
    print*, xVals(i) ,Px
  end do
  close(8)
  close(9)
  close(10)
  stop

  deallocate(xVals,fxVals,dxVals,newtonDivDiff,zxVals,zdxVals,zfxVals)
end program inter

subroutine getNewtonDivDiff(xVals,dxVals,fxVals,newtonDivDiff,n)
    implicit none
    integer ::i,k
    integer, intent(in) ::n
    real(kind=8), allocatable,dimension(:,:) ::table
    real(kind=8),intent(in) ::xVals(0:n),dxVals(0:n),fxVals(0:n)
    real(kind=8), intent(inout) ::newtonDivDiff(0:n)
    allocate(table(0:n,0:n))

    table = 0.0d0

    do i=0,n
        table(i,0) = fxVals(i)
    end do

    do k=1,n
        do i = k,n
            if( k .eq. 1 .and. mod(i,2) .eq. 1) then
                table(i,k) = dxVals(i)
            else
                table(i,k) = (table(i,k-1) - table(i-1,k-1))/(xVals(i) - xVals(i-k))
            end if
        end do 
    end do

    do i=0,n
        newtonDivDiff(i) = table(i,i)
        !print*, newtonDivDiff(i)
    end do
    deallocate(table)
end subroutine getNewtonDivDiff

subroutine evaluatePolynomial(n,newtonDivDiff,x,Px,xVals)
    implicit none
    integer,intent(in) ::n
    real(kind=8),intent(in) ::newtonDivDiff(0:n),xVals(0:n)
    real(kind=8),intent(in) ::x
    real(kind=8), intent(out) ::Px
    integer ::i
    Px = newtonDivDiff(n)
    do i=n,1,-1
        Px  = Px * (x-  xVals(i-1)) + newtonDivDiff(i-1)
    end do
end subroutine evaluatePolynomial

Values

x f(x) f'(x)

1.16, 1.2337, 2.6643

1.32, 1.6879, 2.9989

1.48, 2.1814, 3.1464

1.64, 2.6832, 3.0862

1.8, 3.1553, 2.7697

Output

1.1599999999999999 62.040113431002474

1.3200000000000001 180.40121445431600

1.4800000000000000 212.36319446149312

1.6399999999999999 228.61845650513027

1.8000000000000000 245.11610836104515

Answer 1

You are accessing array newtonDivDiff out of bounds.

You are first allocating it as 0:n (main program's n) then you are passing to subroutine getNewtonDivDiff as 0:n (the subroutine's n) but you pass m (m=n*2+1) to the argument n. That means you tell the subroutine that the array has bounds 0:m which is 0:9, but it has only bounds 0:4.

It is quite difficult to debug the program as it stands, I had to use valgrind. If you move your subroutines to a module and change the dummy arguments to assumed shape arrays (:,:) then the bound checking in gfortran (-fcheck=all) will catch the error.

Other notes:

kind=8 is ugly, 8 can mean different things for different compilers. If you want 64bit variables, you can use kind=real64 (real64 comes from module iso_fortran_env in Fortran 2008) or use selected_real_kind() (Fortran 90 kind parameter)

You do not have to deallocate your local arrays in the subroutines, they are deallocated automatically.

Your deallocate statement in the main program is after the stop statement, it will never be executed. I would just delete the stop, there is no reason to have it.