Data Interpolation For use with Intergration

[VoidSerpent]

I am very new to FORTRAN and have pieced together parts of a code to help process my data outside commercial software. The partial code I will provide below has a very specific purpose. Given data that can be assumed as an M x N matrix, The first row is all string names for the column beneath it and the second row is the units of the data. The only way i can see to read in the data with FORTRAN is as an allocatable array which i conceive to be an M x N matrix once declared. From this matrix i need to construct arrays from the columns of the data. This way i have the proper data with its correct representation i.e. time, distance, temperature and etc... So there should be N arrays constructed from the data matrix. With those arrays i need to select four for interpolation, one will be held as the X value for use against the other three. The interpolation will give a function for each Y value array, that should five me three functions that i then need to multiply together to obtain a single function for my integration code. The problems i have are the construction of my arrays and the interpolation, the interpolation is tricky because i don't know if i should do a spline function for my exponential data, even though they seem to be the only references i can find for interpolation. Taking the log of the Y value arrays would produce a linear system but i cant find linear interpolation code. So this is all i have.

program datain
 real, dimension(:), allocatable :: x, oldx
 real a
 integer io, n
 character(30) :: filename
 allocate( x(0) ) ! size zero to start with?
 n = 0
write(*,*) 'Enter input file name: '
read(*,*) filename
 open( 1, file = filename )
 do
 read(1, *, iostat = io) a
 if (io < 0 ) exit
 n = n + 1
 allocate( oldx( size(x) ) )
 oldx = x 
 deallocate( x )
 allocate( x(n) )
 x = oldx
 x(n) = a ! The nth value of x
 deallocate( oldx )
 end do
end

! Decompose matrix array x into column arrays [1,n]

 

! Interpolate Column arrays, Constant X value will be array ALT with 3 other arrays to compose 3 functions
!Those 3 functions are then multipled to make a single function used during intergration


program integration
 implicit none
 external f1 ! define function to be integrated as external from interpolation
 real :: integral
 real topbound
 real bottombound
 ! integrating the function f1
 write(*,*) 'Enter Top Bound: '
read(*,*) topbound
write(*,*) 'Enter Bottom Bound: '
read(*,*) bottombound
 integral=trapezoid_rule(f1, bottombound, topbound, 100)
 write (*,*) 'Trapezoid rule = ', integral
 integral=simpson_rule(f1, bottombound, topbound, 100)
 write (*,*) "Simpson's rule = ", integral

contains
! Functions
real function trapezoid_rule(f, a, b, n)
  ! Approximate the definite integral of f from a to b by the
  ! composite trapezoidal rule, using N subintervals
  real :: f
  ! [a, b] : the interval of integration
  real, intent(in) ::  a, b
  ! n : number of subintervals used
  integer, intent(in) :: n
  ! temporary variables
  integer :: k
  real :: s
  s = 0
  do k=1, n-1
    s = s + f(a + (b-a)*k/n)
  end do  
  trapezoid_rule = (b-a) * (f(a)/2 + f(b)/2 + s) / n
end function trapezoid_rule

real function simpson_rule(f, a, b, n)
  ! Approximate the definite integral of f from a to b by the
  ! composite Simpson's rule, using N subintervals
  real :: f
  ! [a, b] : the interval of integration difined above
  real, intent(in) ::  a, b
  ! n : number of subintervals used
  integer, intent(in) :: n
  ! temporary variables
  integer :: k
  real :: s
  s = 0
  do k=1, n-1
    s = s + 2**(mod(k,2)+1) * f(a + (b-a)*k/n)
  end do  
  simpson_rule = (b-a) * (f(a) + f(b) + s) / (3*n)
end function simpson_rule

end program integration 

! define here functions to be integrated 
real function f1(x)
  implicit none
  real, intent(in) :: x
  f1 = fintpol !Function defined during interpolation
end function f1

I would appreciate guidance whether in links or example code for obtain my arrays and a usable interpolation method. I hope i was clear on how this code should work when it is done.

 

[mikrom]

Hi VoidSerpent,

Sorry, but I really don't understand what is your problem. In the code you posted are 2 programs merged together. One of them is what I posted back in 2009 here:

http://www.tek-tips.com/viewthread.cfm?qid=1570169

Please post only the relevant code which is problematic.

but i cant find linear interpolation code

There are some examples in this forum - look at the code of measurement_interpolation.f95 here

http://www.tek-tips.com/viewthread.cfm?qid=1542335

or lagrange.f95 here

http://www.tek-tips.com/viewthread.cfm?qid=1694168

 

[gummibaer]

I am very new to FORTRAN

I'd say, this issue with the interpolation formula is one of your minor problems. First you will have top make this thing run.

You should make yourself familiar on how to build a program with one main and a few subroutines or functions and make it run. This is depending on the programming environment.
Then you will have to think about on how to organise your data between your routines, I would recommend you would read about what a module is.
Then you should think about on how to load your data to memory. From what I understand of your post, you should find out about derived data types.

In all, when you start out being new to Fortran, I think this is not a very good example to learn, because you would need a few advanced features to do it efficiently.


Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[FJacq]

  n = n + 1
 allocate( oldx( size(x) ) )
 oldx = x 
 deallocate( x )
 allocate( x(n) )
 x = oldx
 x(n) = a ! The nth value of x
 deallocate( oldx )

This instruction x=oldx is wrong because the two arrays have not the same size. Write instead x(1:n-1)=oldx

François Jacq