Skip to content

Instantly share code, notes, and snippets.

@ivan-pi
Created March 18, 2019 00:59
Show Gist options
  • Save ivan-pi/1b476dabd11b651fc2ea6e4fd0c11289 to your computer and use it in GitHub Desktop.
Save ivan-pi/1b476dabd11b651fc2ea6e4fd0c11289 to your computer and use it in GitHub Desktop.
Code example from book "Chebyshev and Fourier Spectral Methods" by John P. Boyd
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment