Что не так с моей эрмитовой интерполяцией в Фортране?

Question

Эрмит Интерполяционные беды

Я пытаюсь найти разделительную разницу Ньютона для функции и производных значений данного набора х. Я сталкиваюсь с серьезными проблемами с моим кодом, работающим на крошечных примерах, но не на больших. Как ясно видно, мои ответы намного больше, чем исходные значения функций.

Кто-нибудь знает, что я делаю не так?

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,9999

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

Answer 1

Вы получаете доступ к массиву newtonDivDiff за границами.

Вы сначала выделяете его как 0:n (основная программа n) тогда вы переходите к подпрограмме getNewtonDivDiff как 0:n (подпрограмма n) но вы проходите m (m=n*2+1) к аргументу n, Это означает, что вы говорите подпрограмме, что массив имеет границы 0:m который 0:9, но у него есть только границы 0:4,

Отладить программу в ее нынешнем виде довольно сложно, пришлось использовать valgrind, Если вы переместите свои подпрограммы в модуль и измените фиктивные аргументы на массивы предполагаемой формы (:,:) потом перепроверяю проверку в гфортране (-fcheck=all) поймает ошибку.

Другие заметки:

kind=8 уродливо, 8 может означать разные вещи для разных компиляторов. Если вы хотите 64-битные переменные, вы можете использовать kind=real64 (real64 приходит из модуля iso_fortran_env в Fortran 2008) или используйте selected_real_kind() ( Параметр типа Fortran 90)

Вам не нужно освобождать ваши локальные массивы в подпрограммах, они освобождаются автоматически.

Ваш deallocate оператор в основной программе находится после оператора стопа, он никогда не будет выполнен. Я бы просто удалил stopНет причин иметь это.