Lagrange Interpolating Polynomials

A simple Fortran implementation of interpolation by Lagrange polynomials is given below. Though named after Joseph-Louis Lagrange, the formula was first discovered by English mathematician Edward Waring.

subroutine lagrange_interpolation(xx,yy,x,y)
! Interpolate xx=[x1,x2,...xn]
!             yy=[f(x1),f(x2),...,f(xn)]
! Using Lagrange polynomial P(x)
! Returns y = P(x)
use,intrinsic :: iso_fortran_env, wp => real64

implicit none

real(wp),dimension(:),intent(in) :: xx
real(wp),dimension(:),intent(in) :: yy
real(wp),intent(in)              :: x
real(wp),intent(out)             :: y

!local variables:
integer :: j,k,n,m
real(wp) :: p

!check number of points:
n = size(xx)
m = size(yy)
if (n/=m) error stop &
    'Error: vectors must be the same size.'

!sum each of the Pj(x) terms:
y = 0.0_wp
do j=1,n
    !compute Pj(x):
    p = yy(j)
    do k=1,n
        if (k/=j) p = p * (x-xx(k)) / (xx(j)-xx(k))
    end do
    y = y + p
end do

end subroutine lagrange_interpolation

Examples of interpolating polynomials for the \(\sin(x)\) function are shown here for different numbers of points \(x=[0.0, 0.01]\), \(x=[0.0, 0.01, 0.02]\), etc.



Posted in Algorithms, Programming Tagged with: , , ,

Leave a Reply

Your email address will not be published. Required fields are marked *