Interpolation in two dimensions

I am using the **polint **subroutine from **Numerical Recipes in Fortran **to perform interpolation on a two-dimensional function (only increasing the number of columns). The interpolation works quite fine as long as I interpolate between the edges. On the 'real' edges the result is rather zig-zag. My variables are defined as real(kind=0.0) which must be sufficient for the numbers that I am working with. All the functions are monotonic and such a behaviour is not expected. My first thought are: 1.) precision is insufficient or 2.) input point lies slightly outside the boundary. What other possiblities could cause such an behaviour?

Your image is too small for me to see the problem. Are there real data points on the edges? And does the routine's interpolation pass through all the real data points or provide a smoothing fit? Does the routine perhaps offer more than one option for boundary condition and you are using the wrong sort for your particular problem? Perhaps if the routine allows you to specify zero gradient at the real edges you might get a better result (though I think this would only work if there are no data points on the real edges)?

Sorry for the small image: the web-site where I placed it automatically reduce it to that small size. To answer your question on the edges: yes, there are data points on the edges. Below is some simplified figure to explain what I do. I have a MxN matrix and to improve the searching algorithm afterwords I increase the number of columns. On the top and bottom edge I do not get a smooth line (where the * are).

+---+---+---+---+---+      +-*-+-*-*-+-+-+-*-+-*
|   |   |   |   |   |      | | | | | | | | | | | 
+---+---+---+---+---+      +-+-+-+-+-+-+-+-+-+-+
|   |   |   |   |   |  =>  | | | | | | | | | | | 
+---+---+---+---+---+      +-+-+-+-+-+-+-+-+-+-+
|   |   |   |   |   |      | | | | | | | | | | | 
+---+---+---+---+---+      *-+-*-*-+-*-+-*-+-*-*

M (rows) stays the same
N (columns) are increased 

I use the numerical recipes (http://www.nrbook.com/a/bookfpdf.php) interpolation routines. To be precise the explanation in (http://www.nrbook.com/a/bookfpdf/f3-6.pdf). There are no options and I think that I have take a more detailed look myself to see what is really happening. I would have expected the interpolation to behave good by default, since the original function is monotonic and smooth.

You don't say which of the four approaches you are using.

1. If you are using bilinear interpolation, which reduces to linear interpolation on the boundaries, and getting zig-zag behaviour there, then there is something wrong with the implementation.

2. If you're trying to do higher-order for accuracy, then the interpolation makes reference to a block of grid squares rather than just one, and you need to extend your original data set beyond the boundaries in some plausible way, otherwise you'll be introducing garbage into the computation along the boundaries.

3. If you're trying to do higher order for smoothness using bicubic interpolation, you escape the need to reference fictitious additional data points, but you have to specify gradient information. That gives a lot of options, and since you say there are none, I presume you are not doing that.

4. If you're trying to do higher order for smoothness using bicubic splining, the calculation gets the gradient information for you, but you have to provide fictitious additional data points again - so the same caveat as for option 2).

Have you run the code with /checkmate to make sure you aren't accessing beyond the boundaries of your data array? Note the code for option 2) says the ya array is taken for granted, so if you're successfully passing it partial garbage on the boundaries, it won't know.

Andy

I have just been looking at an interpolation I did for a similar grid arrangement and I find that I have a fix to avoid selecting edge values. So I must have had a problem with the edge interpolation. I'm using a four point Lagrange interpolation. All I did was to limit the x value for the interpolation to:

!calculate the width of the grid
xrange = x(max_row) - x(1)
!calculate a small fraction of the width
xlimit = xrange/(max_row*1d6)
x_use = max(x(1)+xlimit ,min(x(max_row)-xlimit ,x_desired))

Then use this slightly modified value in the interpolation.

Then do exactly the same in the z direction and find the Y value.

If you want to see the whole code, please send me a personal message with your e-mail address. Regards Ian

Due to the monotonic behviour I use polynomial interpolation. I think that I first have to find the interval and perform a polinomial interpolation (linear) between the two points. I actually tried to fit some high oder polynomial to some that that are almost linear. If one does this, the result looks something like this:

I did implement the the suggestion from Ian (I think it is a good suggestion). Here is my code:

    id_range = idFI(1,size(idFI,2)) - idFI(1,1)
    id_limit = id_range/(size(idFI,2)*1.0e4) 
!    
    iq_range = iqFI(size(iqFI,1),1) - iqFI(1,1)
    iq_limit = iq_range/(size(iqFI,1)*1.0e4)     
    do i=1,nnew
      id_ref = idFI(1,1)+delta_id_new*(i-1)
      id_use = max(idFI(1,1)+id_limit ,min(idFI(1,size(idFI,2))-id_limit ,id_ref))
      do j=1,matrix%m
!       interpolate along iq      
        iq_ref = iqFI(j,1)
        iq_use = max(iqFI(1,1)+iq_limit ,min(iqFI(size(iqFI,1),1)-iq_limit ,iq_ref))
!       d- and q-current
        iqF(j,i) = iq_ref
        idF(j,i) = id_ref    
!       TORQUE        
        call polin2(iqFI(:,1),idFI(1,:), TFI,  iq_use,id_use ,TF(j,i),    dy)  
      end do
    enddo

I think that the 2D interpolation is not that easy as I thought it to be! Playing with x_limit realised that my polynomial interpolation might not be the best option!

Here the final solution: Keep the polynomial interpolation function. Find the interval where x lies and then use these two points, i.e. linear interpolation. The code:

! Step over the rows
  do j=1,m
!   Select a row
    yntmp=ya(j,:)
    do i=1,ndum-1
!     Find the interval for x2       
      if (x2>=x2a(i).AND.x2<(x2a(i+1))) then
        call polint(x2a(i:i+1),yntmp(i:i+1),x2,ymtmp(j),dy)
        exit
      endif
    enddo
  end do
! The result from stepping over the rows is a column
  do i=1,m-1
    if (x1>=x1a(i).AND.x1<(x1a(i+1))) then
      call polint(x1a(i:i+1),ymtmp(i:i+1),x1,y,dy)
      exit
    endif
  enddo

And here the final result (increased N from 25 to 40):

Maybe some last comment: Here is the original code that I used. As one can see, I used all the points in a row and performed polinomial interpolation (and not linear as I wanted). Anyway 1.) the suggestion from Ian made me play with the limit value and 2.) the questions from Andy made me look at the description of the function in more detail.

In doing so I could understand my problem better and thereby get it to do what I want it to do in the first place. Moreover, the total calculation time has decreased as well!

! Step over the rows 
  do j=1,m 
!   Select a row 
    yntmp=ya(j,:) 
    call polint(x2a,yntmp,x2,ymtmp(j),dy)
  end do 
! The result from stepping over the rows is a column 
  call polint(x1a,ymtmp,x1,y,dy) 

jjgermis,

One of examples i've posted (one which changes perspective) has probably what you are asking for. As an option, when you resize the image, it is doing linear interpolation from one 2D mesh (which is image) to another. Not sure i've paid super attention to boundary problem, please doublecheck.