代码之家 › 专栏 › 技术社区 › gansub

Newton Raphson迭代-无法迭代

newtons-method floating-accuracy precision fortran

gansub · 技术社区 · 7 年前

我不确定这个问题是在这里还是在其他地方(或者根本不是在任何地方)。我继承了Fortran 90代码,该代码执行牛顿-拉斐逊插值,其中温度的对数与压力的对数进行插值。

插值就是这种类型

t =  a ln(p) + b

其中,a,b定义为

a = ln(tup/tdwn)/(alogpu - alogpd)

和

 b = ln T - a * ln P

这是测试程序。它仅针对单个迭代显示。但实际程序在k、j和i上运行三个FOR循环。实际上,pthta是一个3D数组(k、j、i),thta是一个1D数组(k)

 program test

 implicit none
 integer,parameter :: dp = SELECTED_REAL_KIND(12,307)
 real(kind=dp)  kappa,interc,pres,dltdlp,tup,tdwn
 real(kind=dp)  pthta,alogp,alogpd,alogpu,thta,f,dfdp,p1
 real(kind=dp) t1,resid,potdwn,potup,pdwn,pup,epsln,thta1
 integer i,j,kout,n,maxit,nmax,resmax


 kappa = 2./7.
 epsln = 1. 
 potdwn = 259.39996337890625
 potup =  268.41687198359159
 pdwn = 100000.00000000000
 pup =  92500.000000000000
 alogpu =  11.43496392350051
 alogpd =  11.512925464970229
 thta = 260.00000000000000
 alogp = 11.512925464970229
 ! known temperature at lower level
 tdwn = potdwn * (pdwn/100000.)** kappa
 ! known temperature at upper level 
 tup = potup *(pup/100000.)** kappa
 ! linear change of temperature wrt lnP between different levels
 dltdlp = dlog(tup/tdwn)/(alogpu-alogpd)
 ! ln(T) value(intercept) where Pressure is 1 Pa and assume a linear
 ! relationship between P and T
 interc = dlog(tup) - dltdlp*alogpu
 ! Initial guess value for pressure

 pthta = exp((dlog(thta)-interc-kappa*alogp)/(dltdlp-kappa))
 n=0

 1900 continue
 !First guess of temperature at intermediate level
 t1 = exp(dltdlp * dlog(pthta)+interc)
 !Residual error when calculating Newton Raphson iteration(Pascal)
 resid = pthta - 100000.*(t1/thta)**(1./kappa)
 print *, dltdlp,interc,t1,resid,pthta

 if (abs(resid) .gt. epsln) then
   n=n+1
 if (n .le. nmax) then
    ! First guess of potential temperature given T1 and
    ! pressure level guess
    thta1 = t1 * (100000./pthta)**kappa
    f= thta - thta1
    dfdp = (kappa-dltdlp)*(100000./pthta)**kappa*exp(interc + (dltdlp -1.)*dlog(pthta))
    p1 = pthta - f/dfdp
    if (p1 .le. pdwn) then
       if (p1 .ge. pup) then
          pthta = p1
          goto 1900
       else
          n = nmax
       end if
    end if
 else
    if (resid .gt. resmax) resmax = resid
    maxit = maxit+1
    goto 2100
 end if
end if

2100 continue

end program test

使用数据文件中的真实数据运行此程序时,resid的值如下所示

2.7648638933897018E-010

而且在整个执行过程中没有太大区别。大多数值都在范围内

1E-10 and 1E-12

因此,给定这些值,以下IF条件

IF (abs(resid) .gt. epsln)

永远不会被调用,牛顿-拉斐逊迭代也永远不会执行。所以我研究了两种方法来实现这一点。一是在这两个步骤中删除指数调用

pthta = exp((dlog(thta)-interc-kappa*alogp)/(dltdlp-kappa))

t1 = exp(dltdlp * dlog(pthta)+interc)

i、 e.将所有内容保持在对数空间中,并在牛顿-拉斐逊迭代完成后取指数。这一部分确实没有问题。

另一种方法是截断

t1 = exp(dltdlp * dlog(pthta)+interc)

当我将其截断为整数时,resid的值从 1E-10至813。我不理解截断该函数调用如何导致如此大的值更改。截断该结果确实会导致成功完成。因此,我不确定哪一种方式更适合继续下去。

我如何才能决定哪种方法更适合实现这一点?

1 回复 | 直到 7 年前

Scientist 7 年前

从研究的角度来看,我认为您的第一个解决方案可能是更合适的方法。在物理模拟中,应始终使用根据定义始终为正的特性的对数。在上述代码中,这些是温度和压力。无论您使用Fortran或任何其他编程语言,或任何可能的变量类型,严格正定的物理变量通常会导致计算中的上溢和下溢。如果有什么事情可以发生,它就会发生。

其他物理量也是如此,例如,能量(伽马射线暴的典型能量约为10^54 ergs)、任意维物体的体积(半径为10米的100维球体的体积约为10^100),甚至概率(许多统计问题中的似然函数可以取约10^1000}或更小的值)。使用正定变量的对数变换将使您的代码能够处理大约10^10^307的数字(对于双精度变量)。

关于代码中使用的Fortran语法,还有一些注意事项:

变量 RESMAX 在代码中使用,无需初始化。
当给变量赋值时,适当地指定文字常量的类型很重要,否则,程序结果可能会受到影响。例如,以下是在调试模式下使用英特尔Fortran编译器2018编译的原始代码的输出:
```
  -0.152581477302743        7.31503025786548        259.608693509165
  -3.152934473473579E-002   99474.1999921620
```
这是相同代码的输出,但所有文字常量都以kind参数作为后缀 _dp (请参阅以下代码的修订版本):
```
  -0.152580456940175        7.31501855886952        259.608692604963
  -8.731149137020111E-011   99474.2302854451
```
此答案中修订代码的输出与上述问题中原始代码的输出略有不同。
没有必要使用 .gt. , .ge. , .le. , .lt. , ..., 用于比较。据我所知,这些是传统的FORTRAN语法。使用更有吸引力的符号( < , > , <= , >= , == )用于比较。
没有必要使用 GOTO Fortran程序中的语句。这也是遗留的FORTRAN。通常,简单优雅的do循环和if块可以替换 转到 声明,就像下面修订的代码一样。
不再需要在Fortran中使用特定种类的内部函数(例如 dexp , dlog , ... 用于双精度)。几乎所有(或许所有)Fortran内部函数都有泛型名称( exp , log , ...) 在当前Fortran标准中。

以下是该问题中程序的修订版,解决了上述所有过时语法,以及处理非常大或非常小的正定变量的问题(我可能在日志转换一些永远不会导致溢出或下溢的变量时走得太远了,但我在这里的目的只是展示正定变量的日志转换背后的逻辑,以及如何处理它们的算术,而不会潜在地导致溢出/下溢/错误)。

program test

implicit none
integer,parameter :: dp = SELECTED_REAL_KIND(12,307)
real(kind=dp)  kappa,interc,pres,dltdlp,tup,tdwn
real(kind=dp)  pthta,alogp,alogpd,alogpu,thta,f,dfdp,p1
real(kind=dp) t1,resid,potdwn,potup,pdwn,pup,epsln,thta1
integer i,j,kout,n,maxit,nmax,resmax

real(kind=dp) :: log_resmax, log_pthta, log_t1, log_dummy, log_residAbsolute, sign_of_f
real(kind=dp) :: log_epsln, log_pdwn, log_pup, log_thta, log_thta1, log_p1, log_dfdp, log_f
logical :: residIsPositive, resmaxIsPositive, residIsBigger

log_resmax = log(log_resmax)
resmaxIsPositive = .true.

kappa = 2._dp/7._dp
epsln = 1._dp 
potdwn = 259.39996337890625_dp
potup =  268.41687198359159_dp
pdwn = 100000.00000000000_dp
pup =  92500.000000000000_dp
alogpu =  11.43496392350051_dp
alogpd =  11.512925464970229_dp
thta = 260.00000000000000_dp
alogp = 11.512925464970229_dp

log_epsln = log(epsln) 
log_pup =  log(pup)
log_pdwn = log(pdwn)
log_thta = log(thta)

! known temperature at lower level
tdwn = potdwn * (pdwn/1.e5_dp)**kappa
! known temperature at upper level 
tup = potup *(pup/1.e5_dp)** kappa
! linear change of temperature wrt lnP between different levels
dltdlp = log(tup/tdwn)/(alogpu-alogpd)
! ln(T) value(intercept) where Pressure is 1 Pa and assume a linear
! relationship between P and T
interc = log(tup) - dltdlp*alogpu
! Initial guess value for pressure

!pthta = exp( (log(thta)-interc-kappa*alogp) / (dltdlp-kappa) )
log_pthta = ( log_thta - interc - kappa*alogp ) / ( dltdlp - kappa )

n=0

MyDoLoop: do

    !First guess of temperature at intermediate level
    !t1 = exp(dltdlp * log(pthta)+interc)
    log_t1 = dltdlp * log_pthta + interc

    !Residual error when calculating Newton Raphson iteration(Pascal)
    !resid = pthta - 1.e5_dp*(t1/thta)**(1._dp/kappa)
    log_dummy = log(1.e5_dp) + ( log_t1 - log_thta ) / kappa
    if (log_pthta>=log_dummy) then
      residIsPositive = .true.
      log_residAbsolute = log_pthta + log( 1._dp - exp(log_dummy-log_pthta) )
    else
      residIsPositive = .false.
      log_residAbsolute = log_dummy + log( 1._dp - exp(log_pthta-log_dummy) )
    end if

    print *, "log-transformed values:"
    print *, dltdlp,interc,log_t1,log_residAbsolute,log_pthta
    print *, "non-log-transformed values:"
    if (residIsPositive) print *, dltdlp,interc,exp(log_t1),exp(log_residAbsolute),exp(log_pthta)
    if (.not.residIsPositive) print *, dltdlp,interc,exp(log_t1),-exp(log_residAbsolute),exp(log_pthta)

    !if (abs(resid) > epsln) then
    if ( log_residAbsolute > log_epsln ) then
        n=n+1
        if (n <= nmax) then
            ! First guess of potential temperature given T1 and
            ! pressure level guess
            !thta1 = t1 * (1.e5_dp/pthta)**kappa
            log_thta1 = log_t1 + ( log(1.e5_dp)-log_pthta ) * kappa
            !f = thta - thta1
            if ( log_thta>=thta1 ) then
              log_f = log_thta + log( 1._dp - exp( log_thta1 - log_thta ) )
              sign_of_f = 1._dp
            else
              log_f = log_thta + log( 1._dp - exp( log_thta - log_thta1 ) )
              sign_of_f = 1._dp
            end if
            !dfdp = (kappa-dltdlp)*(1.e5_dp/pthta)**kappa*exp(interc + (dltdlp -1._dp)*log(pthta))
            ! assuming kappa-dltdlp>0 is TRUE always:
            log_dfdp = log(kappa-dltdlp) + kappa*(log(1.e5_dp)-log_pthta) + interc + (dltdlp -1._dp)*log_pthta
            !p1 = pthta - f/dfdp
            ! p1 should be, by definition, positive. Therefore:
            log_dummy = log_f - log_dfdp
            if (log_pthta>=log_dummy) then
                log_p1 = log_pthta + log( 1._dp - sign_of_f*exp(log_dummy-log_pthta) )
            else
                log_p1 = log_dummy + log( 1._dp - sign_of_f*exp(log_pthta-log_dummy) )
            end if
            !if (p1 <= pdwn) then
            if (log_p1 <= log_pdwn) then
               !if (p1 >= pup) then
               if (log_p1 >= log_pup) then
                  log_pthta = log_p1
                  cycle MyDoLoop
               else
                  n = nmax
               end if
            end if
        else
            !if (resid > resmax) resmax = resid
            residIsBigger = ( residIsPositive .and. resmaxIsPositive .and. log_residAbsolute>log_resmax ) .or. &
                            ( .not.residIsPositive .and. .not.resmaxIsPositive .and. log_residAbsolute<log_resmax ) .or. &
                            ( residIsPositive .and. .not. resmaxIsPositive )
            if ( residIsBigger ) then
                log_resmax = log_residAbsolute
                resmaxIsPositive = residIsPositive
            end if
            maxit = maxit+1
        end if
    end if

    exit MyDoLoop

end do MyDoLoop

end program test

以下是该程序的示例输出,与原始代码的输出非常一致:

log-transformed values:
 -0.152580456940175        7.31501855886952        5.55917546888014
  -22.4565579499410        11.5076538974964
 non-log-transformed values:
 -0.152580456940175        7.31501855886952        259.608692604963
 -1.767017293116268E-010   99474.2302854451

为了进行比较,以下是原始代码的输出:

 -0.152580456940175        7.31501855886952        259.608692604963
 -8.731149137020111E-011   99474.2302854451