- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I need to do some simple linear interpolation in a Fortran code. I struggle to find the right subroutine and a simple example on how to do it.
For example, Let x(i) = 1.0, 2.0 ,3.0 ,4.0 ,5.0 and y(i) = 1.0, 4.0, 9.0, 16.0, 25.0
I want to plug in x= 3.75 and return a value for y.
Yes, I know the example above is y=x**2, but my real problem is more complicated.
Can someone point me to the correct mkl subroutine and the syntax for calling it?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi William,
Has the information provided helped? Could you please let us know if we could close this case on our end?
Best Regards,
Shanmukh.SS
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Linear interpolation in 1-D is a trivial task, so it is unlikely that anyone would create a library routine for it.
- Step through the sorted array x and find the index i such that x(i) <= t < x(i+1). If no such index is found, either do extrapolation or declare failure. If the array is large, you may use binary search instead of stepping through the x array sequentially.
- Calculate the interpolated value s at x = t using s = y(i) + [y(i+1) - y(i)] / [x(i+1) - x(i)] * [t - x(i)]
You stated that your real problem is more complicated. Please elaborate.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
This is more a question for the MKL Forum but
you can try the MKL function for interpolation. It uses splines, so it's more complicated than a linear interpolation.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I did not know there was a mkl forum. I will look for it.
Elaborating on the complexity of the problem, the equation I am trying to operate on is
y= x/a + b*(x/c)**n
I need to solve for the slope (dy/dx) as a function of y only. I was unable to generate an analytical solution. So my plan was to generate data points that has y as a function of x, then calculate dy/dx as a function of y, then interpolate that. I thought there might be a mkl function to make life easier. But I can just use the routine suggested by the previous respondent. In any event, thanks for the help.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Using interpolation to calculate derivatives amplifies any interpolation errors. You should use a method that guarantees continuous derivatives and discourages oscillation from overfitting. Cubic splines can be a good choice.
Fortunately your expression y= x/a + b*(x/c)**n has an analytic derivative. (I cheated and used the symbolic algebra package maxima.) dy/dx is the final expression below. Checked it with a few random values.
Maxima 5.46.0 https://maxima.sourceforge.io
using Lisp CLISP 2.49+ (2010-07-17)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) ex:x/a + b*(x/c)^n;
x n x
(%o1) b (-) + -
c a
(%i2) diff(ex,x);
x n
b n (-)
c 1
(%o2) -------- + -
x a
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
maxima can generate the Fortran expression, too. Trivial in this case, but useful for longer expressions. The fortran() function generates fixed format. There is also an f90() function that generates free-format.
(%i1) y:x/a + b*(x/c)^n;
x n x
(%o1) b (-) + -
c a
(%i2) dydx:diff(y,x);
x n
b n (-)
c 1
(%o2) -------- + -
x a
(%i3) fortran(dydx);
(b*n*(x/c)**n)/x+1/a
(%o3) done
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Here's the link to the MKL Forum, if you need more help, although there are knowledgeable MKL users who hang out here.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Unless the values of a, b, c and n, and the range of interest in x conspire to make things simple, the expression that you revealed is not amenable to linear interpolation. If the nonlinearity is high, linear interpolation can be quite inaccurate.
The variable p = dy/dx may not be monotonic in x, so you have to be in possession of a rule to select one value based on some external reason such as physics, etc. Similarly, if the range of x contains 0 and n is not an integer, more study is needed.
Much can be learned by making a phase-plane plot of p versus y and inspecting the plot. This should be done before selecting a numerical routine (such as one of the spline interpolation routines available in MKL).
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
You could also use automatic differentiation. It is a fairly simple technique - see for instance https://sourceforge.net/p/flibs/svncode/HEAD/tree/trunk/src/computing/automdiff.f90. But an even simpler technique is to use complex numbers. Here is a small program I wrote to experiment with this technique:
! derivative.f90 --
! Use the complex step method to determine the derivative of a real function
!
! Method described by Nicholas Higham
!
! With simple expressions the result is quite intriguingly accurate. Note though the
! last example - there you get a considerable fluctuation even with h = 1.0e-20 to 1.0e-30
!
program derivative
implicit none
real(kind=kind(1.0d0)) :: h
complex(kind=kind(1.0d0)) :: x, y
integer :: i
do i = 1,10
h = 0.1d0 ** i
x = cmplx(1.0d0,h,kind=kind(1.0d0))
write(*,*) imag( sin(x) ) / h, cos(real(x)), (imag( sin(x) ) / h - cos(real(x))) / cos(real(x))
write(*,*) imag( cos(x) ) / h, -sin(real(x)), (imag( cos(x) ) / h + sin(real(x))) / (-sin(real(x)))
enddo
do i = 1,60
h = 0.1d0 ** (3*i)
x = cmplx(1.0d0,h,kind=kind(1.0d0))
y = x + x**2 + x**3
write(*,*) imag( y ) / h, h
enddo
do i = 1,60
h = 0.1d0 ** (3*i)
x = cmplx(1.0d0,h,kind=kind(1.0d0))
y = cos(x)/(1.0d0 - sin(x))
write(*,*) imag( y ) / h, h
enddo
end program derivative
See also this thread on Fortran discourse: https://fortran-lang.discourse.group/t/fortran-code-snippets/2150
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi William,
Has the information provided helped? Could you please let us know if we could close this case on our end?
Best Regards,
Shanmukh.SS
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
One quick way to work out what function may fit, before you get stuck into Fortran is to look at the data in EXCEL. It has many fine interpolation functions.
Here is your example:
Then you use Fortran. of course if the data is complex, like the COVID death data, you use a FFT routine to look for the waves in the data. A French group of Dr missed a peak in one of their published data sets because they did not use FFT.
It is an interesting field in stats, but you have to know what you are going to program if you want to save a lot of time.
Good luck. PS You got excellent help in the posts above.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thanks. I actually tried that approach and got close, but not close enough as the tail of the data trailed off. Could just not find the right functional form
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
>> not close enough as the tail of the data trailed off.
That may be an artifact of the signal (or simulation) not stabilizing and/or your sampling rate being too slow to catch the higher frequencies that may be present at the time on an event (equivalent to simulation integration step size being to long).
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Give a data sample, you have, excluding myself, some of the best Fortran programmers and excellent applied math people. I just like to hang out with the bright people at the nerds table.
Not all data fits one of the few patterns we understand, most Gaussian data is not really Gaussian it is just close enough that we accept the difference and get on with fun programming.
There is a book - How round is my circle, I was thinking about it in the library the other day as I watched a lady cut our Easter eggs, I suddenly thought - do a linear transformation and then linear regression. For @mecej4 - solving the implicit problem in the title of the book.
Or as my pure math teacher once said if you look at some obscure figure in the 23rd dimension through the 22nd dimension you will see a donut. If that guy had said the moon was made of blue cheese I would have accepted it directly, I failed his exam five times and finally gave up only to be told the exam was out of his textbook which was not on the syllabus. Never do it alone. As I explained to my 16 yr old who failed a French exam at school when she thought I would be upset, Darling I failed exams, it is not the knocking down that is important it is the getting back up.
And I have a theory that @mecej4 is really the Fonz.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi William,
Glad to know that your issue is resolved. If you need any additional information, please post a new question as this thread will no longer be monitored by Intel and further discussion within the forum is considered community only.
Have a great day!
Best Regards,
Shanmukh.SS

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page