1
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

I'm writing a program that solves the root of a given function f(x). At this point, I just want to check if I'm getting right value for each step and I noticed some errors. To make it simple, I have f(x) = x as the function and I expected the program to print f(a) = 1.0 and f(b) = 1.0 but the output is f(a) = 1.8750000000000000 and f(b) = 0.0000000000000000. Plus although I set both a and b equal to x0, it seems to be a = 1.0000002381857485 and b = 1.0000000000000000. Could anyone explain why this is happening? I might be missing something stupid, but I cannot find it. I appreciate your help.

mike
  • 59
  • 1
  • 1
  • 7
  • Replace the line `end program` with one containing the word `contains`. At the end of the source file containing all your code add the line `end program`. Then, immediately after the line beginning `program` (so after that first line) insert the line `implicit none`. Then (re-)compile and see what you get told by the compiler. – High Performance Mark Oct 17 '16 at 17:47
  • I defined x0 = 1.0d0, a = x0, b = x0 in main program and now it works as I intended. Both a and b print exactly 1.0. But I still have a problem with the function. fa and fb print 0.0000000000000000. I'm not sure why the function works this way.. – mike Oct 17 '16 at 17:55
  • Alternatively to the `contains` approach you could just add `implicit none` to the subroutine (and elsewhere, ideally). Doing that will help you conclude that `f` in the subroutine is a function with `real` not `double precision` result. Naturally, it should be the latter. – francescalus Oct 17 '16 at 18:02
  • I tried 'implicit none' and it gives me "f has not implicit type" but I indicated x and f under function f(x) are double precision though.. How can I change f to a function with double precision result? – mike Oct 17 '16 at 18:12
  • See http://stackoverflow.com/questions/24337413/where-to-put-implicit-none-in-fortran – Vladimir F Героям слава Oct 17 '16 at 18:34
  • 1
    In addition to the link about `implicit none` above, we have documentation on [program structure](https://stackoverflow.com/documentation/fortran/2203/program-units-and-file-layout). Essentially, if you don't go for the `contains` approach mentioned first, then `f` in the subroutine needs to be appropriately declared: `double precision f` (which gives it an implicit interface). Without the `contains` containing both the function and subroutine they are both external procedures and know nothing about each other without explicit detail being provided. – francescalus Oct 17 '16 at 19:23

1 Answers1

0

Add these lines to your subroutine

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

then you should get the correct value for f(x).

Edit:

The interface makes the function f(x) visible (and usable) in the subroutine.

Alex
  • 779
  • 7
  • 15