Lagrange Interpolation Function Generator

Specific Domains Modelling & Simulations
1

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?

Thanks

#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
2 Likes
2

Maybe Interpolations.jl?

I don’t think it’s necessary to make the array here. You can just chain(1:i-1,i+1:N).

1 Like
3

I would recommend using the more-stable barycentric interpolation formula rather than the lagrange formula.

Also, I would tend to just return an InterpolatingPolynomial type rather than an anonymous function, and then overload the function-call operation so you can call it like a function. That way you can define additional methods on the interpolating polynomial, e.g. to compute the derivative.

Also, your function is not type-stable. You shouldn’t adopt the Matlab style of having the same function do vectors and scalars, but rather use separate methods. Once you define a scalar function p(x), you can apply it to arrays with p.(x), though in some cases you may want to define a specialized vector method for efficiency. (Your implementation is not particularly efficient because it generates lots of temporary arrays.)

6 Likes
4

Hi! im a little bit late to the topic. What does
** numvalstoevaluate == 1 ? output = sum(LagrangePolynomials.yvals) : output[k] = sum(LagrangePolynomials.yvals)

do? i mean im not sure what the “?” means between the variable numvalstoevaluate and output, nor how does the output function works :cold_sweat: