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.

### References

- Lagrange Interpolating Polynomial [Wolfram Mathworld]
- E. Waring, “Problems Concerning Interpolations“, Philosophical Transactions of the Royal Society of London, Vol. 69 (1779), p. 59-67.

## Leave a Reply