can function subprograms have arrays?.spot the error

baazuka
12-27-2011, 09:43 AM
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

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
12-28-2011, 01:22 PM
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...