! D02UBF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02ubfe_mod ! D02UBF Example Program Module: ! Parameter and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: four = 4.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: three = 3.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: two = 2.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: zero = 0.0_nag_wp INTEGER, PARAMETER :: m = 3, nin = 5, nout = 6 LOGICAL, PARAMETER :: reqerr = .FALSE. ! .. Local Scalars .. REAL (KIND=nag_wp) :: a, b CONTAINS FUNCTION exact(x,q) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: exact ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x INTEGER, INTENT (IN) :: q ! .. Local Scalars .. INTEGER :: lq ! .. Intrinsic Functions .. INTRINSIC cos, mod, sin ! .. Executable Statements .. lq = mod(q,4) SELECT CASE (lq) CASE (0) exact = cos(x) CASE (1) exact = -sin(x) CASE (2) exact = -cos(x) CASE (3) exact = sin(x) END SELECT RETURN END FUNCTION exact SUBROUTINE bndary(m,y,bmat,bvec) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: m ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: bmat(m,m+1), bvec(m), y(m) ! .. Executable Statements .. ! Boundary condition on left side of domain y(1:2) = a y(3) = b ! Set up Dirichlet condition using exact solution at x=a. bmat(1:m,1:m+1) = zero bmat(1:3,1) = one bmat(2:3,2) = two bmat(2:3,3) = three bvec(1) = zero bvec(2) = two bvec(3) = -two RETURN END SUBROUTINE bndary SUBROUTINE pdedef(m,f) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: m ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: f(m+1) ! .. Executable Statements .. f(1) = one f(2) = two f(3) = three f(4) = four RETURN END SUBROUTINE pdedef END MODULE d02ubfe_mod PROGRAM d02ubfe ! D02UBF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02uaf, d02ubf, d02ucf, d02uef, nag_wp, x01aaf, & x02ajf USE d02ubfe_mod, ONLY : a, b, bndary, exact, m, nin, nout, one, pdedef, & reqerr, two, zero ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: pi, resid, teneps INTEGER :: i, ifail, n, q, q1 ! .. Local Arrays .. REAL (KIND=nag_wp) :: bmat(m,m+1), bvec(m), f(m+1), & uerr(m+1), y(m) REAL (KIND=nag_wp), ALLOCATABLE :: c(:), f0(:), u(:,:), uc(:,:), x(:) ! .. Intrinsic Functions .. INTRINSIC abs, cos, int, max, sin ! .. Executable Statements .. WRITE (nout,*) ' D02UBF Example Program Results ' WRITE (nout,*) READ (nin,*) READ (nin,*) n ALLOCATE (u(n+1,m+1),f0(n+1),c(n+1),uc(n+1,m+1),x(n+1)) ! Set up problem boundary conditions and definition pi = x01aaf(one) a = -pi/two b = pi/two CALL bndary(m,y,bmat,bvec) CALL pdedef(m,f) ! Set up solution grid ifail = 0 CALL d02ucf(n,a,b,x,ifail) ! Set up problem right hand sides for grid and transform f0(1:n+1) = two*sin(x(1:n+1)) - two*cos(x(1:n+1)) ifail = 0 CALL d02uaf(n,f0,c,ifail) ! Solve in coefficient space ifail = 0 CALL d02uef(n,a,b,m,c,bmat,y,bvec,f,uc,resid,ifail) ! Transform solution and derivatives back to real space DO q = 0, m ifail = 0 CALL d02ubf(n,a,b,q,uc(1,q+1),u(1:n+1,q+1),ifail) END DO ! Print solution WRITE (nout,*) ' Numerical Solution U and its first three derivatives' WRITE (nout,*) WRITE (nout,99999) DO i = 1, n + 1 WRITE (nout,99998) x(i), u(i,1:m+1) END DO IF (reqerr) THEN uerr(1:m+1) = zero DO i = 1, n + 1 DO q = 0, m q1 = q + 1 uerr(q1) = max(uerr(q1),abs(u(i,q1)-exact(x(i),q))) END DO END DO teneps = 10.0_nag_wp*x02ajf() WRITE (nout,'(//)') DO q = 0, m WRITE (nout,99997) q, 10*(int(uerr(q+1)/teneps)+1) END DO END IF 99999 FORMAT (1X,T8,'X',T18,'U',T28,'Ux',T37,'Uxx',T47,'Uxxx') 99998 FORMAT (1X,5F10.4) 99997 FORMAT (1X,'Error in the order ',I1,' derivative of U is < ',I8, & ' * machine precision.') END PROGRAM d02ubfe