Newton Raphson 迭代 - 无法迭代
Newton Raphson iteration - unable to iterate
我不确定这个问题是否与这里或其他地方的主题相关（或者根本不与主题相关）。
我继承了执行 Newton Raphson 插值的 Fortran 90 代码，其中温度的对数根据压力的对数进行插值。
插值类型为
t = a ln(p) + b
其中a、b定义为
a = ln(tup/tdwn)/(alogpu - alogpd)
和
b = ln T - a * ln P
这是测试程序。它仅针对单次迭代显示。但实际程序 运行s 超过三个 FOR 循环，遍历 k、j 和 i。实际上 pthta 是一个 3D array(k,j,i) 而 thta 是一个 1D array (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)
永远不会被调用，Newton Raphson 迭代也永远不会被执行。所以我研究了两种方法来让它发挥作用。一是去掉这两步中的exponential call
pthta = exp((dlog(thta)-interc-kappa*alogp)/(dltdlp-kappa))
t1 = exp(dltdlp * dlog(pthta)+interc)
即将所有内容保持在对数 space 中，并在 Newton Raphson 迭代完成后取指数。那部分收敛没有问题。
我尝试完成这项工作的另一种方法是 t运行cate
t1 = exp(dltdlp * dlog(pthta)+interc)
当我运行将其转换为整数时，resid 的值从
1E-10 到 813。我不明白 t运行cating 该函数调用如何导致如此大的值变化。 T运行cating 该结果确实导致成功完成。
所以我不确定哪种方法是进一步进行的更好方法。
我如何确定哪种方法更好？
从研究的角度来看，我认为您的第一个解决方案可能是更合适的方法。在物理模拟中，应该始终使用 by-definition 始终为正的属性的对数。在上面的代码中，这些将是温度和压力。严格 positive-definite 物理变量通常会导致计算上溢和下溢，无论您使用 Fortran 还是任何其他编程语言，或任何可能的变量类型。如果某事可能发生，它就会发生。
其他物理量也是如此，例如，能量（Gamma-Ray-Burst 的典型能量是~10^54 尔格），任意维度物体的体积（100 的体积）半径为 10 米的维球体是 ~ 10^100)，甚至是概率（许多统计问题中的似然函数可以取 ~10^{-1000} 或更小的值）。使用 log-transform 个 positive-definite 变量将使您的代码能够处理大至 ~10^10^307 的数字（对于双精度变量）。
还有一些关于代码中使用的 Fortran 语法的注意事项：
您的代码中使用了变量 RESMAX，但未进行初始化。
给变量赋值时，一定要适当地指定字面常量的种类，否则可能会影响程序结果。例如，这是在调试模式下使用英特尔 Fortran 编译器 2018 编译的原始代码的输出：
-0.152581477302743 7.31503025786548 259.608693509165
-3.152934473473579E-002 99474.1999921620
这里是相同代码的输出，但所有文字常量都以种类参数 _dp 为后缀（请参阅下面代码的修订版本）：
-0.152580456940175 7.31501855886952 259.608692604963
-8.731149137020111E-011 99474.2302854451
此答案中修改后的代码输出与上述问题中原始代码的输出略有不同。
不需要用.gt.、.ge.、.le.、.lt.、……来比较。据我所知，这些是遗留的 FORTRAN 语法。使用更具吸引力的符号（ < 、 > 、 <= 、 >= 、 == ）进行比较。
没有必要在 Fortran 程序中使用 GOTO 语句。这又是遗留的 FORTRAN。通常，简单优雅的 do-loops 和 if-blocks 可以替换 GOTO 语句，就像下面修改后的代码一样。
在 Fortran 中不再需要使用 kind-specific 内部函数（例如 dexp、dlog、...双精度）。在当前的 Fortran 标准中，几乎所有（也许是所有）Fortran 内部函数都有通用名称（exp、log、...）。
以下是对本题程序的修改，解决了上述所有过时的语法，以及处理极大或极小positive-definite变量的问题（我可能走得太远了在 log-transforming 一些永远不会导致溢出或下溢的变量中，但我在这里的目的只是展示 log-transformation 背后的逻辑 positive-definite 变量以及如何处理它们的算法而不可能导致 overflow/underflow/error_in_results).
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
我不确定这个问题是否与这里或其他地方的主题相关（或者根本不与主题相关）。 我继承了执行 Newton Raphson 插值的 Fortran 90 代码，其中温度的对数根据压力的对数进行插值。
插值类型为
t = a ln(p) + b
其中a、b定义为
a = ln(tup/tdwn)/(alogpu - alogpd)
和
b = ln T - a * ln P
这是测试程序。它仅针对单次迭代显示。但实际程序 运行s 超过三个 FOR 循环，遍历 k、j 和 i。实际上 pthta 是一个 3D array(k,j,i) 而 thta 是一个 1D array (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)
永远不会被调用，Newton Raphson 迭代也永远不会被执行。所以我研究了两种方法来让它发挥作用。一是去掉这两步中的exponential call
pthta = exp((dlog(thta)-interc-kappa*alogp)/(dltdlp-kappa))
t1 = exp(dltdlp * dlog(pthta)+interc)
即将所有内容保持在对数 space 中，并在 Newton Raphson 迭代完成后取指数。那部分收敛没有问题。
我尝试完成这项工作的另一种方法是 t运行cate
t1 = exp(dltdlp * dlog(pthta)+interc)
当我运行将其转换为整数时，resid 的值从 1E-10 到 813。我不明白 t运行cating 该函数调用如何导致如此大的值变化。 T运行cating 该结果确实导致成功完成。 所以我不确定哪种方法是进一步进行的更好方法。
我如何确定哪种方法更好？
从研究的角度来看，我认为您的第一个解决方案可能是更合适的方法。在物理模拟中，应该始终使用 by-definition 始终为正的属性的对数。在上面的代码中，这些将是温度和压力。严格 positive-definite 物理变量通常会导致计算上溢和下溢，无论您使用 Fortran 还是任何其他编程语言，或任何可能的变量类型。如果某事可能发生，它就会发生。
其他物理量也是如此，例如，能量（Gamma-Ray-Burst 的典型能量是~10^54 尔格），任意维度物体的体积（100 的体积）半径为 10 米的维球体是 ~ 10^100)，甚至是概率（许多统计问题中的似然函数可以取 ~10^{-1000} 或更小的值）。使用 log-transform 个 positive-definite 变量将使您的代码能够处理大至 ~10^10^307 的数字（对于双精度变量）。
还有一些关于代码中使用的 Fortran 语法的注意事项：
您的代码中使用了变量
RESMAX，但未进行初始化。给变量赋值时，一定要适当地指定字面常量的种类，否则可能会影响程序结果。例如，这是在调试模式下使用英特尔 Fortran 编译器 2018 编译的原始代码的输出：
-0.152581477302743 7.31503025786548 259.608693509165 -3.152934473473579E-002 99474.1999921620这里是相同代码的输出，但所有文字常量都以种类参数
_dp为后缀（请参阅下面代码的修订版本）：-0.152580456940175 7.31501855886952 259.608692604963 -8.731149137020111E-011 99474.2302854451此答案中修改后的代码输出与上述问题中原始代码的输出略有不同。
不需要用
.gt.、.ge.、.le.、.lt.、……来比较。据我所知，这些是遗留的 FORTRAN 语法。使用更具吸引力的符号（<、>、<=、>=、==）进行比较。没有必要在 Fortran 程序中使用
GOTO语句。这又是遗留的 FORTRAN。通常，简单优雅的 do-loops 和 if-blocks 可以替换GOTO语句，就像下面修改后的代码一样。在 Fortran 中不再需要使用 kind-specific 内部函数（例如
dexp、dlog、...双精度）。在当前的 Fortran 标准中，几乎所有（也许是所有）Fortran 内部函数都有通用名称（exp、log、...）。
以下是对本题程序的修改，解决了上述所有过时的语法，以及处理极大或极小positive-definite变量的问题（我可能走得太远了在 log-transforming 一些永远不会导致溢出或下溢的变量中，但我在这里的目的只是展示 log-transformation 背后的逻辑 positive-definite 变量以及如何处理它们的算法而不可能导致 overflow/underflow/error_in_results).
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