时间戳

首页 > 解决方案 > 使用两个条件在 Fortran 中编写条件循环

fortran - 使用两个条件在 Fortran 中编写条件循环

问题描述

希望有人可以帮助我。我刚刚被介绍给 Fortran，似乎无法弄清楚为什么我的代码会产生无限循环。

我想编写一个代码，在区间 [2,3] 之间找到函数 f(x)= x^3 - 3x - 4 的根 (c)：

我希望步骤是：初始化a和b。然后计算 c = (a+b)/2。然后如果 f(c) < 0，设置 b=c 并重复上一步。如果 f(c) > 0，则设置 a=c 并重复上一步。

关键是重复这些步骤，直到我们得到接近实际根的 1e-4。

这就是我到目前为止所写的内容，它是否会产生无限循环。

我也对使用两个条件循环是否是一个好主意感到困惑（因为函数必须大于/小于 0 .AND。函数的绝对值必须小于 1e-4）。

任何帮助/提示将不胜感激！

我的代码：

PROGRAM proj

IMPLICIT NONE

REAL :: a=2.0, b=3.0, c, f
INTEGER :: count1

c = (a + b)/2

f = c**3 - 3c - 4


DO
   IF (( f .GT. 0.0) .AND. ( ABS(f) .LT. 1e-4)) EXIT

   c = (a+c)/2
   f = c**3 - 3c - 4
   count1 = count1 + 1

    PRINT*, f, c,count1
END DO

PRINT*, c, f

END PROGRAM proj

我希望能够显示迭代并打印每个步骤（更接近实际的根）。

标签： fortranconditional-statementslogical-operators

解决方案

您所描述的是用于在区间 [a:b] 中定位函数的零的二分法。有三种可能性。

  1. 区间不包含零。
  2. 区间的端点为零。
  3. 区间中有多个零。

该程序实现了对多个子区间进行检查的二分法。还有其他更好的方法，但这对您来说应该是可以理解的。

!
! use bisection to locate the zeros of a function f(x) in the interval
! [a,b].  There are three possibilities to consider:  (1) The interval
! contains no zeros;  (2) One (or both) endpoints is a zero; or (3)
! more than one point is a zero.
!
program proj

   implicit none

   real dx, fl, fr, xl, xr
   real, allocatable :: x(:)
   integer i
   integer, parameter :: n = 1000

   xl = 2                     ! Left endpoint
   xr = 3                     ! Right endpoint
   dx = (xr - xl) / (n - 1)   ! Coarse increment

   allocate(x(n))
   x = xl + dx * [(i, i=0, n-1)]    ! Precompute n x-values
   x(n) = xr                        ! Make sure last point is xr
   !
   ! Check end points for zeros.  Comparison of a floating point variable
   ! against zero is exact.
   !
   fl = f(xl)
   if (fl == 0) then 
      call prn(xl, fl)
      x(1) = x(1) + dx / 10         ! Nudge passed xl
   end if
   fr = f(xr)
   if (fr == 0) then 
      call prn(xr, fr)
      x(n) = x(n) - dx / 10         ! Reduce upper xr
   end if
   !
   ! Now do bisection.  Assumes at most one zero in a subinterval.
   ! Make n above larger for smaller intervals.
   !
   do i = 1, n - 1
      call bisect(x(i), x(i+1))
   end do

   contains
      !
      ! The zero satisfies xl < zero < xr
      !
      subroutine bisect(xl, xr)
         real, intent(in) :: xl, xr
         real a, b, c, fa, fb, fc
         real, parameter :: eps = 1e-5
         a = xl
         b = xr
         do
            c = (a + b) / 2
            fa = f(a)
            fb = f(b)
            fc = f(c)
            if (fa * fc <= 0) then        ! In left interval
               if (fa == 0) then          ! Endpoint is a zero.
                  call prn(a, fa)
                  return
               end if
               if (fc == 0) then          ! Endpoint is a zero.
                  call prn(c, fc)
                  return
               end if
               !
               ! Check for convergence.  The zero satisfies a < zero < c.
               !
               if (abs(c - a) < eps) then
                  c = (a + c) / 2
                  call prn(c, f(c))
                  return
               end if
               !
               ! Contract interval and try again.
               !
               b = c
            else if (fc * fb <= 0) then   ! In right interval
               if (fc == 0) then          ! Endpoint is a zero.
                  call prn(c, fc)
                  return
               end if
               if (fb == 0) then          ! Endpoint is a zero.
                  call prn(b, fb)
                  return
               end if
               !
               ! Check for convergence.  The zero satisfies c < zero < b.
               !
               if (abs(b - c) < eps) then
                  c = (b + c) / 2
                  call prn(c, f(c))
                  return
               end if
               !
               ! Contract interval and try again.
               !
               a = c
            else
               return   ! No zero in this interval.
            end if
         end do
      end subroutine bisect


      elemental function f(x)
         real f
         real, intent(in) :: x
         f = x**3 - 3 * x - 4
      end function f

      subroutine prn(x, f)
         real, intent(in) :: x, f
         write(*,*) x, f
      end subroutine prn
end program proj

推荐阅读