Fortran 90 странное поведение

Question

program main 
call findbracket(x0, a, b) 
end program 

function f(x) 
double precision x,f 
f = x 
end function 

subroutine findbracket(x0,a,b) 
double precision x0, a, b 
double precision fa, fb 
double precision dx 
dx = 0.001d0 
x0 = 1.0d0 
a = x0 
b = x0 
print*, a, b 
print*, f(a) 
print*, f(b) 

do 
fa = f(a) 
fb = f(b) 

print*, "what is fa", fa 
print*, "what is fb", fb 
a = a - dx 

if (fa*fb < 0) then 
     exit 
end if 
print*, b, dx 
b = b + dx 
if (fa*fb < 0) then 
     exit 
end if 

dx = dx*2 

end do 

end subroutine 

Я пишу программу, которая решает корень данной функции f (x). На этом этапе я просто хочу проверить, не получаю ли я правильное значение для каждого шага, и я заметил некоторые ошибки. Чтобы сделать его простым, у меня есть функция f (x) = x, и я ожидал, что программа будет печатать f (a) = 1.0 и f (b) = 1.0, но выход будет f (a) = 1,8750000000000000 и f (b) = 0,0000000000000000. Плюс, хотя я устанавливаю как a, так и b равным x0, он выглядит как a = 1.0000002381857485 и b = 1.0000000000000000. Может ли кто-нибудь объяснить, почему это происходит? Возможно, мне не хватает чего-то глупого, но я не могу найти его. Я ценю вашу помощь.

Answer 1

Добавьте эти строки в подпрограмме

subroutine findbracket(x0,a,b) 
double precision x0, a, b 
double precision fa, fb 
double precision dx 

interface 
    function f(x) 
    double precision x,f 
    end function 
end interface 

то вы должны получить правильное значение для F (х).

Редактировать:

Интерфейс делает функцию F (X) видимый (и годная к употреблению) в подпрограмме.