program main
!=======================================================================
! Chapter 10, code 1: a linear boundary value problem.
!
! Solve y'' + p(x)*y' + q(x)*y = r(x)
! with prescribed endpoint values y(a) = alpha and y(b) = beta.
!
! Method: centered finite differences on a uniform grid.
! Solve the resulting tridiagonal system with the Thomas algorithm.
!
! Example: y'' - y = 0, y(0) = y(1) = 1.
! Exact solution: cosh(x-0.5)/cosh(0.5).
!
! To solve another problem:
!   1. Edit the interval, boundary values, and number of intervals.
!   2. Edit coefficients(x,p,q,r) below the driver.
!   3. Disable check_exact or replace exact_solution.
!=======================================================================
implicit none
integer, parameter :: n = 20, data_unit = 20
! n intervals give n+1 grid points and n-1 interior unknowns.
double precision a, b, alpha, beta, h
double precision x(0:n), y(0:n)
double precision exact, error, max_error, exact_solution
integer i
logical check_exact
external coefficients, exact_solution

!-------------------- Problem and grid: edit here -----------------------
a = 0.0d0
b = 1.0d0
alpha = 1.0d0                 ! y(a)
beta = 1.0d0                  ! y(b)
check_exact = .true.
! Change n in the declaration above to refine the grid.
!-----------------------------------------------------------------------
if (n < 2 .or. b <= a) stop 'Use n >= 2 and b > a.'

h = (b-a)/dble(n)
do i = 0,n
    x(i) = a+dble(i)*h
end do
x(n) = b

call finite_difference(coefficients,n,x,alpha,beta,y)

open(unit=data_unit,file='ch10_finite_difference.dat', &
     status='replace',action='write')

write(*,'(a)') 'Linear boundary value problem: finite differences'
write(*,'(a,i8)') 'Grid intervals:', n
write(*,'(a,f12.6)') 'Grid spacing:  ', h
write(*,'(a)') ''

if (check_exact) then
    write(*,'(a10,2(1x,a16),1x,a12)') &
        'x','numerical y','exact y','abs error'
    write(data_unit,'(a)') '# x numerical_y exact_y abs_error'
else
    write(*,'(a10,1x,a16)') 'x','numerical y'
    write(data_unit,'(a)') '# x numerical_y'
end if

max_error = 0.0d0
do i = 0,n
    if (check_exact) then
        exact = exact_solution(x(i),a,b,alpha,beta)
        error = abs(y(i)-exact)
        max_error = max(max_error,error)
        write(*,'(f10.6,2(1x,es16.8),1x,es12.4)') &
            x(i),y(i),exact,error
        write(data_unit,'(4(1x,es24.16))') x(i),y(i),exact,error
    else
        write(*,'(f10.6,1x,es16.8)') x(i),y(i)
        write(data_unit,'(2(1x,es24.16))') x(i),y(i)
    end if
end do

close(data_unit)
write(*,'(a)') ''
if (check_exact) write(*,'(a,es12.4)') &
    'Largest grid error: ', max_error
write(*,'(a)') 'Data file: ch10_finite_difference.dat'
end program main


subroutine coefficients(x,p,q,r)
!=======================================================================
! Supply the coefficients of the original differential equation:
!   y'' + p(x)*y' + q(x)*y = r(x).
! Default example: y'' - y = 0.
!=======================================================================
implicit none
double precision x, p, q, r
p = 0.0d0
q = -1.0d0
r = 0.0d0
end subroutine coefficients


double precision function exact_solution(x,a,b,alpha,beta)
! Exact solution of y'' - y = 0 with the supplied endpoint values.
! Replace this function when changing the equation, or disable the check.
implicit none
double precision x, a, b, alpha, beta
exact_solution = (alpha*sinh(b-x)+beta*sinh(x-a))/sinh(b-a)
end function exact_solution


subroutine finite_difference(fcn,n,x,alpha,beta,y)
!=======================================================================
! Centered approximations:
!   y''(x_i) = (y_{i-1}-2*y_i+y_{i+1})/h**2
!   y'(x_i)  = (y_{i+1}-y_{i-1})/(2*h).
!
! Multiplying the differential equation by h**2 gives:
!   lower(i)*y_{i-1} + diagonal(i)*y_i + upper(i)*y_{i+1}
!       = h**2*r(x_i).
!
! Only the interior values are unknown. Move the known endpoint
! contributions to the right-hand side before solving.
!=======================================================================
implicit none
integer n, i
double precision x(0:n), alpha, beta, y(0:n)
double precision lower(n-1), diagonal(n-1), upper(n-1)
double precision rhs(n-1), interior(n-1)
double precision h, p, q, r
external fcn

h = (x(n)-x(0))/dble(n)
do i = 1,n-1
    call fcn(x(i),p,q,r)
    lower(i) = 1.0d0-h*p/2.0d0
    diagonal(i) = -2.0d0+h*h*q
    upper(i) = 1.0d0+h*p/2.0d0
    rhs(i) = h*h*r
end do

! Known boundary values contribute to the first and last equations.
! These two operations also work when there is only one interior point.
rhs(1) = rhs(1)-lower(1)*alpha
rhs(n-1) = rhs(n-1)-upper(n-1)*beta
lower(1) = 0.0d0
upper(n-1) = 0.0d0

call thomas(n-1,lower,diagonal,upper,rhs,interior)

y(0) = alpha
y(1:n-1) = interior
y(n) = beta
end subroutine finite_difference


subroutine thomas(n,lower,diagonal,upper,rhs,y)
!=======================================================================
! Solve a tridiagonal linear system by forward elimination and
! backward substitution. diagonal and rhs are overwritten.
!
! No pivoting is used. A near-zero pivot stops the calculation.
! Such a failure can require a solver with pivoting or can indicate
! a singular boundary value problem.
!=======================================================================
implicit none
integer n, i
double precision lower(n), diagonal(n), upper(n), rhs(n), y(n)
double precision multiplier, scale

! Check each pivot before division.
do i = 1,n
    scale = max(abs(lower(i)),abs(diagonal(i)),abs(upper(i)))
    if (abs(diagonal(i)) <= epsilon(1.0d0)*scale) &
        stop 'Thomas algorithm: zero or near-zero pivot.'
    if (i == n) exit
    multiplier = lower(i+1)/diagonal(i)
    diagonal(i+1) = diagonal(i+1)-multiplier*upper(i)
    rhs(i+1) = rhs(i+1)-multiplier*rhs(i)
end do

y(n) = rhs(n)/diagonal(n)
do i = n-1,1,-1
    y(i) = (rhs(i)-upper(i)*y(i+1))/diagonal(i)
end do
end subroutine thomas