fortran - Fortran でのエルミート補間の何が問題になっていますか?

Question

エルミート補間の問題

与えられた x のセットの関数値と導関数値のニュートン除算差を見つけようとしています。私のコードは小さな例では機能しますが、大きな例では失敗します。はっきりとわかるように、私の答えは元の関数値よりもはるかに大きくなっています。

誰かが私が間違っていることを知っていますか?

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

値

xf(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

出力

1.1599999999999999 62.040113431002474

1.3200000000000001 180.40121445431600

1.4800000000000000 212.36319446149312

1.6399999999999999 228.61845650513027

1.8000000000000000 245.11610836104515