Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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