How Can I Correct My Fortran Program for Newton's Method?

Click For Summary
SUMMARY

The discussion centers on correcting a Fortran program implementing Newton's method for solving equations. The original code incorrectly uses the variable x0 in the subroutine FCN instead of xn, leading to infinite values in the output. The suggested correction involves changing the function definitions within FCN to use xn for calculations, which resolves the issue. Additionally, adjusting the tolerance (TOL) and epsilon (EPS) values may improve convergence.

PREREQUISITES
  • Understanding of Newton's method for root finding
  • Familiarity with Fortran programming language
  • Knowledge of function and subroutine definitions in Fortran
  • Basic concepts of numerical tolerance and convergence criteria
NEXT STEPS
  • Implement the corrected subroutine FCN using xn for calculations
  • Research numerical methods for root finding beyond Newton's method
  • Explore debugging techniques in Fortran to identify variable scope issues
  • Learn about numerical stability and convergence criteria in iterative methods
USEFUL FOR

Fortran developers, numerical analysts, and anyone interested in implementing or debugging numerical methods for solving equations.

jjohnson8
Messages
1
Reaction score
0
I'm working on a program for Newton's method for solving equations.

This is my code:

=======================================================
program Newton
implicit real(a-h,o-z)
F(x) = x**2 - 4
!...&---1---------2---------3---------4---------5---------6---------7---

!----All values to be read in/set values----!
PRINT*, 'Please inter initial guess for root: '
READ*, x0
maxIT = 20
TOL = 1.e-14
EPS = 1.e-14
xn = 0
!----Print titles for output----!
PRINT 20, 'x0', 'maxIT', 'TOL', 'EPS', 'n', 'xn', 'F(xn)'

!----Do Loop to go from 0 to max iterations----!
DO n = 0,maxIT
CALL FCN(xn, Fn, DFn)
xn = x0 - (Fn/DFn)
Dx = xn - x0
x0 = xn
IF( ABS(Dx) .LE. TOL) THEN
IF( ABS(Fn) .LE. TOL) THEN
PRINT*, 'DONE: root=',xn,' found in',n,'iterations'
PRINT*, ' residual =', Fn
ELSE
PRINT*, 'Stuck at iteration:', n
PRINT*, ' relDX < EPS but residual=',Fn,' > TOL, exiting'
STOP
ENDIF
ENDIF
PRINT 10, x0, maxIT, TOL, EPS, n, xn, F(xn)
ENDDO
!----Formats for headings----!
10 FORMAT(F10.7,4x,I2,4x,F18.14,13x,F18.14,8x,I2,3x,F10.7,3x,F10.7)
20 FORMAT(A2,10x,A5,8x,A3,35x,A3,10x,A1,6x,A2,11x,A5)
END

subroutine FCN(xn, Fn, DFn)
Fn = (x0**2) - 4
DFn = 2*x0
return
END


Here is my output:

x0 maxIT TOL EPS n xn F(xn)
+Infinity 20 0.00000000000001 0.00000000000001 0 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 1 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 2 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 3 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 4 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 5 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 6 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 7 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 8 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 9 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 10 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 11 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 12 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 13 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 14 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 15 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 16 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 17 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 18 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 19 +Infinity +Infinity
+Infinity 20 0.00000000000001 0.00000000000001 20 +Infinity +Infinity


I've tinkered with this all day, and I can't figure out how to get anything different than infinity for the values. Anyone out there that can correct me?
 
Technology news on Phys.org
Looks like you're not passing xn, the function should be

subroutine FCN(xn, Fn, DFn)
Fn = (xn**2) - 4
DFn = 2*xn
return

also try making tol and eps larger...
 

Similar threads

  • · Replies 4 ·
Replies
4
Views
4K
  • · Replies 11 ·
Replies
11
Views
2K
  • · Replies 6 ·
Replies
6
Views
3K
  • · Replies 6 ·
Replies
6
Views
2K
  • · Replies 4 ·
Replies
4
Views
9K
  • · Replies 6 ·
Replies
6
Views
2K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 6 ·
Replies
6
Views
4K
  • · Replies 15 ·
Replies
15
Views
5K
  • · Replies 6 ·
Replies
6
Views
5K