|
forums.silverfrost.com Welcome to the Silverfrost forums
|
View previous topic :: View next topic |
Author |
Message |
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Fri Feb 05, 2010 12:57 pm Post subject: 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?
|
|
Back to top |
|
|
sparge
Joined: 11 Apr 2005 Posts: 371
|
Posted: Fri Feb 05, 2010 10:44 pm Post subject: |
|
|
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)? |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Sun Feb 07, 2010 11:37 am Post subject: |
|
|
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).
Code: | +---+---+---+---+---+ +-*-+-*-*-+-+-+-*-+-*
| | | | | | | | | | | | | | | | |
+---+---+---+---+---+ +-+-+-+-+-+-+-+-+-+-+
| | | | | | => | | | | | | | | | | |
+---+---+---+---+---+ +-+-+-+-+-+-+-+-+-+-+
| | | | | | | | | | | | | | | | |
+---+---+---+---+---+ *-+-*-*-+-*-+-*-+-*-*
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. |
|
Back to top |
|
|
sparge
Joined: 11 Apr 2005 Posts: 371
|
Posted: Sun Feb 07, 2010 1:12 pm Post subject: Cave boundarem! |
|
|
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 |
|
Back to top |
|
|
IanLambley
Joined: 17 Dec 2006 Posts: 490 Location: Sunderland
|
Posted: Sun Feb 07, 2010 3:41 pm Post subject: |
|
|
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:
Code: |
!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 |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Tue Feb 09, 2010 2:44 pm Post subject: |
|
|
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:
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! |
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Tue Feb 09, 2010 4:06 pm Post subject: |
|
|
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:
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):
|
|
Back to top |
|
|
jjgermis
Joined: 21 Jun 2006 Posts: 404 Location: Nürnberg, Germany
|
Posted: Wed Feb 10, 2010 6:22 am Post subject: |
|
|
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!
Code: |
! 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)
|
|
|
Back to top |
|
|
DanRRight
Joined: 10 Mar 2008 Posts: 2818 Location: South Pole, Antarctica
|
Posted: Thu Feb 11, 2010 10:43 pm Post subject: |
|
|
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. |
|
Back to top |
|
|
DanRRight
Joined: 10 Mar 2008 Posts: 2818 Location: South Pole, Antarctica
|
Posted: Thu Feb 11, 2010 10:48 pm Post subject: |
|
|
|
|
Back to top |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|