满足条件后，While循环不会停止

Posted 2021-04-16

我已经为Fortran中的bisection方法编写了以下代码，但是在满足我在Do While循环中定义的条件后它不会停止。任何解释将不胜感激。

! Bi-section method
program bi_section 
implicit none

!Variables
real :: a, b, p  
real :: fa, fb, fp, fc

! Interval
a = 2.0
b = 3.0

! Value of fa & fb
fa = a*a*a-13
fb = b*b*b-13

    do while (abs(fb-fa) > 0.00001)
    p = (a+b)/2
    fp = p*p*p-13
        if (fp*fa > 0) then
            a = p
        else
            b = p      
        end if
    print*, p
    enddo
end program bi_section

答案

为了找到可能的错误，一个好的方法是打印相关变量的值（例如a，b，fa和fb）并检查事情是否按预期进行。例如，

do while (abs(fb-fa) > 0.00001)
    p = (a+b)/2
    fp = p*p*p-13
    if (fp*fa > 0) then
        a = p
    else
        b = p      
    end if
    !-- print*, p
    print*, p, a, b, fa, fb    !<--- print more variables for check
enddo

以上打印声明给出

   2.50000000       2.00000000       2.50000000      -5.00000000       14.0000000
   2.25000000       2.25000000       2.50000000      -5.00000000       14.0000000
   2.37500000       2.25000000       2.37500000      -5.00000000       14.0000000
   2.31250000       2.31250000       2.37500000      -5.00000000       14.0000000
   2.34375000       2.34375000       2.37500000      -5.00000000       14.0000000
   2.35937500       2.34375000       2.35937500      -5.00000000       14.0000000
   2.35156250       2.34375000       2.35156250      -5.00000000       14.0000000
   2.34765625       2.34765625       2.35156250      -5.00000000       14.0000000
   2.34960938       2.34960938       2.35156250      -5.00000000       14.0000000
   2.35058594       2.35058594       2.35156250      -5.00000000       14.0000000
   2.35107422       2.35107422       2.35156250      -5.00000000       14.0000000
   2.35131836       2.35131836       2.35156250      -5.00000000       14.0000000
   2.35144043       2.35131836       2.35144043      -5.00000000       14.0000000
   2.35137939       2.35131836       2.35137939      -5.00000000       14.0000000
   2.35134888       2.35131836       2.35134888      -5.00000000       14.0000000
   2.35133362       2.35133362       2.35134888      -5.00000000       14.0000000
   2.35134125       2.35133362       2.35134125      -5.00000000       14.0000000
   2.35133743       2.35133362       2.35133743      -5.00000000       14.0000000
   2.35133553       2.35133362       2.35133553      -5.00000000       14.0000000
   2.35133457       2.35133457       2.35133553      -5.00000000       14.0000000
   2.35133505       2.35133457       2.35133505      -5.00000000       14.0000000
   ...

因此，虽然p正在接近正确的解决方案（立方根13 = 2.3513346 ......），但fa和fb都没有改变，并且“条件”（abs(fb-fa) <= 0.00001）从未得到满足（正如@HighPerformanceMark在评论中所建议的那样）。为了避免这种情况，我们需要更新fa和fb，例如，

do while (abs(fb-fa) > 0.00001)
    p = (a+b)/2
    fp = p*p*p-13
    if (fp*fa > 0) then
        a = p
        fa = fp  !<---
    else
        b = p
        fb = fp  !<---
    end if
    print*, p, a, b, fa, fb
enddo

这使

   2.50000000       2.00000000       2.50000000      -5.00000000       2.62500000
   2.25000000       2.25000000       2.50000000      -1.60937500       2.62500000
   2.37500000       2.25000000       2.37500000      -1.60937500      0.396484375
   2.31250000       2.31250000       2.37500000     -0.633544922      0.396484375
   2.34375000       2.34375000       2.37500000     -0.125396729      0.396484375
   2.35937500       2.34375000       2.35937500     -0.125396729      0.133815765
   2.35156250       2.34375000       2.35156250     -0.125396729       3.77845764E-03
   2.34765625       2.34765625       2.35156250      -6.09159470E-02   3.77845764E-03
   2.34960938       2.34960938       2.35156250      -2.85959244E-02   3.77845764E-03
   2.35058594       2.35058594       2.35156250      -1.24149323E-02   3.77845764E-03
   2.35107422       2.35107422       2.35156250      -4.32014465E-03   3.77845764E-03
   2.35131836       2.35131836       2.35156250      -2.70843506E-04   3.77845764E-03
   2.35144043       2.35131836       2.35144043      -2.70843506E-04   1.75380707E-03
   2.35137939       2.35131836       2.35137939      -2.70843506E-04   7.41004944E-04
   2.35134888       2.35131836       2.35134888      -2.70843506E-04   2.35557556E-04
   2.35133362       2.35133362       2.35134888      -1.71661377E-05   2.35557556E-04
   2.35134125       2.35133362       2.35134125      -1.71661377E-05   1.08718872E-04
   2.35133743       2.35133362       2.35133743      -1.71661377E-05   4.48226929E-05
   2.35133553       2.35133362       2.35133553      -1.71661377E-05   1.33514404E-05
   2.35133457       2.35133457       2.35133553      -1.90734863E-06   1.33514404E-05
   2.35133505       2.35133457       2.35133505      -1.90734863E-06   6.67572021E-06

（顺便说一下，我们可以把p*p*p-13写成p**3 - 13）。