Try out NAG Library functions

Explore NAG maths and stats routines with interactive demos
Function ID
D03PCF
Name
nagf_pde_1d_parab_fd_old
Description
General system of parabolic PDEs, method of lines, finite differences, one space variable
Keywords
method of lines | parabolic partial differential equation
We use the example given in Dew and Walsh (1981) which consists of an elliptic-parabolic pair of PDEs. The problem was originally derived from a single third-order in space PDE. The elliptic equation is
1rrr2U1r=4αU2+rU2r  
and the parabolic equation is
1-r2U2t=1rrrU2r-U2U1  
where r,t0,1×0,1. The boundary conditions are given by
U1=U2r=0  at ​r=0,  
and
rrU1=0  and  U2=0  at ​r=1.  
The example data reflects that shown in the "Example" section of the routine documentation. You can change this here to try alternative inputs. The formatting will need to be kept as it is here, otherwise the program is likely to fail to run correctly.

Please note that incompatible data will however cause the example output to display an error message. These error messages are fully explained in the Routine document
!   D03PCF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    Module d03pcfe_mod

!     D03PCF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: bndary, pdedef, uinit
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6, npde = 2
!     .. Local Scalars ..
      Real (Kind=nag_wp), Public, Save :: alpha
    Contains
      Subroutine pdedef(npde,t,x,u,ux,p,q,r,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t, x
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: npde
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: p(npde,npde), q(npde), r(npde)
        Real (Kind=nag_wp), Intent (In) :: u(npde), ux(npde)
!       .. Executable Statements ..
        q(1) = 4.0_nag_wp*alpha*(u(2)+x*ux(2))
        q(2) = 0.0_nag_wp
        r(1) = x*ux(1)
        r(2) = ux(2) - u(1)*u(2)
        p(1,1) = 0.0_nag_wp
        p(1,2) = 0.0_nag_wp
        p(2,1) = 0.0_nag_wp
        p(2,2) = 1.0_nag_wp - x*x
        Return
      End Subroutine pdedef
      Subroutine bndary(npde,t,u,ux,ibnd,beta,gamma,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: ibnd, npde
        Integer, Intent (Inout)        :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: beta(npde), gamma(npde)
        Real (Kind=nag_wp), Intent (In) :: u(npde), ux(npde)
!       .. Executable Statements ..
        If (ibnd==0) Then
          beta(1) = 0.0_nag_wp
          beta(2) = 1.0_nag_wp
          gamma(1) = u(1)
          gamma(2) = -u(1)*u(2)
        Else
          beta(1) = 1.0_nag_wp
          beta(2) = 0.0_nag_wp
          gamma(1) = -u(1)
          gamma(2) = u(2)
        End If
        Return
      End Subroutine bndary
      Subroutine uinit(u,x,npts)

!       Routine for PDE initial condition

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: u(2,npts)
        Real (Kind=nag_wp), Intent (In) :: x(npts)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Executable Statements ..
        Do i = 1, npts
          u(1,i) = 2.0_nag_wp*alpha*x(i)
          u(2,i) = 1.0_nag_wp
        End Do
        Return
      End Subroutine uinit
    End Module d03pcfe_mod
    Program d03pcfe

!     D03PCF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d03pcf, d03pzf, nag_wp, x01aaf
      Use d03pcfe_mod, Only: alpha, bndary, nin, nout, npde, pdedef, uinit
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: acc, hx, pi, piby2, tout, ts
      Integer                          :: i, ifail, ind, intpts, it, itask,    &
                                          itrace, itype, lisave, lrsave, m,    &
                                          neqn, npts, nwk
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: rsave(:), u(:,:), uout(:,:,:), x(:), &
                                          xout(:)
      Integer, Allocatable             :: isave(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real, sin
!     .. Executable Statements ..
      Write (nout,*) 'D03PCF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) intpts, npts, itype
      neqn = npde*npts
      lisave = neqn + 24
      nwk = (10+6*npde)*neqn
      lrsave = nwk + (21+3*npde)*npde + 7*npts + 54
      Allocate (rsave(lrsave),u(npde,npts),uout(npde,intpts,itype),x(npts),    &
        xout(intpts),isave(lisave))

      Read (nin,*) xout(1:intpts)
      Read (nin,*) acc, alpha
      Read (nin,*) m, itrace
      ind = 0
      itask = 1

!     Set spatial mesh points

      piby2 = 0.5_nag_wp*x01aaf(pi)
      hx = piby2/real(npts-1,kind=nag_wp)
      x(1) = 0.0_nag_wp
      x(npts) = 1.0_nag_wp
      Do i = 2, npts - 1
        x(i) = sin(hx*real(i-1,kind=nag_wp))
      End Do

!     Set initial conditions
      Read (nin,*) ts, tout

!     Set the initial values
      Call uinit(u,x,npts)

      Do it = 1, 5
        tout = 10.0_nag_wp*tout

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d03pcf(npde,m,ts,tout,pdedef,bndary,u,npts,x,acc,rsave,lrsave,    &
          isave,lisave,itask,itrace,ind,ifail)

        If (it==1) Then
          Write (nout,99999) acc, alpha
          Write (nout,99998) xout(1:6)
        End If

!       Interpolate at required spatial points

        ifail = 0
        Call d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)

        Write (nout,99996) tout, uout(1,1:intpts,1)
        Write (nout,99995) uout(2,1:intpts,1)
      End Do

!     Print integration statistics

      Write (nout,99997) isave(1), isave(2), isave(3), isave(5)

99999 Format (/,/,' Accuracy requirement  = ',E12.5,/,' Parameter ALPHA =',    &
        '       ',E12.3,/)
99998 Format ('   T  /  X   ',6F8.4,/)
99997 Format (' Number of integration steps in time                  ',I4,/,   &
        ' Number of residual evaluations of resulting ODE system',I4,/,        &
        ' Number of Jacobian evaluations                        ',I4,/,        &
        ' Number of iterations of nonlinear solver              ',I4)
99996 Format (1X,F7.4,' U(1)',6F8.4)
99995 Format (9X,'U(2)',6F8.4,/)
    End Program d03pcfe
The NAG Library
The world’s largest collection of robust, documented, tested and maintained numerical algorithms.
Learn more here or contact us for purchasing information