Numerical Recipes Forum  

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

Reply
 
Thread Tools Display Modes
  #1  
Old 12-27-2011, 09: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




PROGRAM NEWTON

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

theeta(1)=0.1
EPS=0.00001

CALL NEWT(theeta,X,EPS)
STOP
END

!-------------------------------------------------!

SUBROUTINE NEWT(theeta,X,EPS)
real theeta(10),f(100),df(100)



integer L

do L=1,10

X=theeta(L)-(F(theeta(L))/DF(theeta(L)))




DO WHILE(ABS(X-theeta(L)) .GT. EPS)
theeta(L)=X
X=theeta(L)-(F(theeta(L))/DF(theeta(L)))
END DO
end do

RETURN
END




!--------------------------------

function theeta(n)

integer n(10)

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

a=0.1
h=23.0
sigma=46.0


k=1-(a*h/sigma)

do i=1,10

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


end do



return

end

!-------------------------------------------------------------------------


FUNCTION F(theeta)

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

integer m

a=0.1
h=23.0
sigma=46.0




do m=2,10

F(m)=theeta(m)*(1/(1-(a*h/sigma)))
enddo


WRITE(6,10)theeta,F
10 FORMAT(22X,F16.7,4X,F16.7)
RETURN
END

!-------------

FUNCTION DF(theeta)
real df(9),theeta(9)

integer j

do j=2,10

DF(j)=1/(cos(theeta(j)**2))

end do

RETURN
END
Reply With Quote
  #2  
Old 12-28-2011, 01:22 PM
davekw7x davekw7x is offline
Registered User
 
Join Date: Jan 2008
Posts: 453
Quote:
Originally Posted by baazuka View Post
i have been given this fortran
The function for which we are finding zeros is
Code:
    f(x,k) = tan(x) - x / k
where
Code:
    x is "theta"
        and
    k is 1.0 - (a * h / sigma) for given physical parameters a, h, 
sigma
To use the Newton method, we need the derivative of f with respect to x.
I'll call it df:
Code:
    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:
    Code:
         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.

Code:
      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)
        enddo

        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

        stop
      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
Regards,

Dave

Footnote:
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.

Anyhow...

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.
Code:
   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 01:01 PM.
Reply With Quote
Reply

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 04:00 PM.


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