HUR1EY
HUR1EY

Reputation: 1

writing to an array in a loop

the program writes 20 and 200 times the last value to the array. How can you make it so that x is in the range of -2.2 with a step of 0.1 and 0.01, and at the same time the values are correctly written to the array?

do j=0,200
            do x = -2, 2, 0.01
             Yb(j) = exp(-x/5)
            end do
          end do
module integration_methods
    implicit none
    contains
      real function trapezoid_rule(x, n)
        real :: f
        real ::  a, b
        integer, intent(in) :: n
        real, dimension(n), intent(inout) :: x
        integer :: k
        real :: s
        s = 0
        do k=1, n-1
          s = s + ((x(k)+x(k+1))/2*n)
        end do  
        trapezoid_rule = s/N
      end function trapezoid_rule
    end module integration_methods
    
    program integration
      use integration_methods
      implicit none
      real, external :: s
      integer :: j, i, m = 200, n = 20
      Real, Allocatable, Dimension(:) :: Ya
      Real, Allocatable, Dimension(:) :: Yb
      real :: x, integral
      ALLOCATE(Ya(1:n),Yb(1:m))
      do i=0, 20
      do x = -2, 2, 0.1
         Ya(i) = exp(-x/5)
      end do
    end do
        do j=0,200
        do x = -2, 2, 0.01
         Yb(j) = exp(-x/5)
        end do
      end do
      integral=trapezoid_rule(Ya,  20)
      write (*,*) 'Trapezoid rule = ', integral
      integral=trapezoid_rule(Yb, 200)
      write (*,*) 'Trapezoid rule = ', integral
    deallocate(Ya,Yb)
    end program integration 

Upvotes: 0

Views: 273

Answers (3)

Holmz
Holmz

Reputation: 724

Another approach:

Replace this:

do j=0,200
  do x = -2, 2, 0.01
    Yb(j) = exp(-x/5)
  end do
end do

With:

Real :: MinX = -2.0D0
REAL :: MaxX = 2.0D0
REAL :: Step = 0.01D0
INTEGER :: LHS = MinX/Step
INTEGER :: RHS = MaxX/Step
REAL ALLOCATABLE, DIMENSION(:) :: Yb

ALLOCATE(YB(LHS:RHS))

do j= LHS,RHS
  x = i/Step
  Yb(j) = exp(-x/5)
end do

That assumes that the module is modified, like by replacing this:

do k=1, n-1
  s = s + ((x(k)+x(k+1))/2*n)
end do  

With this:

do k= LBOUND(X, DIM=1), UBOUND(X, DIM-1)-1
   s = s + ((x(k)+x(k+1))/2*n)
end do  

And probably some other passing of the “step” to the module, and making sure that the INTENT is allowing the LBOUND and UBOUND to be retrieved correctly.

Upvotes: 0

John Alexiou
John Alexiou

Reputation: 29244

You don't need two nested loops, just one

  do j=0,200
     x = -2.0 + 4.0 * real(j)/200
     Yb(j) = exp(-x/5)
  end do

Upvotes: 0

Steve Lionel
Steve Lionel

Reputation: 7267

First, using a REAL variable as a DO loop index is no longer standard Fortran. The reason is demonstrated by your program in that .1 and 01 are not exactly representable in binary floating point. The semantics of DO loops are that the loop count is computed as (start-end)/increment and with floating point numbers this can result in a higher or lower number of iterations, depending on the values chosen.

But the real problem is that your do i loop starts at zero, but the arrays Ya and Yb are 1-origin, leading to data corruption.

Upvotes: 3

Related Questions