×
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Contact US

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Doubt about interpolation and integration with derivatives

Doubt about interpolation and integration with derivatives

Doubt about interpolation and integration with derivatives

(OP)
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

CODE --> 90

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

CODE --> 90

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

CODE --> data

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!

RE: Doubt about interpolation and integration with derivatives

Quote (R.Silveira)


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.


RE: Doubt about interpolation and integration with derivatives

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


RE: Doubt about interpolation and integration with derivatives

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

RE: Doubt about interpolation and integration with derivatives

(OP)
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.



Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members! Already a Member? Login


Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close