Fortran和Python中的准确性和多项式求值

我将在此处和下面的上下文中进行快速描述.我正在Python和Fortran中评估双变量多项式(2个变量的多项式),并得到不同的结果.我的测试用例的相对误差-4.23e-3-足够大,不会由于精度差异而明显.以下代码片段使用相当原始的类型和相同的算法来尝试使事物尽可能具有可比性.关于差异有什么线索吗？我尝试过改变精度(在Fortran中使用selected_real_kind,在Python中使用numpy.float128),对Fortran进行编译(特别是优化级别),以及算法(霍纳方法,numpy评估).关于差异有什么线索吗？任一版本的代码中都有错误？我已经看过Precision discrepancy between Fortran and Python (sin function),但是还没有机会用不同的编译器完全测试它.

Python：

#!/usr/bin/env python
""" polytest.py
Test calculation of a bivariate polynomial.
"""

# Define polynomial coefficients
coeffs = (
    (101.34274313967416, 100015.695367145, -2544.5765420363),
    (5.9057834791235253,-270.983805184062,1455.0364540468),
    (-12357.785933039,1455.0364540468,-756.558385769359),
    (736.741204151612,-672.50778314507,499.360390819152)
)
nx = len(coeffs)
ny = len(coeffs[0])

# Values of variables
x0 = 0.0002500000000011937
y0 = -0.0010071334522899211

# Calculate polynomial by looping over powers of x, y
z = 0.
xj = 1.
for j in range(nx):
    yk = 1.
    for k in range(ny):
        curr = coeffs[j][k] * xj * yk
        z += curr
        yk *= y0
    xj *= x0

print(z)   # 0.611782174444

Fortran：

! polytest.F90
! Test calculation of a bivariate polynomial.

program main

  implicit none
  integer, parameter :: dbl = kind(1.d0)
  integer, parameter :: nx = 3, ny = 2
  real(dbl), parameter :: x0 = 0.0002500000000011937, &
                          y0 = -0.0010071334522899211
  real(dbl), dimension(0:nx,0:ny) :: coeffs
  real(dbl) :: z, xj, yk, curr
  integer :: j, k

  ! Define polynomial coefficients
  coeffs(0,0) = 101.34274313967416d0
  coeffs(0,1) = 100015.695367145d0
  coeffs(0,2) = -2544.5765420363d0
  coeffs(1,0) = 5.9057834791235253d0
  coeffs(1,1) = -270.983805184062d0
  coeffs(1,2) = 1455.0364540468d0
  coeffs(2,0) = -12357.785933039d0
  coeffs(2,1) = 1455.0364540468d0
  coeffs(2,2) = -756.558385769359d0
  coeffs(3,0) = 736.741204151612d0
  coeffs(3,1) = -672.50778314507d0
  coeffs(3,2) = 499.360390819152d0

  ! Calculate polynomial by looping over powers of x, y
  z = 0d0
  xj = 1d0
  do j = 0, nx-1
    yk = 1d0
    do k = 0, ny-1
      curr = coeffs(j,k) * xj * yk
      z = z + curr
      yk = yk * y0
    enddo
    xj = xj * x0
  enddo

  ! Print result
  WRITE(*,*) z   ! 0.61436839888538231

end program

编译：gfortran -O0 -o polytest.o polytest.F90

上下文：我正在编写现有Fortran库的纯Python实现,主要是作为练习,但同时还增加了灵活性.我正在以Fortran为基准对我的结果进行基准测试,并且能够在1e-10左右的范围内获得几乎所有结果,但是我无法掌握这一点.其他函数也更加复杂,这使得简单多项式的分歧令人困惑.

特定的系数和测试变量来自该库.实际的多项式实际上在(x,y)中具有度(7,6),因此这里有许多其他系数未包括在内.该算法直接取自Fortran,因此,如果有误,我应该联系原始开发人员.通用函数还可以计算导数,这就是为什么此实现可能不是最佳的原因的一部分-我知道我仍然应该只编写Horner的方法版本,但这并没有改变差异.我只是在以y的较大值计算导数时才注意到这些错误,但该错误的确存在于此更简单的设置中.

解决方法:

为了使结果与Python版本和Fortran版本相匹配,应更正Fortran代码中的两件事.

1.完成后,将特定的双精度类型声明为：

integer, parameter :: dbl = kind(0.d0)

然后,应通过将种类指示符附加为以下内容来定义变量：

real(dbl) :: z
z = 1.0_dbl

例如,在fortran90.org gotchas上对此进行了讨论.语法可能不方便,但是,我没有制定规则.

2. Fortran Do循环迭代由nx和ny控制.您打算访问coeffs的每个元素,但是索引使迭代变得更短了.将nx-1和ny-1分别更改为nx和ny.更好的是,使用Fortran固有的ubound来确定沿所需维度的范围,例如：

do j = 0, ubound(coeffs, dim=1)

下面显示的更新代码可纠正这些问题并打印结果,该结果与python代码产生的结果匹配.

program main
    implicit none
    integer, parameter :: dbl = kind(1.d0)
    integer, parameter :: nx = 3, ny = 2
    real(dbl), parameter :: x0 = 0.0002500000000011937_dbl, &
                          y0 = -0.0010071334522899211_dbl
    real(dbl), dimension(0:nx,0:ny) :: coeffs
    real(dbl) :: z, xj, yk, curr
    integer :: j, k

    ! Define polynomial coefficients
    coeffs(0,0) = 101.34274313967416_dbl
    coeffs(0,1) = 100015.695367145_dbl
    coeffs(0,2) = -2544.5765420363_dbl
    coeffs(1,0) = 5.9057834791235253_dbl
    coeffs(1,1) = -270.983805184062_dbl
    coeffs(1,2) = 1455.0364540468_dbl
    coeffs(2,0) = -12357.785933039_dbl
    coeffs(2,1) = 1455.0364540468_dbl
    coeffs(2,2) = -756.558385769359_dbl
    coeffs(3,0) = 736.741204151612_dbl
    coeffs(3,1) = -672.50778314507_dbl
    coeffs(3,2) = 499.360390819152_dbl

    ! Calculate polynomial by looping over powers of x, y
    z = 0.0_dbl
    xj = 1.0_dbl
    do j = 0, ubound(coeffs, dim=1)
        yk = 1.0_dbl
        do k = 0, ubound(coeffs, dim=2)
            print "(a,i0,a,i0,a)", "COEFF(",j,",",k,")="
            print *, coeffs(j,k)
            curr = coeffs(j,k) * xj * yk
            z = z + curr
            yk = yk * y0
        enddo
        xj = xj * x0
    enddo

    ! Print result
    WRITE(*,*) z   ! Result: 0.611782174443735
end program  

标签：precision,floating-accuracy,fortran,polynomial-math,python
来源： https://codeday.me/bug/20191111/2017404.html