Problem with multiplying functions to create Lagrange interpolating polynomial basis

General Usage
1

I am trying to write some functions to generate a Lagrange interpolating polynomial basis. That is, I want to create a vector where each element is a lagrange polynomial basis function. I have some code below but it does not work. And I am not sure how to set this up to work correctly. I know that there is a package in julia for Lagrange polynomials, but I wanted to write this myself so that I understand them better.

what I would like is to have some vector of functions such that each element corresponds to a lagrange basis function. So if the vector is \psi, then \psi1should apply the first lagrange basis function at the value ofx`. Here is some code I am working with, but multiplying the functions to create the polynomial basis functions is proving harder than I through. Some subtle errors in how to handle functions and multiply them, etc. I was hoping someone could help. Here is some code that I am starting with.

function lagrange_polynomial(x, i, points)
    p(x) = 1;
    N = length(points)
    for k=1:N
        if k != i
            p = x -> p(x)*(x - points[k])/(points[i] - points[k])
            @show p(x)
        end
    end
    return p
end

function lagrange_polynomials_01(x, N)
    h = 1.0/N

    points = [(i-1)*h for i=1:N+1]
    ψ_array = Array{Function}(undef, N+1);
    for i=1:N+1
        ψ_array[i] = x -> lagrange_polynomial(x, i, points)
    end
    return ψ_array, points
end

So I should be able to do something like.

julia>psi, points = lagrange_polynomials_01(0.4, 2)
julia>psi[1](0.4)

But then I get the error:

ERROR: StackOverflowError:
Stacktrace:
 [1] (::var"#9#11"{Int64, Vector{Float64}, Int64})(x::Float64) (repeats 10823 times)
   @ Main /media/hayagriva/Dropbox/sandbox/julia/julia_pde/julia_finite_elements/langetangen/approx1d.jl:145
 [2] macro expansion
   @ ./show.jl:1128 [inlined]
 [3] lagrange_polynomial(x::Float64, i::Int64, points::Vector{Float64})
   @ Main /media/hayagriva/Dropbox/sandbox/julia/julia_pde/julia_finite_elements/langetangen/approx1d.jl:146

So not sure how to make this run smoothly. Seems like there are some simple function tricks that I am missing. Any help is appreciated.

2

How about this approach:

julia> lagrange_polynomial(x, i, points) = prod((x-points[k])/(points[i] - points[k]) for k in eachindex(points) if k ≠ i)
lagrange_polynomial (generic function with 1 method)

julia> points = 0:10
0:10

julia> lagrange_polynomials = [x -> lagrange_polynomial(x, i, points) for i in eachindex(points)]
11-element Vector{var"#42#44"{Int64}}:
 #42 (generic function with 1 method)
 #42 (generic function with 1 method)
 ⋮
 #42 (generic function with 1 method)

julia> lagrange_polynomials[3].(0:10)
11-element Vector{Float64}:
 -0.0
  0.0
  1.0
 -0.0
  0.0
 -0.0
  0.0
 -0.0
  0.0
 -0.0
  0.0