program main
!====================================================================
! Runge phenomenon and Chebyshev nodes
!
! The Runge function f(x) is sampled using two sets of n base points:
! equally spaced points and Chebyshev points. A global polynomial
! interpolant is then evaluated for both cases at nint points across
! the interval. The program also calculates the average absolute
! interpolation error for each choice of nodes.
!
! The interpolating polynomial is evaluated in Lagrange form.
!
! This example illustrates that increasing the polynomial order with
! equally spaced points can produce large oscillations near the ends
! of the interval, while Chebyshev nodes greatly reduce this effect.
!
! Alexander Godunov
! Prepared for the companion website, 2026.
!====================================================================
implicit none

integer, parameter :: n=13      ! base points for interpolation
integer, parameter :: nint=25   ! compute interpolation in nint points

double precision, parameter :: pi=3.14159265358979323846d0

double precision xmin, xmax
double precision, dimension(n) :: xe, ye, xc, yc
double precision x, y, step
double precision yse, ysc, erre, errc
double precision errav_e, errav_c
integer i
double precision f, lagrange_global

xmin = -1.0d0
xmax =  1.0d0

! Step 1: generate equally spaced interpolation points
step = (xmax-xmin)/dble(n-1)

do i=1,n
  xe(i) = xmin + step*dble(i-1)
  ye(i) = f(xe(i))
end do

! Step 2: generate Chebyshev interpolation points
! The nodes are written in increasing order from left to right.
do i=1,n
  xc(i) = cos((2.0d0*dble(n-i)+1.0d0)*pi/(2.0d0*dble(n)))
  yc(i) = f(xc(i))
end do

! Step 3: evaluate both global interpolants at nint points
errav_e = 0.0d0
errav_c = 0.0d0
step = (xmax-xmin)/dble(nint-1)

write(*,*) '     Runge Function: Equally Spaced and Chebyshev Nodes'
write(*,203) n
write(*,201)

do i=1,nint
  x = xmin + step*dble(i-1)
  y = f(x)

  yse = lagrange_global(x, xe, ye, n)
  ysc = lagrange_global(x, xc, yc, n)

  erre = yse-y
  errc = ysc-y

  write(*,200) x, yse, erre, ysc, errc

! Step 4: calculate the average absolute interpolation errors
  errav_e = errav_e + abs(erre)/dble(nint)
  errav_c = errav_c + abs(errc)/dble(nint)
end do

write(*,202) errav_e, errav_c

200 format (5f13.6)
201 format ('        x    equal poly   equal error', &
            '    Cheb poly    Cheb error')
202 format (/,' Average absolute error:',/, &
            '   equally spaced = ',f12.6,/, &
            '   Chebyshev      = ',f12.6)
203 format ('     number of interpolation points = ',i2)

end program main


!
! Function f(x): Runge function
!
function f(x)
implicit none
double precision f, x

f = 1.0d0/(1.0d0 + 25.0d0*x*x)

end function f


function lagrange_global(xx, xi, yi, n)
!====================================================================
! Global polynomial interpolation in Lagrange form.
!
! All n data points are used to construct and evaluate the
! interpolating polynomial.
!
! Input:
!   xx   = abscissa at which the interpolation is evaluated
!   xi   = array of data abscissas
!   yi   = array of data ordinates
!   n    = number of interpolation points
!
! Output:
!   lagrange_global = interpolated value at xx
!
! Comments:
!   The routine works for both equally and unequally spaced points.
!   Unlike the local interpolation routine used elsewhere in this
!   chapter, all available points are used here so that the effect
!   of node placement on a high-order global polynomial is visible.
!
! Alexander Godunov
! Prepared for the companion website, 2026.
!====================================================================
implicit none

double precision lagrange_global, xx
integer n
double precision xi(n), yi(n)

integer i, j
double precision lambda, y

y = 0.0d0

do i=1,n
  lambda = 1.0d0

  do j=1,n
    if (j /= i) then
      lambda = lambda*(xx-xi(j))/(xi(i)-xi(j))
    end if
  end do

  y = y + yi(i)*lambda
end do

lagrange_global = y

end function lagrange_global