Блог

Надеялся, что кто-нибудь сможет мне помочь. Я только что познакомился с Fortran и, похоже, не могу понять, почему мой код создает бесконечный цикл.

Я хочу написать код, который находит корень (c) функции f(x)= x^3 — 3x — 4 между интервалами [2,3]:

Я хочу, чтобы шаги были следующими: инициализируйте a и b. Затем вычислите c = (a b)/2. Тогда, если f(c) Если f(c) > 0, то установите a=c и повторите предыдущий шаг.

Смысл в том, чтобы повторять эти шаги до тех пор, пока мы не приблизимся на 1e-4 к фактическому корню.

Это то, что я написал до сих пор, и это создает бесконечный цикл.

Я также не понимаю, стоит ли использовать цикл с двумя условиями (так как функция должна быть больше/меньше 0 .И. абсолютное значение функции должно быть меньше 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
 

Я хочу иметь возможность показывать итерации и печатать каждый шаг (приближаясь к фактическому корню).

То, что вы описали, — это метод деления пополам для локализации нуля функции в интервале [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