c=========================================================== c HISTORY: bvp1d.f. Refer to that code for more c extensive documentation of the basic solution c procedure c c Solves 1-d linear boundary value problem c c u''(x) = f(x) on x = [0,1]; u(0) = u0, u(1) = u1 c c using base second-order finite difference technique c and LAPACK tridiagonal solver DGTSV. Solves system c using two levels of discretization then uses c Richardson extrapolation to produce O(h^4) solution. c=========================================================== program rexbvp1d implicit none integer i4arg real*8 xmin, xmax parameter ( xmin = 0.0d0, xmax = 1.0d0 ) integer maxn parameter ( maxn = 1 048 577 ) c----------------------------------------------------------- c Extrapolation constants c----------------------------------------------------------- real*8 c13, c43 parameter ( c13 = -0.3333 3333 3333 3333 d0, & c43 = 1.3333 3333 3333 3333 d0 ) c----------------------------------------------------------- c Storage for discrete x-values, unknowns, exact c solution and right hand side values. Note that c these arrays now have an additional dimension since c we want to generate solutions on two distinct grids. c Also note that we declare them (maxn,2), rather than c (2,maxn) so that values on a specific mesh remain c contiguous. c----------------------------------------------------------- real*8 x(maxn,2), u(maxn,2), & uexact(maxn,2), f(maxn,2) 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) integer nrhs, info c----------------------------------------------------------- c Base discretization level, size of system, output c option and locals. c----------------------------------------------------------- integer blevel, n, option, & level, j, lev c----------------------------------------------------------- c Mesh spacing and related constants (1/h**2, -2/h**2) c----------------------------------------------------------- real*8 h, hm2, m2hm2 real*8 rmserr c----------------------------------------------------------- c Argument parsing. c----------------------------------------------------------- blevel = i4arg(1,-1) c----------------------------------------------------------- c Base level must be .ge. 1 since we also need to c compute a solution at level (blevel - 1). c----------------------------------------------------------- if( blevel .lt. 1 ) go to 900 n = 2 ** blevel + 1 if( n .gt. maxn ) then write(0,*) 'Insufficient internal storage' stop end if option = i4arg(2,0) c----------------------------------------------------------- c OUTER LOOP: Compute solution for levels blevel and c blevel - 1 c----------------------------------------------------------- do lev = 1 , 2 c----------------------------------------------------------- c Set up finite-difference 'mesh' (discrete x-values) c and define some useful constants. c----------------------------------------------------------- level = blevel + 1 - lev n = 2 ** level + 1 h = 1.0d0 / (n - 1) do j = 1 , n x(j,lev) = xmin + (j - 1) * h end do hm2 = 1.0d0 / (h * h) m2hm2 = -2.0d0 / (h * h) x(n,lev) = xmax c----------------------------------------------------------- c Set up exact solution and right hand side vector. c----------------------------------------------------------- call exact(uexact(1,lev),f(1,lev),x(1,lev),n) c----------------------------------------------------------- c Set up tridiagonal system. c----------------------------------------------------------- d(1) = 1.0d0 du(1) = 0.0d0 rhs(1) = uexact(1,lev) do j = 2 , n - 1 dl(j-1) = hm2 d(j) = m2hm2 du(j) = hm2 rhs(j) = f(j,lev) end do dl(n-1) = 0.0d0 d(n) = 1.0d0 rhs(n) = uexact(n,lev) 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, save solution for subsequent c processing. c----------------------------------------------------------- do j = 1 , n u(j,lev) = rhs(j) end do else c----------------------------------------------------------- c Solver failed: print error message and exit. c----------------------------------------------------------- write(0,*) 'bvp1d: dgtsv() failed, info = ',info stop end if c----------------------------------------------------------- c END OF OUTER LOOP c----------------------------------------------------------- end do c----------------------------------------------------------- c Richardson extrapolate. Note that without inter- c polation, solution can only be computed on coarse c grid (blevel - 1). c----------------------------------------------------------- n = 2 ** (blevel - 1 ) + 1 do j = 1 , n u(j,2) = c43 * u(2*(j - 1) + 1,1) + c13 * u(j,2) end do c----------------------------------------------------------- c Compute RMS error and output either solution or c error (on coarse grid) to standard output, rms error c to standard error. c----------------------------------------------------------- rmserr = 0.0d0 do j = 1 , n if( option .eq. 0 ) then write(*,*) x(j,2), u(j,2) else write(*,*) x(j,2), (uexact(j,2) - u(j,2)) end if rmserr = rmserr + (uexact(j,2) - u(j,2)) ** 2 end do rmserr = sqrt(rmserr / n) write(0,*) 'rmserr = ', rmserr stop 900 continue write(0,*) 'usage: rexbvp1d [