I couldn’t find an implementation of Lagrange Interpolation in Julia, so built my own. This function simply takes a set of points, as stored in the two vectors xvals and yvals, and spits out the Lagrange Polynomial that passes through these points. I’ve attached the code below, along with an example. I’m sure it needs some tidying up but thought I’d share in case it’s useful - unlike the spectral decompositions in ApproxFun, lower-order Lagrange polynomials can often be used (with care) for small extrapolations. Are there any basic interpolation packages who might appreciate this (tidied-up) code as a pull request?
#This function returns another function, which is the Lagrange Interpolant of the values xvals and yvals. function LagrangeInterpolantGenerator(xvals,yvals) function LagrangeInterpolant(x) numvalstoevaluate = length(x) numvalstoevaluate == 1 ? output = 0 : output = zeros(numvalstoevaluate) for k = 1:numvalstoevaluate N = length(xvals) LagrangePolynomials = ones(N) for i in 1:N for j in [1:i-1;i+1:N] #Surprisingly, this works even in the i=1 and i=N cases. LagrangePolynomials[i] = LagrangePolynomials[i].*(x[k]-xvals[j])./(xvals[i]-xvals[j]) end end numvalstoevaluate == 1 ? output = sum(LagrangePolynomials.*yvals) : output[k] = sum(LagrangePolynomials.*yvals) end return output end return LagrangeInterpolant end #Examples interpolantfunc = LagrangeInterpolantGenerator([1.0,2.0,3.0],[4.0,2.0,8.0]) interpolantfunc(2.0) #returns 2.0 interpolantfunc(3.0) #returns 8.0 interpolantfunc(2.5) #returns 5.375