我的Hermite插值在Fortran中有什么问题?

Question

Hermite插值问题

我试图找到一个给定的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

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

输出

1.159999999999999999 62.040113431002474

1.3200000000000001 180.40121445431600

1.4800000000000000 212.36319446149312

1.639999999999999999 228.61845650513027

1.8000000000000000 245.11610836104515

Answer 1

您正在越界访问数组newtonDivDiff

您首先将其分配为0:n(主程序的n)，然后将其传递给子程序getNewtonDivDiff作为0:n(子程序的n)，但将m (m=n*2+1)传递给参数n。这意味着你告诉子程序数组的边界0:m是0:9，但它只有0:4的边界。

这个程序很难调试，我不得不使用valgrind。如果您将子例程移动到模块中，并将虚拟参数更改为假定的形状数组(:,:)，那么gfortran (-fcheck=all)中的绑定检查将捕获错误。

其他笔记:

kind=8很难看，8对于不同的编译器意味着不同的东西。如果您想要64位变量，您可以使用kind=real64 (real64来自Fortran 2008中的iso_fortran_env模块)或使用selected_real_kind() (Fortran 90类型参数)

您不必在子例程中释放本地数组，它们会自动释放。

主程序中的deallocate语句在stop语句之后，它永远不会被执行。我只会删除stop，没有理由有它。