Code example from book "Chebyshev and Fourier Spectral Methods" by John P. Boyd

boyd_bvp_example.f90

program boyd_bvp

    dimension xi(20),aphi(20),g(20),h(20,20),ugraph(101)

    common/pbasis/phi(20),phix(20),phixx(20)
    common/thiele/th

    pi = 4.*atan(1.0)

    ! specify parameters
    th = 12.
    alpha = 1.
    beta = 1.
    n = 22
    nbasis = n - 2

    ! compute the interior collocation points and the forcing vector g
    do i = 1, nbasis
        xi(i) = cos(pi*i/(nbasis+1))
        x = xi(i)
        b = alpha*(1-x)/2. + beta*(1+x)/2.
        bx = (-alpha + beta)/2.
        g(i) = f(x) - d0(x)*b - d1(x)*bx
    end do

    ! compute the square matrix
    do i = 1, nbasis
        x = xi(i)
        call basis(x,nbasis)
        dd0 = d0(x)
        dd1 = d1(x)
        dd2 = d2(x)
        do j = 1, nbasis
            h(i,j) = dd2*phixx(j) + dd1*phix(j) + dd0*phi(j)
        end do
    end do

    ! call a linear algebra routine to solve the matrix equation h * aphi = g

    call linsolv(h,g,aphi,nbasis)

    ! the array aphi now contains the nbasis coefficients of the
    ! expansion v(x) in terms of the basis functions phi_j(x).
    ! We may convert these coefficients into those of the ordinary
    ! Chebyshev series as described in Sec. 5, but this is optional.

    ! Make a graph of u(x) at 101 points by summing the basis function
    ! [not Chebyshev] series.

    ugraph(1) = alpha
    ugraph(101) = beta

    do i = 2, 100
        x = -1 + (i-1)*0.02
        call basis(x,nbasis)
        ! Recall that u(x) = v(x) + B(x) where v(x) is given by the
        ! computed series and B(x) is the R.H.S. of the next line.
        ugraph(i) = (alpha*(1-x) + beta*(1+x))/2.
        do j = 1, nbasis
            ugraph(i) = ugraph(i) + aphi(j)*phi(j)
        end do
    end do

    write(*,*) -1., ugraph(1)
    do i = 2, 100
        x = -1 + (i-1)*0.02
        write(*,*) x, ugraph(i)
    end do
    write(*,*) 1., ugraph(101)

end program

subroutine basis(x,nbasis)
common/pbasis/phi(20),phix(20),phixx(20)
! After this call, the arrays phi, phix, and phixx contain the
! values of the basis functions (and their first two derivatives,
! respectively, at x.)
if (abs(x) < 1) then
    t = acos(x)
    c = cos(t)
    s = sin(t)
    do i = 1, nbasis
        n = i + 1
        tn = cos(n*t)
        tnt = -n*sin(n*t)
        tntt = -n*n*tn

        ! Convert t-derivatives into x-derivatives
        
        tnx = -tnt/s
        tnxx = tntt/(s*s) - tnt*c/(s*s*s)

        ! Final step: convert T_Ns into the basis functions phi_Ns.
        ! We subtract 1 from the even degree T_Ns, x from the odd degree
        ! T_N and 1 from the first derivative of the odd polynomials only

        if (mod(n,2) == 0) then
            phi(i) = tn - 1.
            phix(i) = tnx
        else
            phi(i) = tn - x
            phix(i) = tnx - 1.
        end if
        phixx(i) = tnxx
    end do
else
    do i = 1, nbasis
        phi(i) = 0.
        n = i + 1
        if (mod(n,2) == 0) then
            phix(i) = sign(x,1.)*n*n
        else
            phix(i) = n*n - 1
        end if
        phixx(i) = (sign(x,1.))**n * n*n * (n*n-1.)/3.
    end do
end if
end subroutine

real function d2(x)
    d2 = 1.
end function

real function d1(x)
    d1 = 0.
end function

real function d0(x)
    common/thiele/th

    d0 = -th*th
end function

real function f(x)
    f = 0.0
end function

subroutine linsolv(a,b,x,n)
    implicit none
    integer, intent(in) :: n
    real, intent(inout) :: a(n,n)
    real, intent(in) :: b(n)
    real, intent(out) :: x(n)

    integer :: ipiv(n)
    integer :: info

    call sgetrf(n,n,a,n,ipiv,info)
    if (info /= 0) then
        write(*,*) "Factorization failed with flag ", info
        stop 1
    end if
        
    x = b ! copy right-hand side
    call sgetrs('N',n,1,a,n,ipiv,x,n,info)
    if (info /= 0) then
        write(*,*) "Solution failed with flag ", info
        stop 1
    end if
end subroutine