Eu escrevi um algoritmo rudimentar em Fortran 95 para calcular o gradiente de uma função (um exemplo do qual é prescrito no código) usando diferenças centrais aumentadas com um procedimento conhecido como extrapolação de Richardson.

function f(n,x)
! The scalar multivariable function to be differentiated

integer :: n
real(kind = kind(1d0)) :: x(n), f

f = x(1)**5.d0 + cos(x(2)) + log(x(3)) - sqrt(x(4))

end function f
!=====!
!=====!
!=====!

program gradient
!==============================================================================!
! Calculates the gradient of the scalar function f at x=0using a finite        !
! difference approximation, with a low order Richardson extrapolation.         !
!==============================================================================!

parameter (n = 4, M = 25)
real(kind = kind(1d0)) :: x(n), xhup(n), xhdown(n), d(M), r(M), dfdxi, h0, h, gradf(n)

h0 = 1.d0
x  = 3.d0

! Loop through each component of the vector x and calculate the appropriate
! derivative
do i = 1,n
! Reset step size
h = h0

    ! Carry out M successive central difference approximations of the derivative
do j = 1,M
        xhup = x
        xhdown = x
        xhup(i) = xhup(i) + h
        xhdown(i) = xhdown(i) - h
        d(j) = ( f(n,xhup) - f(n,xhdown) ) / (2.d0*h)
    h = h / 2.d0
    end do

    r = 0.d0
    do k = 3,M      r(k) = ( 64.d0*d(k) - 20.d0*d(k-1) + d(k-2) ) / 45.d0
        if ( abs(r(k) - r(k-1)) < 0.0001d0 ) then
        dfdxi = r(k)
            exit
        end if
    end do

    gradf(i) = dfdxi
end do

! Print out the gradient
write(*,*) " "
write(*,*) " Grad(f(x)) = "
write(*,*) " "
do i = 1,n
    write(*,*) gradf(i)
end do

end program gradient

Em precisão única, ele funciona bem e me dá resultados decentes. Mas quando tento mudar para a precisão dupla, como mostrado no código, recebo um erro ao tentar compilar a afirmação de que a instrução de atribuição

d(j) = ( f(n,xhup) - f(n,xhdown) ) / (2.d0*h)

está produzindo um tipo de incompatibilidadereal(4)/real(8). Eu tentei várias declarações diferentes de precisão dupla, anexadas a cada constante de precisão dupla apropriada no código comd0, e recebo o mesmo erro todas as vezes. Estou um pouco confuso sobre como a funçãof está possivelmente produzindo um único número de precisão.