Problem with qromb subroutine/nrerror
Author Message
Problem with qromb subroutine/nrerror

Dear all ;

I'm trying to compute an integral with Romberg's Method.
I used qromb subroutine  which calls polint and trapzd.
but I cannot get the answer correctly.

the program stops the running with nrerror .
I checked all parts but I didn't find the mistake.
result which can be printed :

nrerror: qromb:too many steps
STOP program terminated by nrerror

program test2
IMPLICIT NONE
integer, parameter ::  dp = 8
real
(kind=dp),parameter::PI=3.141592653589793238462643383279502884197_dp
u=0
do while(u.lt.20.D0)
t = 1.D0
frac=u/(2*t)
if (u==0.D0) then
else

u = u + 0.1D0
end do
end program test2
!-------------------------------------------------------
!-----------------------------------------------

SUBROUTINE NewTrapzd(a, b, Resu, N, yy)
!use nrtype
implicit none
!real(sp)::resu
REAL(Kind=8) :: A, B, X, YY, Del,func,Resu
INTEGER :: I, N
if(N==1) then
Resu = 0.5D0*(func(A,YY) + func(B,YY))
else
Resu = 0.5D0*(func(A,YY) + func(B,YY))
Del = (B-A)/DBLE(N-1) !h=(b-a)/(n-1)
DO I=1, N-2
X = A + I*Del  !x_i=x0+ih
Resu = Resu + Del*func(X,YY) !f(x_i)*h--> ...
END DO
end if
END SUBROUTINE NewTrapzd

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

real(kind=8)  function func(x,yy)
!use nrtype
implicit none
real(kind=8)::x ,yy
if (abs(yy)<1.D-6)yy=0.1D0
if (abs(x)<1.D-6)then
func=-1/(2*yy)
else if (abs(x-1)<1.D-6)then
func=-1/(4*yy)
else
func=(besj0(-log(x)/yy)*besj1(-log(x)/yy))/((x+1)*log(x))
end if
end function func

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

!FUNCTION qromb(func,a,b)
real(kind=8) FUNCTION qromb(a,b,yy)
USE nrtype; USE nrutil, ONLY :nrerror
!USE nr, ONLY :polint,trapzd
USE nr, ONLY :polint,newtrapzd
IMPLICIT NONE
REAL(kind=8), INTENT(IN) :: a,b
REAL(sp) :: qromb

REAL(kind=8) ::func
REAL(kind=8), INTENT(IN) ::yy

!INTERFACE
!FUNCTION func(x)
!USE nrtype
!REAL(SP), DIMENSION(:), INTENT(IN) :: x
!REAL(SP), DIMENSION(size(x)) :: func
!END FUNCTION func
!END INTERFACE
INTEGER(I4B), PARAMETER :: JMAX=20,JMAXP=JMAX+1,K=5,KM=K-1
REAL(SP), PARAMETER :: EPS=1.0e-6_sp
!Returns the integral of the function func from a to b.
Integration is performed by Romberg's
!method of order 2K, where, e.g., K=2 is Simpson's rule.
!Parameters: EPS is the fractional accuracy desired, as
determined by the extrapolation error
!estimate; JMAX limits the total number of steps; K is the
number of points used in the
!extrapolation.
REAL(sp), DIMENSION(JMAXP) :: h,s !These store the successive
trapezoidal approximations
!and
REAL(SP) :: dqromb !their relative stepsizes.
INTEGER(I4B) :: j
real(kind=8)::ss
h(1)=1.0
do j=1,JMAX
!call trapzd(func,a,b,s(j),j)
! call trapzd(a,b,s(j),j,yy)
!call NewTrapzd(a,b,s(j),j,yy)
call NewTrapzd(a,b,ss,j,yy)
if (j >= K) then
call polint(h(j-KM:j),s(j-KM:j),0.0_sp,qromb,dqromb)
if (abs(dqromb) <= EPS*abs(qromb)) RETURN
end if
s(j+1)=s(j)
h(j+1)=0.25_sp*h(j) !This is a key step: The factor is 0.25 even
!though the stepsize is decreased by only
!0.5. This makes the extrapolation a polynomial
!in h2 as allowed by equation (4.2.1),
!not just a polynomial in h.
end do
call nrerror('qromb:too many steps')
END FUNCTION qromb

!-------------------------------------
SUBROUTINE polint(xa,ya,x,y,dy)
USE nrtype; USE nrutil, ONLY : assert_eq,iminloc,nrerror
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya
REAL(SP), INTENT(IN) :: x
REAL(SP), INTENT(OUT) :: y,dy
!Given arrays xa and ya of length N, and given a value x, this
routine returns a value y,
!and an error estimate dy. If P(x) is the polynomial of degree N
- 1 such that P(xai) =
!yai, i = 1, . . . ,N, then the returned value y = P(x).
INTEGER(I4B) :: m,n,ns
REAL(SP), DIMENSION(size(xa)) :: c,d,den,ho
n=assert_eq(size(xa),size(ya),'polint')
c=ya !Initialize the tableau of c's and d's.
d=ya
ho=xa-x
ns=iminloc(abs(x-xa)) !Find index ns of closest table entry.
y=ya(ns) !This is the initial approximation to y.
ns=ns-1
do m=1,n-1 !For each column of the tableau,
den(1:n-m)=ho(1:n-m)-ho(1+m:n) !we loop over the current c's and
d's and up-
if(any(den(1:n-m) == 0.0)) &   !date them.
call nrerror('polint calculation failure')
!This error can occur only if two input xa's are (to within
roundoff) identical.
den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
d(1:n-m)=ho(1+m:n)*den(1:n-m) !Here the c's and d's are updated.
c(1:n-m)=ho(1:n-m)*den(1:n-m)
if (2*ns < n-m) then !After each column in the tableau is
completed, we decide
!which correction, c or d, we want to add to our accumulating
!value of y, i.e., which path to take through
!the tableau--forking up or down. We do this in such a
!way as to take the most "straight line" route through the
!tableau to its apex, updating ns accordingly to keep track
!of where we are. This route keeps the partial approximations
!centered (insofar as possible) on the target x. The
!last dy added is thus the error indication.
dy=c(ns+1)
else
dy=d(ns)
ns=ns-1
end if
y=y+dy
end do
END SUBROUTINE polint

Fri, 09 Sep 2011 19:28:09 GMT
Problem with qromb subroutine/nrerror

Quote:

>    Dear all ;

> I'm trying to compute an integral with Romberg's Method.
> I used qromb subroutine  which calls polint and trapzd.
> but I cannot get the answer correctly.

> the program stops the running with nrerror .
> me? I checked all parts but I didn't find the mistake.

A little knowledge is a dangerous thing. You have mangled the NR code to such an extent that it cannot do anything useful.

You need to learn how to pass function/subroutines as arguments. This is covered in textbooks/manuals.

Quote:

>    program test2
>       IMPLICIT NONE
>       integer, parameter ::  dp = 8
>       real
> (kind=dp),parameter::PI=3.141592653589793238462643383279502884197_dp
>       u=0
>       do while(u.lt.20.D0)
>         t = 1.D0
>         frac=u/(2*t)
>         if (u==0.D0) then
>       else

>         u = u + 0.1D0
>       end do
>       end program test2

<--CUT rest of code-->

An actual argument of type function/subroutine must be declared 'external' or you need to provide an explicit interface.

-- mecej4

Fri, 09 Sep 2011 20:32:13 GMT
Problem with qromb subroutine/nrerror

Quote:

> > ? ?Dear all ;

> > I'm trying to compute an integral with Romberg's Method.
> > I used qromb subroutine ?which calls polint and trapzd.
> > but I cannot get the answer correctly.

> > the program stops the running with nrerror .
> > me? I checked all parts but I didn't find the mistake.

> A little knowledge is a dangerous thing. You have mangled the NR code to such an extent that it cannot do anything useful.

> You need to learn how to pass function/subroutines as arguments. This is covered in textbooks/manuals.
> > ? ? ? answer=qromb(0.D0,1.D0,frac)
> An actual argument of type function/subroutine must be declared 'external' or you need to provide an explicit interface.

Except that Fatemeh is NOT passing his function as an argument to any
routine. That is the usual way of using a general purpose integraion
routine, but the OP has decided not to do that.

His problem is that his integrand is a function of two variables, one
of which is a "parameter" (as it is used in math).
them, as you say. It would have been far better to add the "parameter"
as an argument to the argument list of existing routines, and pass it
down as needed.

Is it NR or a lack of experience that encourages students to get in

--- e

Fri, 09 Sep 2011 22:53:00 GMT
Problem with qromb subroutine/nrerror

Quote:

> > > ? ?Dear all ;

> > > I'm trying to compute an integral with Romberg's Method.
> > > I used qromb subroutine ?which calls polint and trapzd.
> > > but I cannot get the answer correctly.

> > > the program stops the running with nrerror .
> > > me? I checked all parts but I didn't find the mistake.

> > A little knowledge is a dangerous thing. You have mangled the NR code to such an extent that it cannot do anything useful.

> > You need to learn how to pass function/subroutines as arguments. This is covered in textbooks/manuals.
> > > ? ? ? answer=qromb(0.D0,1.D0,frac)
> > An actual argument of type function/subroutine must be declared 'external' or you need to provide an explicit interface.

> Except that Fatemeh is NOT passing his function as an argument to any
> routine. That is the usual way of using a general purpose integraion
> routine, but the OP has decided not to do that.

> His problem is that his integrand is a function of two variables, one
> of which is a "parameter" (as it is used in math).
> Instead of adapting existing routines that do work, he has mangled
> them, as you say. It would have been far better to add the "parameter"
> as an argument to the argument list of existing routines, and pass it
> down as needed.

Other alternatives would allow testing and compiling the numerical
subroutines to a library, as written.

1. put the parameter in a common block
2. put the parameter in a module

Then modify the code of the function passed as an argument as needed.

[Been there, done that, screwed it up myself!]

---- e

Fri, 09 Sep 2011 23:28:20 GMT
Problem with qromb subroutine/nrerror

Quote:

>> > Dear all ;

>> > I'm trying to compute an integral with Romberg's Method.
>> > I used qromb subroutine  which calls polint and trapzd.
>> > but I cannot get the answer correctly.

>> > the program stops the running with nrerror .
>> > guide me? I checked all parts but I didn't find the mistake.

>> A little knowledge is a dangerous thing. You have mangled the
>> NR code to such an extent that it cannot do anything useful.

>> You need to learn how to pass function/subroutines as
>> arguments. This is covered in textbooks/manuals.

>> An actual argument of type function/subroutine must be declared
>> 'external' or you need to provide an explicit interface.

> Except that Fatemeh is NOT passing his function as an argument
> to any routine. That is the usual way of using a general purpose
> integraion routine, but the OP has decided not to do that.

> His problem is that his integrand is a function of two
> variables, one of which is a "parameter" (as it is used in
> has mangled them, as you say.

He did so without stating his reasons for modifying the NR routines; furthermore, he did not make it explicit that he had changed the NR routines, but asks "Is there anyone who has any experience about this?".

NR, in an earlier edition, gloated about a "bookkeeping swindle" (p.412 of the fortran 77 edition; Google for the phrase) as a way of using a 1-D minimization code for a unidirectional minimization of an n-D function, using a COMMON block variable.

Quote:
> It would have been far better to
> add the "parameter" as an argument to the argument list of
> existing routines, and pass it down as needed.

> Is it NR or a lack of experience that encourages students to get

The environment into which NR was born had a lot to do with both the erstwhile popularity and the risks associated with NR. It arrived at a time when few had network access. You had to go to the library, photocopy pages from journals such as TOMS, and transcribe the Algol-like algorithm listings into working software. NR sailed in, seductively offering an easy path to a dream world of numerical bliss. Many of us took the trip.

-- mecej4

Sat, 10 Sep 2011 00:13:16 GMT
Problem with qromb subroutine/nrerror

Quote:

> >> > Dear all ;

> >> > I'm trying to compute an integral with Romberg's Method.
> >> > I used qromb subroutine ?which calls polint and trapzd.
> >> > but I cannot get the answer correctly.

> >> > the program stops the running with nrerror .
> >> > guide me? I checked all parts but I didn't find the mistake.

> >> A little knowledge is a dangerous thing. You have mangled the
> >> NR code to such an extent that it cannot do anything useful.

> NR, in an earlier edition, gloated about a "bookkeeping swindle" (p.412 of the Fortran 77 edition; Google for the phrase) as a way of using a 1-D minimization code for a unidirectional minimization of an n-D function, using a COMMON block variable.
> > Is it NR or a lack of experience that encourages students to get
> > in over their heads???
> The environment into which NR was born had a lot to do with both the erstwhile popularity and the risks associated with NR. It arrived at a time when few had network access. You had to go to the library, photocopy pages from journals such as TOMS, and transcribe the Algol-like algorithm listings into working software. NR sailed in, seductively offering an easy path to a dream world of numerical bliss. Many of us took the trip.

[OT?] Sometimes it is difficult to understand how otherwise competent
science and engineering students have such trouble with Fortran. There
are *some* obscure rules, but more often the students seem to trip on
the fundamentals. Are there any good textbooks that cover Fortran,
general programming principles and some of the numerical techniques
that would be of use to those students?

---- e

Sat, 10 Sep 2011 02:00:56 GMT
Problem with qromb subroutine/nrerror

Quote:
> I'm trying to compute an integral with Romberg's Method.
> I used qromb subroutine  which calls polint and trapzd.
> but I cannot get the answer correctly.

I haven't looked over your code carefully.  One thing that jumps out
at me is way you are handling KINDs.  SP sounds like the KIND of
single-precision REALs, like SP = KIND(1.0).  Why then are you using
single-precision variables all over the place when your core function
func is double precision?  Also you define double precision constants
such as 1.D-6.  It's much more flexible to define them as 1.0e-6_dp.
That way if you want to change everything between single and double
precision you just have to change the definition of dp in the module
where it's defined (should be nrtype) and your whole program changes
in one easy step.  Also get rid of declarations like real(kind=8) x;
uncomment the use nrtype line and change them to real(kind=dp) x.
If the KIND of an actual argument and the associated dummy argument
don't agree, you get unpredictable results.  Setting the KIND of
variables throughout your program in a consistent way can give more
confidence that you aren't making this mistake.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Sat, 10 Sep 2011 02:21:27 GMT
Problem with qromb subroutine/nrerror
(someone wrote)

Quote:
>> Is it NR or a lack of experience that encourages
>> students to get in over their heads???
> The environment into which NR was born had a lot to do
> with both the erstwhile popularity and the risks associated
> with NR. It arrived at a time when few had network access.
> You had to go to the library, photocopy pages from journals
> such as TOMS, and transcribe the Algol-like algorithm
> listings into working software. NR sailed in, seductively
> offering an easy path to a dream world of numerical bliss.
> Many of us took the trip.

Well, there is that.  But it was also common to have
"black box" routines such that you didn't know what was
inside but called them and hoped for the correct answer.
(There are still many closed source libraries.)

NR routines were open.  You were supposed to look inside,
see how they work, and possibly modify them to do
what you really needed.

The book explains the theory somewhere between the beginner
and expert level.  Readable by someone with some experience,
but not up to what an expert would need.  Experts should

Also, not so bad for someone expert in one part of numerical
analysis, but needing a quick introduction to a different
part of the subject.

-- glen

Sat, 10 Sep 2011 02:49:15 GMT
Problem with qromb subroutine/nrerror

Quote:

>      SUBROUTINE NewTrapzd(a, b, Resu, N, yy)
>      !use nrtype
>      implicit none
>      !real(sp)::resu
>      REAL(Kind=8) :: A, B, X, YY, Del,func,Resu
>      INTEGER :: I, N
>      if(N==1) then
>      Resu = 0.5D0*(func(A,YY) + func(B,YY))
>      else
>      Resu = 0.5D0*(func(A,YY) + func(B,YY))
>      Del = (B-A)/DBLE(N-1) !h=(b-a)/(n-1)
>      DO I=1, N-2
>       X = A + I*Del  !x_i=x0+ih
>       Resu = Resu + Del*func(X,YY) !f(x_i)*h--> ...
>      END DO
>      end if
>      END SUBROUTINE NewTrapzd

I didn't look all the way through, but this looks wrong.
Assuming that func() is a constant, C, if N==1 then
resu=c.  If N is not 1 then Resu starts out C,
but then has more terms added to it.  An integration
routine should at least get the right answer for a
constant, and that would be a good first test.

-- glen

Sat, 10 Sep 2011 02:54:05 GMT
Problem with qromb subroutine/nrerror

Quote:

> (someone wrote)
>>> Is it NR or a lack of experience that encourages
>>> students to get in over their heads???

>> The environment into which NR was born had a lot to do
>> with both the erstwhile popularity and the risks associated
>> with NR. It arrived at a time when few had network access.
>> You had to go to the library, photocopy pages from journals
>> such as TOMS, and transcribe the Algol-like algorithm
>> listings into working software. NR sailed in, seductively
>> offering an easy path to a dream world of numerical bliss.
>> Many of us took the trip.

> Well, there is that.  But it was also common to have
> "black box" routines such that you didn't know what was
> inside but called them and hoped for the correct answer.
> (There are still many closed source libraries.)

> NR routines were open.  You were supposed to look inside,
> see how they work, and possibly modify them to do
> what you really needed.

Yes, and that encouraged tinkering and learning from doing so. However, this statement from the preface of the 2007 edition qualifies "open":

"..amount of code from this book for use in their own applications. If you personally
keyboard no more than 10 routines from this book into your computer, then we au-
thorize you (and only you) to use those routines (and only those routines) on that
single computer. You are not authorized to transfer or distribute the routines to any
..."

I wondered if the authors considered that people were going to put PDF versions of the book on the Web, from which one could "mouse" instead of "keyboard" routines.

Quote:
> The book explains the theory somewhere between the beginner
> and expert level.  Readable by someone with some experience,
> but not up to what an expert would need.  Experts should

> Also, not so bad for someone expert in one part of numerical
> analysis, but needing a quick introduction to a different
> part of the subject.

> -- glen

I agree with most of your comments. In my case, the book served only as a stepping-stone. The current edition comes only in C++, if I am not mistaken.

-- mecej4

Sat, 10 Sep 2011 04:03:49 GMT
Problem with qromb subroutine/nrerror

Quote:

>>    Dear all ;

>> I'm trying to compute an integral with Romberg's Method.
>> I used qromb subroutine  which calls polint and trapzd.
>> but I cannot get the answer correctly.

>> the program stops the running with nrerror .
>> me? I checked all parts but I didn't find the mistake.

> A little knowledge is a dangerous thing.

"A little learning is a dangerous thing
Drink deep, or taste not the Pierian spring"

Alexander Pope

;-)

Sat, 10 Sep 2011 06:23:25 GMT
Problem with qromb subroutine/nrerror

Quote:

> [OT?] Sometimes it is difficult to understand how otherwise competent
> science and engineering students have such trouble with Fortran. There
> are *some* obscure rules, but more often the students seem to trip on
> the fundamentals. Are there any good textbooks that cover Fortran,
> general programming principles and some of the numerical techniques
> that would be of use to those students?

It generally doesn't take much experimentation to figure out what the
language is doing.  The skill needed is a kind of general
problem-solving ability.  Is it learned or innate?  I've been amazed to
discover how unpractical in simple domestic or vehicle repairs and
maintenance some very bright physicists are.  Something about theory and
practice ...

Sat, 10 Sep 2011 06:31:11 GMT
Problem with qromb subroutine/nrerror

Quote:

>> A little knowledge is a dangerous thing.

>"A little learning is a dangerous thing
>Drink deep, or taste not the Pierian spring"

>Alexander Pope

>;-)

A professional pedant pontificates:

A little learning is a dang'rous thing;
Drink deep, or taste not the Pierian spring:
There shallow draughts intoxicate the brain,
And drinking largely sobers us again.

Inter alia, note the use of the colon: something that is sadly
out of fashion.

Regards,
Nick Maclaren.

Sat, 10 Sep 2011 05:44:47 GMT

 Page 1 of 2 [ 28 post ] Go to page: [1] [2]

Relevant Pages