! D02UEF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02uefe_mod ! D02UEF Example Program Module: ! Parameters 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) ! .. Function Return Value .. Real (Kind=nag_wp) :: exact ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: x Integer, Intent (In) :: q ! .. Intrinsic Procedures .. Intrinsic :: cos, sin ! .. Executable Statements .. Select Case (q) Case (0) exact = cos(x) Case (1) exact = -sin(x) Case (2) exact = -cos(x) Case (3) exact = sin(x) End Select End Function exact Subroutine bndary(m,y,bmat,bvec) ! .. 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 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) ! .. 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 d02uefe_mod Program d02uefe ! D02UEF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02uaf, d02ubf, d02ucf, d02uef, nag_wp, x01aaf, & x02ajf Use d02uefe_mod, Only: a, b, bndary, exact, m, nin, nout, 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 Procedures .. Intrinsic :: abs, cos, int, max, sin ! .. Executable Statements .. Write (nout,*) ' D02UEF 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 domain, boundary conditions and definition pi = x01aaf(zero) 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) ! Evaluate solution and derivatives on Chebyshev grid. Do q = 0, m ifail = 0 Call d02ubf(n,a,b,q,uc(1,q+1),u(1,q+1),ifail) End Do ! Print solution Write (nout,*) ' Numerical Solution U and its first three derivatives' Write (nout,*) Write (nout,99999) Write (nout,99998)(x(i),u(i,1:m+1),i=1,n+1) 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,'(//)') Write (nout,99997)(q,10*(int(uerr(q+1)/teneps)+1),q=0,m) 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 d02uefe