can function subprograms have arrays?.spot the error

[baazuka]

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

[davekw7x]

Originally Posted by baazuka

i have been given this fortran

The function for which we are finding zeros is

    f(x,k) = tan(x) - x / k

where

    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:

    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)
        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.

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