Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Fortran 90 Programming with NR

Thread Tools Display Modes
Old 12-27-2011, 10:43 AM
baazuka baazuka is offline
Registered User
Join Date: May 2011
Posts: 2
can function subprograms have arrays?.spot the error

i have been given this fortran assignment..but geting errors
please do the corrections for me..i have tried allot but failed...

i haave attached the question paper for your reference

heres the coode

!hello..this is a questions of newton iteration
!in this i need to generate the first 10 roots of the given equation F=theeta*(1/(1-(a*h/sigma)))

!as you can see that the equation is a function of theeta and i also been given an equation for theeta

! theeta=(1-1/(((n-0.5)**2)*(3.142**2)*(k**2)))*3.142*(n-0.5)

!annd now you can see that theeta is a function of n..that is the number of roots..for for the first root.i have n=1 and so on..10

!also i need to have the derivative of the function because it will be needed for the newton iteration

!DF=1/(cos(theeta**2)) as u c that it depends on theeta as well

! so u see that i need to have the roots(theeta) for each value of n..and all these values of rooots need to be stored in an array

!then having done all this..these values of roots will be needed to calculate the values of the fnal equation which i have ommited for now

! i have made three function subprograms for these equations that i have mentioned.and this subroutine..but i am getting errors and warnings

!please make corrections and add comments for identification so that i can know my mistakes


real theeta(10),f(100),df(100),eps




real theeta(10),f(100),df(100)

integer L

do L=1,10


DO WHILE(ABS(X-theeta(L)) .GT. EPS)
end do



function theeta(n)

integer n(10)

real k,a,h,sigma,theeta
dimension theeta(10)



do i=1,10


end do




FUNCTION F(theeta)

real a,h,sigma,f(9),theeta(9)

integer m


do m=2,10


10 FORMAT(22X,F16.7,4X,F16.7)


real df(9),theeta(9)

integer j

do j=2,10


end do

Reply With Quote
Old 12-28-2011, 02:22 PM
davekw7x davekw7x is offline
Registered User
Join Date: Jan 2008
Posts: 453
Originally Posted by baazuka View Post
i have been given this fortran
The function for which we are finding zeros is
    f(x,k) = tan(x) - x / k
    x is "theta"
    k is 1.0 - (a * h / sigma) for given physical parameters a, h, 
To use the Newton method, we need the derivative of f with respect to x.
I'll call it df:
    df(x, k) = sec(x)**2 - 1.0/k
For each root, here's the drill:
  • Get an estimate of the root
  • Set x equal to the estimate
  • Apply the following iteration for a certain number of times:
         x = x - f(x, k)/df(x, k)
    If the absolute value of f(x,k) is less than some specified
    value (epsilon), quit iterating and accept the value of x
    as a reasonable estimate of the root.

You don't actually need an array of estimates or an array of roots, but you can do it if you want to. Similarly, you can put the roots in an array if you want to.

Here's the skeleton of a program that puts the estimates in an array. The subroutine that calculates estimates fills an array that was one of its parameters.

      program NewtonDemo
        ! I like to make sure that all variables are declared
        implicit none

        double precision estimates(10)
        double precision a, h, sigma
        double precision est, x
        double precision f
        double precision eps
        double precision k
        integer n

        ! Upper limit on absolute error used to terminate
        ! Newton iteration
        eps = 0.0000001

        ! Problem physical parameters
        a = 0.1
        h = 23.0
        sigma = 46.0

        ! To keep from passing a, h, sigma to functions and subroutines
        ! define this constant
        k = 1.0 - (a*h/sigma)

        ! For convenience, calculate ten estimates and put them
        ! into an array
        call CalcTenEstimates(estimates, k)

        write(*, '("First estimates of the roots to be used by Newton are")')
        do n=1,10
          write(*, &
          '("   Estimate for root number ",I2," is x = ", 1pe14.7)') &
           n, estimates(n)

        print *

        do n = 1,10
          est = estimates(n)
          call Newton(est, x, eps, k)
          write(*, &
          '("   Root number ", I2, ": x = ", 1pe14.7, ", f(x) = ", 1pe14.7)') &
          n, x, f(x, k)
        end do

      end program NewtonDemo


      ! This is the way that the book suggested to calculate
      ! estimates of the roots.  Each estimate will be "refined"
      ! by calling the Newton subroutine
      subroutine CalcTenEstimates(th, k)
        implicit none

        ! Declare types of the parameter variables
        double precision th(10), k

        ! Local constant
        double precision Pi
        data Pi  /3.14159265358979323846/

        ! Local variable
        integer n

        ! The book suggests a special procedure for root number 1.
        ! Heck, I'll just use the formula he gave for all of them.
        ! Maybe I'll get lucky.

        ! The book suggests this formula for estimates 2, 3, ...
        do n = 1,10
            th(n)=(n -0.5) * Pi * (1.0-1.0/(((n-0.5)**2)*(Pi**2)*(k**2)))
        end do

        ! The subroutine returns by falling off of the end...
      end subroutine CalcTenEstimates


      subroutine Newton(est, x, eps, k)
        implicit none

        ! Declare types of parameter variables
        double precision est, x, eps, k

        ! Type declarations for functions called from this subroutine
        double precision df, f
        ! Set x = est, then
        ! do the iteration x = x-f(x,k)/df(x,k) for a maximum of,
        ! say, 20 times.  Break out of the loop if the
        ! absolute value of f(x,k) gets smaller than eps

      end subroutine Newton

      ! Derivative of the function for which we are getting roots
      function df(th, k)
        implicit none

        double precision df, th, k

        ! Calculate the derivative for a given value of theta and k
        df = whatever...
      end function df

      ! The function for which we are getting roots
      function f(th, k)
        implicit none

        double precision f, th, k

        ! Calculate the value of the function for a given value of
        ! theta and k
        f = whatever...
      end function f


The book states that the first root can be seen to be close to 0.4. By "inspection," the first zero is at zero. (Plug in zero and see what the function value is.) As a programming exercise, try the method he suggested for the next one, instead of just plowing ahead with his formula as my program did.


Here are the results I got from the program I showed (didn't bother with the root at zero).

Also, I can't see any reason to use single precision arithmetic, so I used doubles everywhere.

Maybe someone can check the values I got.
   Estimate for root number  1 is x =  8.6540052E-01
   Estimate for root number  2 is x =  4.4772572E+00
   Estimate for root number  3 is x =  7.7129027E+00
   Estimate for root number  4 is x =  1.0894804E+01
   Estimate for root number  5 is x =  1.4058790E+01
   Estimate for root number  6 is x =  1.7214633E+01
   Estimate for root number  7 is x =  2.0366092E+01
   Estimate for root number  8 is x =  2.3514919E+01
   Estimate for root number  9 is x =  2.6662044E+01
   Estimate for root number 10 is x =  2.9808005E+01

   Root number  1: x =  3.8536818E-01, f(x) =  8.8982981E-09
   Root number  2: x =  4.5045364E+00, f(x) =  6.8410989E-14
   Root number  3: x =  7.7317240E+00, f(x) =  9.4918951E-13
   Root number  4: x =  1.0908707E+01, f(x) =  2.4237921E-12
   Root number  5: x =  1.4069749E+01, f(x) =  3.8825670E-12
   Root number  6: x =  1.7223659E+01, f(x) =  4.6440386E-12
   Root number  7: x =  2.0373757E+01, f(x) =  7.0303624E-12
   Root number  8: x =  2.3521578E+01, f(x) =  8.5806171E-12
   Root number  9: x =  2.6667929E+01, f(x) =  8.4786241E-12
   Root number 10: x =  2.9813276E+01, f(x) =  8.6695685E-12
I could be wrong, you know. It wouldn't be the first time. Not even the first time today.

The real bottom line is that Newton is not guaranteed to converge to a particular root (or even to converge at all). Coming up with a reasonable estimate for a given root is a major part of the exercise...

Last edited by davekw7x; 12-29-2011 at 02:01 PM.
Reply With Quote

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump

All times are GMT -5. The time now is 11:15 AM.

Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.