c----------------------------------------------------------- c NOTE: The majority of this source code is derived from c code written by Dr. Matthew Choptuik. The original source c code can be found on wwwrel.ph.utexas.edu in the c ~phy329/linsys/ex2 directory. c----------------------------------------------------------- c=========================================================== c Solves 1-d linear boundary value problem c c u"(x) + a(x)u'(x) + b(x)u(x) = f(x) 0 <= x <= 1 c c using second-order finite difference technique and c LAPACK tridiagonal solver DGTSV. c=========================================================== program gbvp1d implicit none integer i4arg c----------------------------------------------------------- c Extrema of problem domain; note that this approach c of defining extrema as parameters makes it easier c to generalize program to arbitrary domains. c----------------------------------------------------------- real*8 xmin, xmax parameter ( xmin = 0.0d0, xmax = 1.0d0 ) c----------------------------------------------------------- c Maximum problem size (2**15 + 1) c----------------------------------------------------------- integer maxn parameter ( maxn = 32769 ) c----------------------------------------------------------- c Storage for discrete x-values, unknowns, exact c solution and right hand side values. a(maxn) and c b(maxn) are the constants in the O.D.E. c----------------------------------------------------------- real*8 x(maxn), u(maxn), & uexact(maxn), f(maxn), & a(maxn), b(maxn) c----------------------------------------------------------- c Storage for main, upper and lower diagonals of c tridiagonal system, and right-hand-side vector c for use with LAPACK routine DGTSV c----------------------------------------------------------- real*8 d(maxn), du(maxn), & dl(maxn), rhs(maxn), & cnst1(maxn), cnst2(maxn), cnst3(maxn) integer nrhs, info c----------------------------------------------------------- c Discretization level and size of system (# of discrete c unknowns) c----------------------------------------------------------- integer level, n, j, & option c----------------------------------------------------------- c Mesh spacing and related constants c----------------------------------------------------------- real*8 rmserr, h c----------------------------------------------------------- c Argument parsing. c----------------------------------------------------------- level = i4arg(1,-1) if( level .lt. 0 ) go to 900 n = 2 ** level + 1 if( n .gt. maxn ) then write(0,*) 'Insufficient internal storage' stop end if option = i4arg(2,0) c----------------------------------------------------------- c Set up finite-difference 'mesh' (discrete x-values) c and define some useful constants. c----------------------------------------------------------- h = 1.0d0 / (n - 1) do j = 1 , n x(j) = xmin + (j - 1) * h end do c----------------------------------------------------------- c Set up constants a(x) and b(x) c----------------------------------------------------------- do j = 1 , n a(j) = 16.0d0 * x(j)**2 b(j) = 8.0d0 * x(j)**3 cnst1(j) = (-a(j)/(2*h) + 1/(h*h)) cnst2(j) = (b(j) - 2/(h*h)) cnst3(j) = (a(j)/(2*h) + 1/(h*h)) enddo c----------------------------------------------------------- c This only ensures that x(n) = xmax EXACTLY and is not c essential. c----------------------------------------------------------- x(n) = xmax c----------------------------------------------------------- c Set up exact solution and right hand side vector. c----------------------------------------------------------- call exact(uexact,f,x,n) c=========================================================== c Set up tridiagonal system. Note that indexing on c lower diagonal is always (j-1) when implementing the c j'th equation. c=========================================================== c----------------------------------------------------------- c Left boundary: u(1) = u_0 c----------------------------------------------------------- d(1) = 1.0d0 du(1) = 0.0d0 rhs(1) = uexact(1) c----------------------------------------------------------- c Interior: Second order FDA of ODE. c----------------------------------------------------------- do j = 2 , n - 1 dl(j-1) = cnst1(j) d(j) = cnst2(j) du(j) = cnst3(j) rhs(j) = f(j) end do c----------------------------------------------------------- c Right boundary: u(n) = u_1 c----------------------------------------------------------- dl(n-1) = 0.0d0 d(n) = 1.0d0 rhs(n) = uexact(n) c=========================================================== c Solve tridiagonal system. c=========================================================== nrhs = 1 call dgtsv( n, nrhs, dl, d, du, rhs, n, info ) if( info .eq. 0 ) then c----------------------------------------------------------- c Solver successful, output either (x_j, u_j) or c (x_j, error_j) to stdout. Also compute rms error c and output to standard error. c----------------------------------------------------------- rmserr = 0.0d0 do j = 1 , n if( option .eq. 0 ) then write(*,*) x(j), rhs(j) else if (option .eq. 3) then write(*,*) x(j), uexact(j) else write(*,*) x(j), (uexact(j) - rhs(j)) end if rmserr = rmserr + (uexact(j) - rhs(j)) ** 2 end do rmserr = sqrt(rmserr / n) write(0,*) 'rmserr = ', rmserr else c----------------------------------------------------------- c Solver failed. c----------------------------------------------------------- write(0,*) 'gbvp1d: dgtsv() failed, info = ', info end if stop 900 continue write(0,*) 'usage: gbvp1d [