Problem with qromb subroutine/nrerror
Author 
Message 
Fateme #1 / 28

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 . 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 = (BA)/DBLE(N1) !h=(ba)/(n1) DO I=1, N2 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.D6)yy=0.1D0 if (abs(x)<1.D6)then func=1/(2*yy) else if (abs(x1)<1.D6)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=K1 REAL(SP), PARAMETER :: EPS=1.0e6_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(jKM:j),s(jKM: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=xax ns=iminloc(abs(xxa)) !Find index ns of closest table entry. y=ya(ns) !This is the initial approximation to y. ns=ns1 do m=1,n1 !For each column of the tableau, den(1:nm)=ho(1:nm)ho(1+m:n) !we loop over the current c's and d's and up if(any(den(1:nm) == 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:nm)=(c(2:nm+1)d(1:nm))/den(1:nm) d(1:nm)=ho(1+m:n)*den(1:nm) !Here the c's and d's are updated. c(1:nm)=ho(1:nm)*den(1:nm) if (2*ns < nm) 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 tableauforking 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=ns1 end if y=y+dy end do END SUBROUTINE polint

Fri, 09 Sep 2011 19:28:09 GMT 


mecej #2 / 28

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

Fri, 09 Sep 2011 20:32:13 GMT 


e p chandle #3 / 28

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

Fri, 09 Sep 2011 22:53:00 GMT 


e p chandle #4 / 28

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


mecej #5 / 28

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 . >> > 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.
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 1D minimization code for a unidirectional minimization of an nD 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 > 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 Algollike 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 


e p chandle #6 / 28

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 . > >> > 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 1D minimization code for a unidirectional minimization of an nD 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 Algollike 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 


James Van Buskir #7 / 28

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 singleprecision REALs, like SP = KIND(1.0). Why then are you using singleprecision variables all over the place when your core function func is double precision? Also you define double precision constants such as 1.D6. It's much more flexible to define them as 1.0e6_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.5794487871554595D85, & 6.0134700243160014d154/),(/'x'/)); end

Sat, 10 Sep 2011 02:21:27 GMT 


glen herrmannsfeld #8 / 28

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

Sat, 10 Sep 2011 02:49:15 GMT 


glen herrmannsfeld #9 / 28

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 = (BA)/DBLE(N1) !h=(ba)/(n1) > DO I=1, N2 > 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 


mecej #10 / 28

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 Algollike 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 > 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 steppingstone. The current edition comes only in C++, if I am not mistaken.  mecej4

Sat, 10 Sep 2011 04:03:49 GMT 


Gib Bogl #11 / 28

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

Sat, 10 Sep 2011 06:23:25 GMT 


Gib Bogl #12 / 28

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


n.. #13 / 28

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] 
