1

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

价值观

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

4

1 回答 1

4

您正在newtonDivDiff越界访问数组。

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

就目前的情况调试程序非常困难，我不得不使用valgrind. 如果您将子例程移至模块并将虚拟参数更改为假定的形状数组(:,:)，则 gfortran ( ) 中的边界检查-fcheck=all将捕获错误。

其他注意事项：

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

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

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

于 2015-10-17T00:11:16.707 回答