Programming succesive trapezoidal integrator

[Telemachus]

Hi, I was trying to write a code in fortran that takes an array 'fx' of length 'm', and integrates it between a and b using the successive trapezoidal rule. Off course, the length of the array, m, has to be a power of 2, and the underlying grid uniform. So, I'm not sure on how to call the vector indices, in order that, if for example, n=2, in rsum it should be calling f(m/2). If n=4, it should call f(m/4) and f(3*m/4), when n=8, it should call for f(m/8), f(3*m/8), f(5*m/8), f(7*m/8), etc. In general $n=2^l$ with l some exponent, and it determines how many partitions will be done for the interval (a,b). More over, the program is running too slow, in my opinion should go way much much faster. It takes several seconds if for example I choose m=n=50.

So, this is the subroutine I have written. The result is written in the variable 's', and iterates 'n' times on the successive trapezoidal rule:

      SUBROUTINE trapzoidal(fx,m,a,b,s,n)
  
  
      integer m
      real*8 fx
      dimension fx(m)

      REAL*8 a,b,s,h
      integer n,j,k

      REAL*8 T(100)

      j=1
      h=(b-a)
      T(1)=h*(fx(1)+fx(m))/2.d0
      do i=1,n
        j=2*j
        rsum=0.d0
        h=h/2.d0
        do k=1,j,2
          rsum=rsum+fx(k*m/j)
        enddo
        T(i+1)=T(i)/2.d0+h*rsum
      enddo

      s=T(n+1)      return
      END

 

[DrClaude]

I don't get it. If fx is an array, meaning that you have the value of the function at discrete points, why can't you simply loop over the points? And why keep track of the progress of the integral by storing all intermediate values of T is you are not going to use them?

 

[Mark44]

Hi, I was trying to write a code in fortran that takes an array 'fx' of length 'm', and integrates it between a and b using the successive trapezoidal rule. Off course, the length of the array, m, has to be even, and the underlying grid uniform. So, I'm not sure on how to call the vector indices, in order that, if for example, n=2, in rsum it should be calling f(m/2). If n=4, it should call f(m/4) and f(3*m/4), when n=8, it should call for f(m/8), f(3*m/8), f(5*m/8), f(7*m/8), etc. In general $n=2^l$ with l some exponent, and it determines how many partitions will be done for the interval (a,b). More over, the program is running too slow, in my opinion should go way much much faster. It takes several seconds if for example I choose m=n=50.

So, this is the subroutine I have written. The result is written in the variable 's', and iterates 'n' times on the successive trapezoidal rule:

      SUBROUTINE trapzoidal(fx,m,a,b,s,n)
 
 
      integer m
      real*8 fx
      dimension fx(m)

      REAL*8 a,b,s,h
      integer n,j,k

      REAL*8 T(100)

      j=1
      h=(b-a)
      T(1)=h*(fx(1)+fx(m))/2.d0
      do i=1,n
        j=2*j
        rsum=0.d0
        h=h/2.d0
        do k=1,j,2
          rsum=rsum+fx(k*m/j)
        enddo
        T(i+1)=T(i)/2.d0+h*rsum
      enddo

      s=T(n+1)      return
      END

The first problem I see is that you probably aren't declaring your array parameter correctly. I would replace these two lines

real*8 fx
dimension fx(m)

with this line

real*8 fx(*)

Possibly this is what is slowing down your code, as I believe what might be happening is that you are passing what the subroutine believes is a double-precision real when in fact the actual argument is the address of an array.

Also, I don't understand what you are doing in your code. h would typically be the size of one subinterval, but in your code it is just the length of the entire interval. A subinterval in your grid would be (b - a)/n.

Try to see where you code is going wrong by choosing an interval divided into a small number (say 4) subintervals, and hand-simulating what your code is doing. Compare what the code does to what it should be doing.

It would also be very helpful to add comments...

 

[Telemachus]

Hi, thank you for your post. The idea is to use this to apply romberg integration later. I want to integrate an array that comes from the solution of a differential equation, that is the reason why I am using an array. However, you are right about the storage, I think it is unnecessary. But I need to loop this way in order to do the succesive trapezoidal rule, which starts with a rough $\int_a^bf(x)dx\sim \frac{1}{2}(f(a)-f(b))$ and then subdivides the interval, and applies the extended trapezoidal rule successively. In this way, then one can use these results to extrapolate the solution with very good accuracy. That is the reason why I'm trying to work it this way.

 

[DrClaude]

The first problem I see is that you probably aren't declaring your array parameter correctly. I would replace these two lines

real*8 fx
dimension fx(m)

with this line

real*8 fx(*)

Possibly this is what is slowing down your code, as I believe what might be happening is that you are passing what the subroutine believes is a double-precision real when in fact the actual argument is the address of an array.

I don't see anything wrong with that. It is perfectly standard Fortran 77. And Fortran is always pass by reference.

 

[Telemachus]

The first problem I see is that you probably aren't declaring your array parameter correctly. I would replace these two lines

real*8 fx
dimension fx(m)

with this line

real*8 fx(*)

Possibly this is what is slowing down your code, as I believe what might be happening is that you are passing what the subroutine believes is a double-precision real when in fact the actual argument is the address of an array.

Hi Mark, thank you for your post. It is declared that way because the array fx is allocated in order to vary m without recompilation. However, I will check if there is a problem with it.

Also, I don't understand what you are doing in your code. h would typically be the size of one subinterval, but in your code it is just the length of the entire interval. A subinterval in your grid would be (b - a)/n.

Try to see where you code is going wrong by choosing an interval divided into a small number (say 4) subintervals, and hand-simulating what your code is doing. Compare what the code does to what it should be doing.

It would also be very helpful to add comments...

Here there is an explanation for the successive trapezoidal rule: http://mathfaculty.fullerton.edu/mathews/n2003/rombergquad/RombergProof.pdf

 

[DrClaude]

Hi, thank you for your post. The idea is to use this to apply romberg integration later. I want to integrate an array that comes from the solution of a differential equation, that is the reason why I am using an array. However, you are right about the storage, I think it is unnecessary. But I need to loop this way in order to do the succesive trapezoidal rule, which starts with a rough $\int_a^bf(x)dx\sim \frac{1}{2}(f(a)-f(b))$ and then subdivides the interval, and applies the extended trapezoidal rule successively. In this way, then one can use these results to extrapolate the solution with very good accuracy. That is the reason why I'm trying to work it this way.

Ok, now I get it.

To correct what you wrote initially, this means that m has to be of size $2^l + 1$. m = 50 will definitely not work!

Otherwise, I think that what you wrote should work.

 

[Telemachus]

Ok, now I get it.

To correct what you wrote initially, this means that m has to be of size $2^l + 1$. m = 50 will definitely not work!

Otherwise, I think that what you wrote should work.

Yes, I realized of that later. That was a silly mistake. Do you think the code is right? I've tried to use it with a subroutine that does romberg integration, and it doesn't work as it should.

I am using this array to test it:

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!... Create array to be evaluated. This is done in order to simulate
!... a realistic numerical situation where f is not know anallytically

        subroutine createarray(fx,m)
        implicit none
        real*8 y1,y2
        integer mxpts
        common/bkinput/ y1,y2,mxpts

        integer m
        real*8 fx
        dimension fx(m)

        integer i
        real*8 dx,x
  
        dx=(y2-y1)/mxpts
        do i=1,mxpts+1
         x=y1+(i-1)*dx
         fx(i)=4.d0*dsqrt(1.d0-x**2)
        enddo

        return
        end

Which integrated between 0 and 1 should equal $\pi$.

 

[Mark44]

I don't see anything wrong with that. It is perfectly standard Fortran 77.

I was going by this example, from a class at Stanford using Fortran 77 - http://web.stanford.edu/class/me200c/tutorial_77/12_arrays2.html.

      subroutine saxpy (n, alpha, x, y)
      integer n
      real alpha, x(*), y(*)

 

[Telemachus]

Ok, I'll post the entire code, if you want to try it. Romberg integration should work much better than this, I think.

        program integrator

!... This program integrates a function f(x) by means of Romberg integration.

!... Version 0.0.1

        implicit none

        REAL*8, allocatable :: fx(:)

!... Input parameters        
        real*8 y1,y2
        integer mxpts
        common/bkinput/ y1,y2,mxpts

        integer m
        integer ierr,stat
!... Input file
        open (unit=10,file='romberg.inp',status='old')

!... call input      
        call readinput(m)        allocate(fx(m),stat=ierr)
        if(ierr.ne.0)stop 'allocation fails - fx'

!... Create array of data
        call createarray(fx,m)

!... Romberg integration
        call integrateRomb(fx,m)

!... Close files
     deallocate(fx)

        close (unit=10)
        
        end
      
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        subroutine readinput(m)
        implicit none    
      
!... Input parameters        
        real*8 y1,y2
        integer mxpts
        common/bkinput/ y1,y2,mxpts

        integer m

        real*8 rzero, one, two, five
        data rzero,one,two,five/0.0d0,1.0d0,2.0d0,5.0d0/
      
      
        namelist /griddata/mxpts,y1,y2
              
        mxpts=0
        y1=rzero
        y2=rzero
      
!... read griddata

        read(10,griddata)
        m=mxpts+1
        write(*,111) y2,y1,m
111     format(/5x,'y2=',f9.3,5x,'y1=',f9.3,5x,'m=',i8,5x)
        return
        end
      

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!... Create array to be evaluated. This is done in order to simulate
!... a realistic numerical situation where f is not know anallytically

        subroutine createarray(fx,m)
        implicit none  
        real*8 y1,y2  
        integer mxpts
        common/bkinput/ y1,y2,mxpts

        integer m
        real*8 fx
        dimension fx(m)

        integer i
        real*8 dx,x
    
        dx=(y2-y1)/mxpts
        do i=1,mxpts+1
         x=y1+(i-1)*dx
         fx(i)=4.d0*dsqrt(1.d0-x**2)
        enddo  

        return
        end
      

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE trapzd(fx,m,a,b,s,n) 
    
    
      integer m
   real*8 fx
      dimension fx(m)
      REAL*8 a,b,s,h
      integer n,j,k

      REAL*8 T(200)

      j=1
      h=(b-a)
      T(1)=h*(fx(1)+fx(m))/2.d0
      do i=1,n
        j=2*j
        rsum=0.d0
        h=h/2.d0
        do k=1,j-1,2
          rsum=rsum+fx(k*m/j)
        enddo
        T(i+1)=T(i)/2.d0+h*rsum
      enddo
      s=T(n+1)

      return 
      END 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      SUBROUTINE polint(xa,ya,n,x,y,dy) 

!... From numerical recipes in Fortran77

      INTEGER n,NMAX 
      REAL*8 dy,x,y,xa(n),ya(n) 
      PARAMETER (NMAX=10) 
      INTEGER i,m,ns 
      REAL*8 den,dif,dift,ho,hp,w,c(NMAX),d(NMAX) 
      ns=1 
      dif=abs(x-xa(1)) 
      do 32 i=1,n 
        dift=abs(x-xa(i)) 
        if (dift.lt.dif) then 
          ns=i 
          dif=dift 
        endif 
        c(i)=ya(i) 
        d(i)=ya(i) 
32    continue 
      y=ya(ns) 
      ns=ns-1 
      do 35 m=1,n-1 
        do 34 i=1,n-m 
          ho=xa(i)-x 
          hp=xa(i+m)-x 
          w=c(i+1)-d(i) 
          den=ho-hp 
          if(den.eq.0.)stop 'failure in polint' 
          den=w/den 
          d(i)=hp*den 
          c(i)=ho*den 
34      continue 
        if (2*ns.lt.n-m)then 
          dy=c(ns+1) 
        else 
          dy=d(ns) 
          ns=ns-1 
        endif 
        y=y+dy 
35    continue 
      return 
      END 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE qromb(fx,m,a,b,ss) 

!... From numerical recipes in Fortran77

      integer m
   real*8 fx
      dimension fx(m)

      real*8 y1,y2
      integer mxpts
      common/bkinput/ y1,y2,mxpts

      INTEGER JMAX,JMAXP,K,KM 
      REAL*8 a,b,ss,EPS,dss

      PARAMETER (EPS=1.e-6, K=2, KM=K-1) 
CU    USES polint,trapzd 
      INTEGER j 

      REAL*8, allocatable :: h(:),s(:)

      JMAX=mxpts
      JMAXP=JMAX+1

      ALLOCATE( h(JMAXP), s(JMAXP) )

      h(1)=1.0d0 
      do 33 j=1,JMAX
        call trapzd(fx,m,a,b,s(j),j) 
        if (j.ge.K) then 
          call polint(h(j-KM),s(j-KM),K,0.0d0,ss,dss) 
        endif 
        s(j+1)=s(j) 
        h(j+1)=0.25d0*h(j) 
33    continue 

   deallocate(h,s)

      END 

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        subroutine integrateRomb(fx,m)

!... This routine performs the numerical integration of the function
!... f(x) using Romberg integration.

        implicit none 

        real*8 y1,y2
        integer mxpts
        common/bkinput/ y1,y2,mxpts

        integer m
     real*8 fx
        dimension fx(m)        integer i
        real*8 rint,intaux,ss

        call qromb(fx,m,y1,y2,ss)        print*,'Integral value by Romberg integration is=',ss

        print*,'Romberg-Exact=',
     +  3.14159265358979323846264338327950288419716d0-ss
!     +  7.954926521012845274513219665329394328161342771816638573400-ss 

        return
        end

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This is the input file, which should be called "romberg.inp"

!... Input file for romberg.f

!... griddata
!      y1: beginning of the integration interval
!      y2: endpoint for the integration interval
!      mang: number of points for the Gauss-Legendre quadrature

 &griddata mxpts=16 y1=0.d0 y2=1.d0 &end

 

[Mark44]

Which integrated between 0 and 1 should equal $\pi$.

Aren't you finding the area between the curve $y = \sqrt{1 - x^2}$ and the x-axis, for $0 \le x \le 1$? In that case, the area would be $\frac \pi 4$, not $\pi$.

 

[Telemachus]

No, it's 4 times that integrand.

 

[DrClaude]

What value of n are you using?

I was going by this example, from a class at Stanford using Fortran 77 - http://web.stanford.edu/class/me200c/tutorial_77/12_arrays2.html.

      subroutine saxpy (n, alpha, x, y)
      integer n
      real alpha, x(*), y(*)

These are called assumed-sized arrays. Both approaches are valid. Passing the size of the array as an argument and declaring its exact dimension allows for out-of-bound checks.

 

[DrClaude]

The value of n used is incorrect. For $m = 2^l + 1$, n must be at most l.

In trapeze, replacing fx(k*m/j) by fx(k*m/j+1) appears to work.

 

[Telemachus]

The value of n used is incorrect. For $m = 2^l + 1$, n must be at most l.

In trapeze, replacing fx(k*m/j) by fx(k*m/j+1) appears to work.

Right! thank you verymuch.

 

[Tom.G]

A minor quibble, the compiler may already optimize these.

SUBROUTINE trapzd (
!
!
!
!
do i=1,n
j=2*j
rsum=0.d0
h=h/2.d0
do k=1,,j-1,2
rsum=rsum+fx(k*m/j)
enddo
T(i+1)=T(i)/2.d0+h*rsum
enddo
s=T(n+1)

return
END

Recode and move these out of their For Loops:
j=2*j → j=j+j
rsum=rsum+fx(k*m/j) → move outside the loop variable=m/j, and in the loop change to k*variable

 

[Telemachus]

Alright, thanks. However, I think I won't use this integration, it isn't as good as I thought, and I have other integration methods that achieve similar accuracy with much less effort.

I think that before the code was running slow because I have defined wrong the proper values for the iteration number as DrClaude noted.

I've been testing many integration methods, and I wanted to know, how long can the step length be for a quadrature to work? I don't know what the issue is, but if I take, for example only two points, the thing works fine, with some degradation in the accuracy of course. With two points I have a step h=1/2 in this case. But if I increase the interval to, let's say b=10, it doesn't work, even if I take a lot of points, and at this point I'm not seeing why that is. It gives me non arithmetic number, and I don't know why, probably some bug in the code.

 

[DrClaude]

I've been testing many integration methods, and I wanted to know, how long can the step length be for a quadrature to work? I don't know what the issue is, but if I take, for example only two points, the thing works fine, with some degradation in the accuracy of course. With two points I have a step h=1/2 in this case. But if I increase the interval to, let's say b=10, it doesn't work, even if I take a lot of points, and at this point I'm not seeing why that is. It gives me non arithmetic number, and I don't know why, probably some bug in the code.

What function are you testing it on? I suggest you use polynomials, and progressively increase the order. The trapezoidal rule is exact for first order polynomials, so you can check that and then check how much of an error you get when you go to a second-order polynomial.

 

[Telemachus]

I was using the same integral than before: $\int_a^b \sqrt{1-x^2}dx$. First I was using a=0 and b=1. If I decrease b, it's ok. But if I increase it, I get NaN.

 

[DrClaude]

I was using the same integral than before: $\int_a^b \sqrt{1-x^2}dx$. First I was using a=0 and b=1. If I decrease b, it's ok. But if I increase it, I get NaN.

What is the value of $\sqrt{1-x^2}$ when $x>1$?

 

[Telemachus]

Oh! I'm such an idiot! sorry.