Bisection метод не сходится

Question

Я пытаюсь реализовать Bisection метод с Fortran 90, чтобы получить решение, с точностью до 10^-5 для 3x - е^х = 0 для 1 < = х < = 2

Это код, который я придумал, но когда я запускаю код, он просто перечисляет .1.5000000000000000 100 раз.

Как исправить этот код, чтобы я мог правильно применять метод биссекции до тех пор, пока не получу число вокруг 10^-5?

program main 
implicit none 
double precision a,b, TOL 
integer maxN 
a = 1.d0 
b = 2.d0 
TOL = .000001d0 
maxN = 100 
call bisection(a, b, TOL, maxN) 
end program 

function f(x) 
double precision x,f 
f = 3*x - e**x 
end function 

subroutine bisection(a, b, TOL,maxN) 
implicit none 
double precision a, b, TOL 
integer maxN 

double precision p, fp, fa, fb, f 
integer j 

p = (a+b)/2.d0 
do j =1, maxN 
fp = f(p) 
fa = f(a) 
fb = f(b) 

    if (0.5d0 * (b-a) < TOL) then 
      print *, "Reach desired tolerance",p 
      return 

    end if 

    if (fa*fp <0) then 
      b = p 
    else if (fb*fp < 0) then 
      a = p 
    end if 

    p = (a+b)/2.d0 
    print *,j,p 
end do 
end subroutine bisection 

Answer 1

Проблема на самом деле не с алгоритмом, а в том, как рассчитать f. Поскольку вы не указали implicit none в функции, компилятор разрешил e**x проскользнуть, хотя Fortran не определяет e, как вам бы хотелось.

При коррекции функции следующим образом, программа работает отлично:

function f(x) 
    implicit none 
    double precision x,f 
    f = 3*x - exp(x) 
end function f 

Это хороший урок, чтобы использовать неявное ни везде.