Problem with qromb subroutine/nrerror
Fatemeh
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 .
Is there anyone who has any experience about this and can guide me?
I checked all parts but I didn't find the mistake.
result which can be printed :
u= 0.00000000000000 answer= -0.900316316157106
u=0.100000000000000 answer= -7.64440298080444
nrerror: qromb:too many steps
STOP program terminated by nrerror
program test2
IMPLICIT NONE
REAL(Kind=8)::frac,u,t,answer,qromb
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
answer=-4/PI*sin(PI/4)
else
answer=qromb(0.D0,1.D0,frac)
print*,u,answer,-4*t*answer
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)
s(j)=ss !added by Fatemeh
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
mecej4
> 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 .
> Is there anyone who has any experience about this and can 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.
>
>
>
> program test2
> IMPLICIT NONE
> REAL(Kind=8)::frac,u,t,answer,qromb
> 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
> answer=-4/PI*sin(PI/4)
> else
> answer=qromb(0.D0,1.D0,frac)
> print*,u,answer,-4*t*answer
>
> 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
e p chandler
> Fatemeh wrote:
> > 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 .
> > Is there anyone who has any experience about this and can 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.
> > answer=qromb(0.D0,1.D0,frac)
> > print*,u,answer,-4*t*answer
> 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.
Is it NR or a lack of experience that encourages students to get in
over their heads???
--- e
e p chandler
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
mecej4
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.
> 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 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.
-- mecej4
e p chandler
> e p chandler wrote:
> > On Mar 23, 8:32 am, mecej4 <mec...@operamail.invalid> wrote:
> >> Fatemeh wrote:
> >> > 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 .
> >> > Is there anyone who has any experience about this and can
> >> > 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
James Van Buskirk
news:f7c64083-e7d0-494a...@d19g2000yqb.googlegroups.com...
> 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
glen herrmannsfeldt
(someone wrote)
>> 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
read the journals instead.
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
glen herrmannsfeldt
> 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
mecej4
> mecej4 <mec...@operamail.invalid> wrote:
> (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.
> 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
> read the journals instead.
>
> 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
Gib Bogle
> Fatemeh wrote:
>
>> 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 .
>> Is there anyone who has any experience about this and can guide
>> 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
;-)
Gib Bogle
> [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 ...
nm...@cam.ac.uk
Gib Bogle <bo...@ihug.too.much.spam.co.nz> wrote:
>mecej4 wrote:
>>
>> 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.
glen herrmannsfeldt
> 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.
Fortran does have some places where things are less
obvious than they should be. Many are the way they
are for historical reasons.
> 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 ...
Google for "pauli effect", most likely the first link.
-- glen
mecej4
Indeed; and compartmentalization of knowledge.
"Q: How many programmers does it take to change a light bulb? A: None; that's a hardware problem."
-- mecej4
e p chandler
> 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.
> answer=qromb(0.D0,1.D0,frac)
[snip]
> real(kind=8) function func(x,yy)
[snip]
Here's a short demonstration program.
The integration technique used is a two point Gaussian quadrature
method. Unlike Romberg/Trapezoid, it is not adaptive. N is the number
of intervals. This may be varied. Gaussian quadrature has a few
problems but it is easy to code and remarkably accurate for low values
of N for some integrands. It is difficult to estimate the errror. It
is difficult to determine what value of N to choose (except by
experiment).
----- start text ----
module funcs
implicit none
integer p
contains
function f(x)
real x,f
f=x**p
end function f
function g(x)
real x,g
g=1/x
end function g
end module funcs
module integral
implicit none
contains
subroutine gauss(f,a,b,s,n)
implicit none
real f,a,b,s
integer n
real x,dx,hx,gx
integer i
dx=(b-a)/n
hx=dx/2
gx=hx/sqrt(3.0)
s=0
do i=0,n-1
x=a+(i+.5)*dx
s=s+f(x-gx)+f(x+gx)
end do
s=s*hx
end subroutine gauss
end module integral
program test
use funcs
use integral
implicit none
real q
integer n
print *,'p,integral [0,1] of x**p'
do p=0,3 ! module parameter
call gauss(f,0.0,1.0,q,16)
print *,p,q
end do
print *
print *,'n,integral [1,2] of 1/x'
do n=1,10
call gauss(g,1.0,2.0,q,n)
print *,n,q
end do
end program test
----- end text ----
Sample output:
C:\Users\epc\temp>g95 x.f90
C:\Users\epc\temp>a
p,integral [0,1] of x**p
0 1.
1 0.5
2 0.3333333
3 0.25
n,integral [1,2] of 1/x
1 0.6923077
2 0.6930766
3 0.69313216
4 0.6931423
5 0.69314516
6 0.69314617
7 0.69314665
8 0.6931468
9 0.693147
10 0.69314706
C:\Users\epc\temp>
---- end output ----
Even if you don't like Gaussian integration, this program shows how to
put an integration routine in a module, how to pass functions as
arguments to that routine and how to pass additional arguments to the
integrand without changing the integration module's code.
-- e
Gib Bogle
Yes. Maybe while I was making excellent use of my time as a young man, learning how to take my
motorcycle apart and put it together again, they were wasting theirs in understanding theoretical
physics.
Q: How many psychotherapists does it take to change a light bulb? A: Only one, but the bulb has to
really want to change.
robin
> 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.
????
TOMS routines are mostly already working software.
mecej4
I wrote about the early 1980-s, at which time hardly any of us had network access.
-- mecej4
nm...@cam.ac.uk
mecej4 <mec...@operamail.invalid> wrote:
>robin wrote:
>
>> TOMS routines are mostly already working software.
>
>I wrote about the early 1980-s, at which time hardly any of us had
> network access.
Which backward country was that? :-)
You are correct, in that TOMS was available primarily from paper at
that date - so were most of the other sources, such as Applied
Statistics Algorithms, too.
Regards,
Nick Maclaren.
mecej4
> In article <gqafu3$a42$1...@news.motzarella.org>,
> mecej4 <mec...@operamail.invalid> wrote:
>>robin wrote:
>>
>>> TOMS routines are mostly already working software.
>>
>>I wrote about the early 1980-s, at which time hardly any of us
>>had
>> network access.
>
> Which backward country was that? :-)
The same country in which everyone was waiting their turn for these men to set them up:
http://upload.wikimedia.org/wikipedia/commons/thumb/3/35/Phoc96v1.jpg/200px-Phoc96v1.jpg
> You are correct, in that TOMS was available primarily from paper
> at that date - so were most of the other sources, such as
> Applied Statistics Algorithms, too.
>
>
> Regards,
> Nick Maclaren.
--
-- mecej4
Jugoslav Dujic
> e p chandler wrote:
>
> 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?".
She, to be quite precise.
While that does not excuse Fatemeh's blunders, it might remind
chiefly male population of clf to be more gentle towards fellow
posters, especially if they're apparent newcomers.
--
Jugoslav
www.xeffort.com
Please reply to the newsgroup.
You can find my real e-mail on my home page above.