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

Join the discussion
Registration is free. Ask a follow-up in this thread, or start your own.
1 reply · 12K views
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?
 
Physics 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...