What's causing my Fortran bisection method to give incorrect results?

Click For Summary
SUMMARY

The forum discussion centers on a Fortran implementation of the bisection method to solve a transcendental equation for temperature in thermodynamics. The user encountered an issue where the program incorrectly reported "There's an even number of root(s) in that interval" due to the function not properly utilizing the temperature variables t_1 and t_2. After incorporating these variables as arguments in the function definition, the problem was resolved, although the user initially received an incorrect solution of "3.02481484E-39". The final solution was achieved by correcting a variable reference in the output statement.

PREREQUISITES
  • Understanding of Fortran programming language
  • Familiarity with numerical methods, specifically the bisection method
  • Knowledge of transcendental equations and their solutions
  • Ability to debug and troubleshoot code
NEXT STEPS
  • Learn how to pass arguments to functions in Fortran
  • Explore debugging techniques in Fortran to identify variable scope issues
  • Study the implementation of numerical methods in Fortran for solving equations
  • Investigate the behavior of transcendental functions and their roots
USEFUL FOR

This discussion is beneficial for Fortran programmers, students learning numerical methods, and anyone involved in computational thermodynamics or debugging scientific code.

fluidistic
Gold Member
Messages
3,932
Reaction score
283
I'm rusty with fortran and programming in general. I can't see my "error" in a code that I wrote from scratch.
Basically I wanted to get some fun and solve for a temperature in thermodynamics where I must get "T_f" which appear in a transcendental equation: ##A\ln \left ( \frac{T_f^2}{T_1T_2} \right )+2B[2T_f-(T_1T_2)]=0##.
Where T_1 is worth 100, T_2 is worth 200, A=2 and B=0.005.
Here is my code:
Code:
program rootfinderbisection
implicit none

real(4)::a,b,xn,fa,fb,fn,f,e,x_n,t_1,t_2


write(*,*)"Enter the temperatures t_1 and t_2"
read(*,*)t_1, t_2
write(*,*)"Enter your interval guess, [a,b]"
read(*,*)a,b
write(*,*)"Enter the precision you want"
read(*,*)e
fa=f(a)
fb=f(b)

if (abs(fa)<e) then
write(*,*) "The solution is",a
stop
else if (abs(fb)<e) then
write(*,*)"The solution is",b
stop
end if

if (f(a)*f(b)>0) then
write(*,*)"There's an even number of root(s) in that interval"
stop
endif

xn=0.5*(a+b)
fn=f(xn)

do while (abs(fn)>e) 
if (fa*fn<0) then
b=xn
fb=fn
else
a=xn
fa=fn
endif

xn=0.5*(a+b)
fn=f(xn)
end do
write(*,*)"The solution is",x_n


end program
that I compile along with
Code:
function f(x)
implicit none

real(4)::f,x, t_1,t_2
f=2.*log((x**2.)/(t_1*t_2))+2.*0.005*(2.*x-(t_1+t_2))

end function
The problem is that for the values I said above, I get the message "There's an even number of root(s) in that interval". Thus apparently the condition f(a)*f(b)>0 is satisfied. However when I enter the numbers in my pocket calculator I get f(a)=-2.3862... and f(b)=2.3862... where I truncated both values. Thus the condition is not satisfied.
I really don't see what's wrong. Any idea?

Edit: I tried with many different values for a and b. I even changed t_1 and t_2, same message.

Edit2: I've just tested the values of fa and fb. It thinks they are worth infinity. Going to check this out.
Edit3: No idea why it thinks they are worth infinity. o_0.
Edit4: ok it does not plug in t_1 and t_2 inside my function.
 
Last edited:
Technology news on Phys.org
Try including t_1 and t_2 as arguments in the function f and see if this resolves the calculation errors.
 
Your function "f(x)" does not agree with your formula at the top of your post. Is the last term in the bracket term of your equation "T_1 * T_2" or "T_1 + T_2".
 
fluidistic said:
Edit4: ok it does not plug in t_1 and t_2 inside my function.

t_1 and t_2 in the function are completely separate variables from t_1 and t_2 in the main program. If you want t_1 and t_2 in the function to have the same values as in main, you need to make them arguments to the function along with x.
 
Thanks all for the tips. To TheoMcCloskey, you're right, I made a latex typo here for T_f. The one in my program is ok.
So it got rid of the problem! However I now get "The solution is 3.02481484E-39" which is basically 0. I'll investigate and post if stuck.

Edit: Worked, problem solved. I had to change the line "write(*,*)"The solution is",x_n" to "write(*,*)"The solution is",xn".
Thanks a lot guys!
 
Last edited:

Similar threads

  • · Replies 8 ·
Replies
8
Views
2K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 11 ·
Replies
11
Views
2K
  • · Replies 12 ·
Replies
12
Views
2K
  • · Replies 6 ·
Replies
6
Views
3K
  • · Replies 11 ·
Replies
11
Views
2K
  • · Replies 13 ·
Replies
13
Views
3K
  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 17 ·
Replies
17
Views
23K
  • · Replies 19 ·
Replies
19
Views
3K