Doubt about interpolation and integration with derivatives

[R.Silveira]

Hello , I'm doing a calculate with Fortran which involves an integral of data set and a numerical derivative.

So first i calculate de derivative of the data set using a three-point numerical differentiation

 subroutine derivation(np,Pdata,Mdata,dMdata)
      implicit none
      INTEGER :: i,np    
      REAL(8),DIMENSION (np) :: Mdata,dMdata,Pdata
C      
      open(unit=10,file="M(p).out")  
      DO i= 1, np
         READ(10,*)Pdata(i), Mdata(i)
      ENDDO   
      DO i = 2, np-1 
         dMdata(i) = (Mdata(i+1) - Mdata(i-1)) 
     $             / (Pdata(i+1) - Pdata(i-1))  
      ENDDO     
      CLOSE(10)
      return 
      end

i do a linear interpolation

      subroutine intlin1D(np,p,M,dM)
      implicit none
      integer    :: np,i,j
      double precision :: pdata(np),Mdata(np),dMdata(np),p,M,dM      
C      
      CALL derivation(np,pdata,Mdata,dMdata)     
      do  i = 1, np-1
         if((pdata(i).le.p).and.(p.lt.pdata(i+1))) then
            M = Mdata(i)+(Mdata(i+1)-Mdata(i))/
     .        (pdata(i+1)-pdata(i))*(p-pdata(i)) 
            dM = dMdata(i)+(dMdata(i+1)-dMdata(i))/
     .        (pdata(i+1)-pdata(i))*(p-pdata(i))     
         endif
      end do 
      if(p2data(np).le.p) then
         M = Mdata(np-1)
         dM = dMdata(np-1)       
      elseif(pdata(1).ge.p2)then
         M = Mdata(1)+(Mdata(2)-Mdata(1))/
     .        (pdata(2)-pdata(1))*(p-pdata(1))  
         dM = dMdata(1)+(dMdata(2)-dMdata(1))/
     .        (pdata(2)-pdata(1))*(p-pdata(1))     
      endif      
      return
      end

so a integrate a function which evolves M and dM (derivative of M) with Cuba library.

My doubts are, is there any problem in interpolating the data and the numeric derivative together? (because i calculate the same function with Numerical Integration - Gaussian-Legendre Quadrature too and the result in de end is different)

is it possible to calculate the numerical derivative of M that comes out of the interpolation instead of making the numerical derivative of the data, thus avoiding having to interpolate it?



this is de data i'm using , first column is pdata and second is Mdata

 4.0000000000000002E-004  0.48155697366644118     
   4.6243651255342965E-004  0.48153136433247140     
   5.3461882035644613E-004  0.48150204844084604     
   6.1806815707765898E-004  0.48146842554432528     
   7.1454320769829508E-004  0.48142994965416580     
   8.2607717259185133E-004  0.48138586131403838     
   9.5502061698343342E-004  0.48133520837718807     
   1.1040910088361102E-003  0.48127687706687472     
   1.2764299891694217E-003  0.48120954609567218     
   1.4756695817752976E-003  0.48113186058869717     
   1.7060087376933671E-003  0.48104237973500269     
   1.9723018276114985E-003  0.48093917687271143     
   2.2801609471585450E-003  0.48082006854785625     
   2.6360741911613077E-003  0.48068269063131219     
   3.0475423894818407E-003  0.48052435714576774     
   3.5232371861268251E-003  0.48034167202798700     
   4.0731837931276176E-003  0.48013105261563060     
   4.7089722707077154E-003  0.47988913299446051     
   5.4440017864422069E-003  0.47960896083928828     
   6.2937630011424316E-003  0.47928607847267535     
   7.2761645327152954E-003  0.47891358380159188     
   8.4119103781845408E-003  0.47848457742783890     
   9.7249362479991455E-003  0.47798976391773773     
   1.1242914008322915E-002  0.47742051560898680     
   1.2997834862367379E-002  0.47676024745597523     
   1.5026683561246378E-002  0.47600722992710653     
   1.7372217853266828E-002  0.47514024701332058     
   2.0083869598457850E-002  0.47414015578615443     
   2.3218786539221749E-002  0.47299164594570997     
   2.6843036682300574E-002  0.47167473102402246     
   3.1033000674267154E-002  0.47016216858316834     
   3.5876981514690841E-002  0.46843271407334797     
   4.1477065531493731E-002  0.46645324962418350     
   4.7951273838335098E-002  0.46418328068246484     
   5.5436049615735492E-002  0.46161726687750571     
   6.4089133635099022E-002  0.45867789287072369     
   7.4092888626964837E-002  0.45534071782346125     
   8.5658142554158115E-002  0.45155739426560604     
   9.9028631786373611E-002  0.44725238147602842     
  0.11448613781557095       0.44240351911969994     
  0.13235642576785989       0.43692163355841096     
  0.15301610986531489       0.43077689646306305     
  0.17690059052652182       0.42386144294190370     
  0.20451323037931773       0.41616042194768510     
  0.23643596256911972       0.40753306413305990     
  0.27334155493169215       0.39797881761856313     
  0.31600778849635835       0.38742670779277022     
  0.36533384912994471       0.37583338813469724     
  0.42235927777343146       0.36313509225157431     
  0.48828587864532691       0.34935823671631466     
  0.56450304712458088       0.33447802423824147     
  0.65261705109518930       0.31856335840222438     
  0.75448488285340631       0.30163733534600001     
  0.87225339500253274       0.28388199917252699     
   1.0084045451196499       0.26541457505731875     
   1.1658077027203950       0.24634088304621429     
   1.3477801208848612       0.22710169947721609     
   1.5581568469770919       0.20763767530329702     
   1.8013715458183344       0.18852868803826955     
   2.0825499386530293       0.17001184462689403

Thanks!

 

[mikrom]

is it possible to calculate the numerical derivative of M that comes out of the interpolation instead of making the numerical derivative of the data, thus avoiding having to interpolate it?

(pdata, Mdata) is your (x,y) data.
In your procedure
derivation(np,pdata,Mdata,dMdata)
you compute the derivative using the two step formula
f'(x) = (f(x+h) - f(x-h))/(2h)
where h is the step.
In the procedure
intlin1D(np,p,M,dM)
you interpolate the data linearly with the construction of a straight line, i.e. with a tangent where you compute the derivative using the one step formula
f'(x) = (f(x+h) - f(x))/h
Although the two step formula for computing derivative is numerically more stable, in my opinion using it in this case has not sense, because for the construction of the tangential line you are computing the derivative using the one step formula.
So, I would compute the derivative using one step formula like this:
dMdata(i) = (Mdata(i+1)-Mdata(i))/(pdata(i+1)-pdata(i))
M = Mdata(i)+ dMdata(i)*(p-pdata(i))

But it's only my opinion, it is question about numerical mathematics and not about fortran programming.

 

[R.Silveira]

I see ,
in my function i need both the numerical data and derivative of the numerical data, how can i interpolate the derivative dM ?

 

[mikrom]

your function [x, f(x)] is [pdata, Mdata]
the linear interpolation of the function in point p, where p(i) < p < p(i+1), is M
your derivative [x, f'(x)] is [pdata, dMdata]
the linear interpolation of the derivative in point p, where p(i) < p < p(i+1), is dM

 

[R.Silveira]

Yes that's right , so can i interpolate both at the same time ? or i need first interpolate M and later i call again to interpolate dM ?
sorry for this doubt, but I had never interpolated data before.