hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_pde_2d_gen_order2_rectilinear (d03rb)

Purpose

nag_pde_2d_gen_order2_rectilinear (d03rb) integrates a system of linear or nonlinear, time-dependent partial differential equations (PDEs) in two space dimensions on a rectilinear domain. The method of lines is employed to reduce the PDEs to a system of ordinary differential equations (ODEs) which are solved using a backward differentiation formula (BDF) method. The resulting system of nonlinear equations is solved using a modified Newton method and a Bi-CGSTAB iterative linear solver with ILU preconditioning. Local uniform grid refinement is used to improve the accuracy of the solution. nag_pde_2d_gen_order2_rectilinear (d03rb) originates from the VLUGR2 package (see Blom and Verwer (1993) and Blom et al. (1996)).

Syntax

[ts, dt, rwk, iwk, ind, ifail] = d03rb(ts, tout, dt, tols, tolt, inidom, pdedef, bndary, pdeiv, monitr, opti, optr, rwk, iwk, lenlwk, itrace, ind, 'npde', npde, 'lenrwk', lenrwk, 'leniwk', leniwk)
[ts, dt, rwk, iwk, ind, ifail] = nag_pde_2d_gen_order2_rectilinear(ts, tout, dt, tols, tolt, inidom, pdedef, bndary, pdeiv, monitr, opti, optr, rwk, iwk, lenlwk, itrace, ind, 'npde', npde, 'lenrwk', lenrwk, 'leniwk', leniwk)

Description

nag_pde_2d_gen_order2_rectilinear (d03rb) integrates the system of PDEs:
Fj (t,x,y,u,ut,ux,uy,uxx,uxy,uyy) = 0 ,   j = 1,2,,npde ,   (x,y) Ω ,   t0 t tout ,
Fj (t,x,y,u,ut,ux,uy,uxx,uxy,uyy) = 0 ,   j=1,2,,npde ,   (x,y) Ω ,   t0 t tout ,
(1)
where ΩΩ is an arbitrary rectilinear domain, i.e., a domain bounded by perpendicular straight lines. If the domain is rectangular then it is recommended that nag_pde_2d_gen_order2_rectangle (d03ra) is used.
The vector uu is the set of solution values
u(x,y,t) = [u1(x,y,t),,unpde(x,y,t)]T,
u(x,y,t)=[u1(x,y,t),,unpde(x,y,t)]T,
and utut denotes partial differentiation with respect to tt, and similarly for uxux, etc.
The functions FjFj must be supplied by you in pdedef. Similarly the initial values of the functions u(x,y,t)u(x,y,t) for (x,y)Ω(x,y)Ω must be specified at t = t0t=t0 in pdeiv.
Note that whilst complete generality is offered by the master equations (1), nag_pde_2d_gen_order2_rectilinear (d03rb) is not appropriate for all PDEs. In particular, hyperbolic systems should not be solved using this function. Also, at least one component of utut must appear in the system of PDEs.
The boundary conditions must be supplied by you in bndary in the form
Gj (t,x,y,u,ut,ux,uy) = 0 ,   j = 1,2,,npde ,   (x,y)Ω ,   t0ttout .
Gj (t,x,y,u,ut,ux,uy) = 0 ,   j=1,2,,npde ,   (x,y)Ω ,   t0ttout .
(2)
The domain is covered by a uniform coarse base grid specified by you, and nested finer uniform subgrids are subsequently created in regions with high spatial activity. The refinement is controlled using a space monitor which is computed from the current solution and a user-supplied space tolerance tols. A number of optional parameters, e.g., the maximum number of grid levels at any time, and some weighting factors, can be specified in the arrays opti and optr. Further details of the refinement strategy can be found in Section [Further Comments].
The system of PDEs and the boundary conditions are discretized in space on each grid using a standard second-order finite difference scheme (centred on the internal domain and one-sided at the boundaries), and the resulting system of ODEs is integrated in time using a second-order, two-step, implicit BDF method with variable step size. The time integration is controlled using a time monitor computed at each grid level from the current solution and a user-supplied time tolerance tolt, and some further optional user-specified weighting factors held in optr (see Section [Further Comments] for details). The time monitor is used to compute a new step size, subject to restrictions on the size of the change between steps, and (optional) user-specified maximum and minimum step sizes held in dt. The step size is adjusted so that the remaining integration interval is an integer number times ΔtΔt. In this way a solution is obtained at t = toutt=tout.
A modified Newton method is used to solve the nonlinear equations arising from the time integration. You may specify (in opti) the maximum number of Newton iterations to be attempted. A Jacobian matrix is calculated at the beginning of each time step. If the Newton process diverges or the maximum number of iterations is exceeded, a new Jacobian is calculated using the most recent iterates and the Newton process is restarted. If convergence is not achieved after the (optional) user-specified maximum number of new Jacobian evaluations, the time step is retried with Δt = Δt / 4Δt=Δt/4. The linear systems arising from the Newton iteration are solved using a Bi-CGSTAB iterative method, in combination with ILU preconditioning. The maximum number of iterations can be specified by you in opti.
In order to define the base grid you must first specify a virtual uniform rectangular grid which contains the entire base grid. The position of the virtual grid in physical (x,y)(x,y) space is given by the (x,y)(x,y) coordinates of its boundaries. The number of points nxnx and nyny in the xx and yy directions must also be given, corresponding to the number of columns and rows respectively. This is sufficient to determine precisely the (x,y)(x,y) coordinates of all virtual grid points. Each virtual grid point is then referred to by integer coordinates (vx,vy)(vx,vy), where (0,0)(0,0) corresponds to the lower-left corner and (nx1,ny1)(nx-1,ny-1) corresponds to the upper-right corner. vxvx and vyvy are also referred to as the virtual column and row indices respectively.
The base grid is then specified with respect to the virtual grid, with each base grid point coinciding with a virtual grid point. Each base grid point must be given an index, starting from 11, and incrementing row-wise from the leftmost point of the lowest row. Also, each base grid row must be numbered consecutively from the lowest row in the grid, so that row 11 contains grid point 11.
As an example, consider the domain consisting of the two separate squares shown in Figure 1. The left-hand diagram shows the virtual grid and its integer coordinates (i.e., its column and row indices), and the right-hand diagram shows the base grid point indices and the base row indices (in brackets).
Figure 1
Figure 1
Hence the base grid point with index 66 say is in base row 22, virtual column 4, and virtual row 11, i.e., virtual grid integer coordinates (4,1)(4,1); and the base grid point with index 1919 say is in base row 55, virtual column 2, and virtual row 55, i.e., virtual grid integer coordinates (2,5)(2,5).
The base grid must then be defined in inidom by specifying the number of base grid rows, the number of base grid points, the number of boundaries, the number of boundary points, and the following integer arrays:
Finally, ilbnd contains the types of the boundaries and corners, as follows:
Boundaries:
External corners (90 ° 90°):
Internal corners (270 ° 270°):
Figure 2 shows the boundary types of a domain with a hole. Notice the logic behind the labelling of the corners: each one includes the types of the two adjacent boundary edges, in a clockwise fashion (outside the domain).
Figure 2
Figure 2
As an example, consider the domain shown in Figure 3. The left-hand diagram shows the physical domain and the right-hand diagram shows the base and virtual grids. The numbers outside the base grid are the indices of the left and rightmost base grid points, and the numbers inside the base grid are the boundary or corner numbers, indicating the order in which the boundaries are stored in lbnd.
Figure 3
Figure 3
For this example we have
NROWS = 11
NPTS = 105
NBNDS = 28
NBPTS = 72

LROW = (1,4,15,26,37,46,57,68,79,88,97)

IROW = (0,1,2,3,4,5,6,7,8,9,10)

ICOL = (0,1,2,
        0,1,2,3,4,5,6,7,8,9,10,
        0,1,2,3,4,5,6,7,8,9,10,
        0,1,2,3,4,5,6,7,8,9,10,
        0,1,2,3,4,5,8,9,10,
        0,1,2,3,4,5,6,7,8,9,10,
        0,1,2,3,4,5,6,7,8,9,10,
        0,1,2,3,4,5,6,7,8,9,10,
        0,1,2,3,4,5,6,7,8,
        0,1,2,3,4,5,6,7,8,
        0,1,2,3,4,5,6,7,8)

LBND = (2,
        4,15,26,37,46,57,68,79,88,
        98,99,100,101,102,103,104,
        96,
        86,85,84,83,82,
        70,59,48,39,28,17,6,
        8,9,10,11,12,13,
        18,29,40,49,60,
        72,73,74,75,76,77,
        67,56,45,36,25,
        33,32,
        42,
        52,53,
        43,
        1,97,105,87,81,3,7,71,78,14,31,51,54,34)

LLBND = (1,2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,
         61,62,63,64,65,66,67,68,69,70,71,72)

ILBND = (1,2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,
         12,23,34,41,43,14,21,32)
This particular domain is used in the example in Section [Example], and data statements are used to define the above arrays in that example program. For less complicated domains it is simpler to assign the values of the arrays in do-loops. This also allows flexibility in the number of base grid points.
The function nag_pde_2d_gen_order2_checkgrid (d03ry) can be called from inidom to obtain a simple graphical representation of the base grid, and to verify the data that you have specified in inidom.
Subgrids are stored internally using the same data structure, and solution information is communicated to you in pdeiv, pdedef and bndary in arrays according to the grid index on the particular level, e.g., x(i)xi and y(i)yi contain the (x,y)(x,y) coordinates of grid point ii, and u(i,j)uij contains the jjth solution component ujuj at grid point ii.
The grid data and the solutions at all grid levels are stored in the workspace arrays, along with other information needed for a restart (i.e., a continuation call). It is not intended that you extract the solution from these arrays, indeed the necessary information regarding these arrays is not provided. The user-supplied monitor (monitr) should be used to obtain the solution at particular levels and times. monitr is called at the end of every time step, with the last step being identified via the input parameter tlast. The function nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) should be called from monitr to obtain grid information at a particular level.
Further details of the underlying algorithm can be found in Section [Further Comments] and in Blom and Verwer (1993) and Blom et al. (1996) and the references therein.

References

Blom J G, Trompert R A and Verwer J G (1996) Algorithm 758. VLUGR2: A vectorizable adaptive grid solver for PDEs in 2D Trans. Math. Software 22 302–328
Blom J G and Verwer J G (1993) VLUGR2: A vectorized local uniform grid refinement code for PDEs in 2D Report NM-R9306 CWI, Amsterdam
Trompert R A (1993) Local uniform grid refinement and systems of coupled partial differential equations Appl. Numer. Maths 12 331–355
Trompert R A and Verwer J G (1993) Analysis of the implicit Euler local uniform grid refinement method SIAM J. Sci. Comput. 14 259–278

Parameters

Compulsory Input Parameters

1:     ts – double scalar
The initial value of the independent variable tt.
Constraint: ts < toutts<tout.
2:     tout – double scalar
The final value of tt to which the integration is to be carried out.
3:     dt(33) – double array
The initial, minimum and maximum time step sizes respectively.
dt(1)dt1
Specifies the initial time step size to be used on the first entry, i.e., when ind = 0ind=0. If dt(1) = 0.0dt1=0.0 then the default value dt(1) = 0.01 × (toutts)dt1=0.01×(tout-ts) is used. On subsequent entries (ind = 1ind=1), the value of dt(1)dt1 is not referenced.
dt(2)dt2
Specifies the minimum time step size to be attempted by the integrator. If dt(2) = 0.0dt2=0.0 the default value dt(2) = 10.0 × machine precisiondt2=10.0×machine precision is used.
dt(3)dt3
Specifies the maximum time step size to be attempted by the integrator. If dt(3) = 0.0dt3=0.0 the default value dt(3) = touttsdt3=tout-ts is used.
Constraints:
  • if ind = 0ind=0, dt(1)0.0dt10.0;
  • if ind = 0ind=0 and dt(1) > 0.0dt1>0.0,
    10.0 × machine precision × max (|ts|,|tout|)dt(1)toutts10.0×machine precision×max(|ts|,|tout|)dt1tout-ts and
    dt(2)dt(1)dt(3)dt2dt1dt3, where the values of dt(2)dt2 and dt(3)dt3 will have been reset to their default values if zero on entry;
  • 0 dt(2) dt(3) 0 dt2 dt3 .
4:     tols – double scalar
The space tolerance used in the grid refinement strategy (σσ in equation (4)). See Section [Refinement Strategy].
Constraint: tols > 0.0tols>0.0.
5:     tolt – double scalar
The time tolerance used to determine the time step size (ττ in equation (7)). See Section [Time Integration ].
Constraint: tolt > 0.0tolt>0.0.
6:     inidom – function handle or string containing name of m-file
inidom must specify the base grid in terms of the data structure described in Section [Description]. inidom is not referenced if, on entry, ind = 1ind=1. nag_pde_2d_gen_order2_checkgrid (d03ry) can be called from inidom to obtain a simple graphical representation of the base grid, and to verify the data that you have specified in inidom. nag_pde_2d_gen_order2_rectilinear (d03rb) also checks the validity of the data, but you are strongly advised to call nag_pde_2d_gen_order2_checkgrid (d03ry) to ensure that the base grid is exactly as required.
Note: the boundaries of the base grid should consist of as many points as are necessary to employ second-order space discretization, i.e., a boundary enclosing the internal part of the domain must include at least 33 grid points including the corners. If Neumann boundary conditions are to be applied the minimum is 44.
[xmin, xmax, ymin, ymax, nx, ny, npts, nrows, nbnds, nbpts, lrow, irow, icol, llbnd, ilbnd, lbnd, ierr] = inidom(maxpts, ierr)

Input Parameters

1:     maxpts – int64int32nag_int scalar
The maximum number of base grid points allowed by the available workspace.
2:     ierr – int64int32nag_int scalar
Will be initialized by nag_pde_2d_gen_order2_rectilinear (d03rb) to some value prior to internal calls to inidom.

Output Parameters

1:     xmin – double scalar
2:     xmax – double scalar
The extents of the virtual grid in the xx-direction, i.e., the xx coordinates of the left and right boundaries respectively.
3:     ymin – double scalar
4:     ymax – double scalar
The extents of the virtual grid in the yy-direction, i.e., the yy coordinates of the left and right boundaries respectively.
5:     nx – int64int32nag_int scalar
6:     ny – int64int32nag_int scalar
The number of virtual grid points in the xx- and yy-direction respectively (including the boundary points).
7:     npts – int64int32nag_int scalar
The total number of points in the base grid. If the required number of points is greater than maxpts then inidom must be exited immediately with ierr set to 1-1 to avoid overwriting memory.
8:     nrows – int64int32nag_int scalar
The total number of rows of the virtual grid that contain base grid points. This is the maximum base row index.
9:     nbnds – int64int32nag_int scalar
The total number of physical boundaries and corners in the base grid.
10:   nbpts – int64int32nag_int scalar
The total number of boundary points in the base grid.
11:   lrow( : :) – int64int32nag_int array
lrow(i)lrowi, for i = 1,2,,nrowsi=1,2,,nrows, must contain the base grid index of the first grid point in base grid row ii.
12:   irow( : :) – int64int32nag_int array
irow(i)irowi, for i = 1,2,,nrowsi=1,2,,nrows, must contain the virtual row number vyvy that corresponds to base grid row ii.
13:   icol( : :) – int64int32nag_int array
icol(i)icoli, for i = 1,2,,nptsi=1,2,,npts, must contain the virtual column number vxvx that contains base grid point ii.
14:   llbnd( : :) – int64int32nag_int array
llbnd(i)llbndi, for i = 1,2,,nbndsi=1,2,,nbnds, must contain the element of lbnd corresponding to the start of the iith boundary or corner.
Note: the order of the boundaries and corners in llbnd must be first all the boundaries and then all the corners. The end points of a boundary (i.e., the adjacent corner points) must not be included in the list of points on that boundary. Also, if a corner is shared by two pairs of physical boundaries then it has two types and must therefore be treated as two corners.
15:   ilbnd( : :) – int64int32nag_int array
ilbnd(i)ilbndi, for i = 1,2,,nbndsi=1,2,,nbnds, must contain the type of the iith boundary (or corner), as given in Section [Description].
16:   lbnd( : :) – int64int32nag_int array
lbnd(i)lbndi, for i = 1,2,,nbptsi=1,2,,nbpts, must contain the grid index of the iith boundary point. The order of the boundaries is as specified in llbnd, but within this restriction the order of the points in lbnd is arbitrary.
17:   ierr – int64int32nag_int scalar
If the required number of grid points is larger than maxpts, ierr must be set to 1-1 to force a termination of the integration and an immediate return to the calling program with ifail = 3ifail=3. Otherwise, ierr should remain unchanged.
7:     pdedef – function handle or string containing name of m-file
pdedef must evaluate the functions FjFj, for j = 1,2,,npdej=1,2,,npde, in equation (1) which define the system of PDEs (i.e., the residuals of the resulting ODE system) at all interior points of the domain. Values at points on the boundaries of the domain are ignored and will be overwritten by bndary. pdedef is called for each subgrid in turn.
[res] = pdedef(npts, npde, t, x, y, u, ut, ux, uy, uxx, uxy, uyy)

Input Parameters

1:     npts – int64int32nag_int scalar
The number of grid points in the current grid.
2:     npde – int64int32nag_int scalar
The number of PDEs in the system.
3:     t – double scalar
The current value of the independent variable tt.
4:     x(npts) – double array
x(i)xi contains the xx coordinate of the iith grid point, for i = 1,2,,nptsi=1,2,,npts.
5:     y(npts) – double array
y(i)yi contains the yy coordinate of the iith grid point, for i = 1,2,,nptsi=1,2,,npts.
6:     u(npts,npde) – double array
u(i,j)uij contains the value of the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
7:     ut(npts,npde) – double array
ut(i,j)utij contains the value of (u)/(t) u t  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
8:     ux(npts,npde) – double array
ux(i,j)uxij contains the value of (u)/(x) u x  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
9:     uy(npts,npde) – double array
uy(i,j)uyij contains the value of (u)/(y) u y  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
10:   uxx(npts,npde) – double array
uxx(i,j)uxxij contains the value of (2u)/(x2) 2u x2  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
11:   uxy(npts,npde) – double array
uxy(i,j)uxyij contains the value of (2u)/(xy) 2u xy  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
12:   uyy(npts,npde) – double array
uyy(i,j)uyyij contains the value of (2u)/(y2) 2u y2  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.

Output Parameters

1:     res(npts,npde) – double array
res(i,j)resij must contain the value of FjFj, for j = 1,2,,npdej=1,2,,npde, at the iith grid point, for i = 1,2,,nptsi=1,2,,npts, although the residuals at boundary points will be ignored (and overwritten later on) and so they need not be specified here.
8:     bndary – function handle or string containing name of m-file
bndary must evaluate the functions GjGj, for j = 1,2,,npdej=1,2,,npde, in equation (2) which define the boundary conditions at all boundary points of the domain. Residuals at interior points must not be altered by this function.
[res] = bndary(npts, npde, t, x, y, u, ut, ux, uy, nbnds, nbpts, llbnd, ilbnd, lbnd, res)

Input Parameters

1:     npts – int64int32nag_int scalar
The number of grid points in the current grid.
2:     npde – int64int32nag_int scalar
The number of PDEs in the system.
3:     t – double scalar
The current value of the independent variable tt.
4:     x(npts) – double array
x(i)xi contains the xx coordinate of the iith grid point, for i = 1,2,,nptsi=1,2,,npts.
5:     y(npts) – double array
y(i)yi contains the yy coordinate of the iith grid point, for i = 1,2,,nptsi=1,2,,npts.
6:     u(npts,npde) – double array
u(i,j)uij contains the value of the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
7:     ut(npts,npde) – double array
ut(i,j)utij contains the value of (u)/(t) u t  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
8:     ux(npts,npde) – double array
ux(i,j)uxij contains the value of (u)/(x) u x  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
9:     uy(npts,npde) – double array
uy(i,j)uyij contains the value of (u)/(y) u y  for the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
10:   nbnds – int64int32nag_int scalar
The total number of physical boundaries and corners in the grid.
11:   nbpts – int64int32nag_int scalar
The total number of boundary points in the grid.
12:   llbnd(nbnds) – int64int32nag_int array
llbnd(i)llbndi, for i = 1,2,,nbndsi=1,2,,nbnds, contains the element of lbnd corresponding to the start of the iith boundary (or corner).
13:   ilbnd(nbnds) – int64int32nag_int array
ilbnd(i)ilbndi, for i = 1,2,,nbndsi=1,2,,nbnds, contains the type of the iith boundary, as given in Section [Description].
14:   lbnd(nbpts) – int64int32nag_int array
lbnd(i)lbndi, contains the grid index of the iith boundary point, where the order of the boundaries is as specified in llbnd. Hence the iith boundary point has coordinates x(lbnd(i))xlbndi and y(lbnd(i))ylbndi, and the corresponding solution values are u(lbnd(i),j)ulbndij, for i = 1,2,,nbptsi=1,2,,nbpts and j = 1,2,,npdej=1,2,,npde.
15:   res(npts,npde) – double array
Contains function values returned by pdedef.

Output Parameters

1:     res(npts,npde) – double array
res(lbnd(i),j)reslbndij must contain the value of GjGj, for j = 1,2,,npdej=1,2,,npde, at the iith boundary point, for i = 1,2,,nbptsi=1,2,,nbpts.
Note:  elements of res corresponding to interior points, i.e., points not included in lbnd, must not be altered.
9:     pdeiv – function handle or string containing name of m-file
pdeiv must specify the initial values of the PDE components uu at all points in the base grid. pdeiv is not referenced if, on entry, ind = 1ind=1.
[u] = pdeiv(npts, npde, t, x, y)

Input Parameters

1:     npts – int64int32nag_int scalar
The number of grid points in the base grid.
2:     npde – int64int32nag_int scalar
The number of PDEs in the system.
3:     t – double scalar
The (initial) value of the independent variable tt.
4:     x(npts) – double array
x(i)xi contains the xx coordinate of the iith grid point, for i = 1,2,,nptsi=1,2,,npts.
5:     y(npts) – double array
y(i)yi contains the yy coordinate of the iith grid point, for i = 1,2,,nptsi=1,2,,npts.

Output Parameters

1:     u(npts,npde) – double array
u(i,j)uij must contain the value of the jjth PDE component at the iith grid point, for i = 1,2,,nptsi=1,2,,npts and j = 1,2,,npdej=1,2,,npde.
10:   monitr – function handle or string containing name of m-file
monitr is called by nag_pde_2d_gen_order2_rectilinear (d03rb) at the end of every successful time step, and may be used to examine or print the solution or perform other tasks such as error calculations, particularly at the final time step, indicated by the parameter tlast.
The input arguments contain information about the grid and solution at all grid levels used. nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) should be called from monitr in order to extract the number of points and their (x,y)(x,y) coordinates on a particular grid.
monitr can also be used to force an immediate tidy termination of the solution process and return to the calling program.
[ierr] = monitr(npde, t, dt, dtnew, tlast, nlev, xmin, ymin, dxb, dyb, lgrid, istruc, lsol, sol, ierr)

Input Parameters

1:     npde – int64int32nag_int scalar
The number of PDEs in the system.
2:     t – double scalar
The current value of the independent variable tt, i.e., the time at the end of the integration step just completed.
3:     dt – double scalar
The current time step size ΔtΔt, i.e., the time step size used for the integration step just completed.
4:     dtnew – double scalar
The time step size that will be used for the next time step.
5:     tlast – logical scalar
Indicates if intermediate or final time step. tlast = falsetlast=false for an intermediate step, tlast = truetlast=true for the last call to monitr before returning to your program.
6:     nlev – int64int32nag_int scalar
The number of grid levels used at time t.
7:     xmin – double scalar
8:     ymin – double scalar
The (x,y)(x,y) coordinates of the lower-left corner of the virtual grid.
9:     dxb – double scalar
10:   dyb – double scalar
The sizes of the base grid spacing in the xx- and yy-direction respectively.
11:   lgrid( : :) – int64int32nag_int array
Contains pointers to the start of the grid structures in istruc, and must be passed unchanged to nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) in order to extract the grid information.
12:   istruc( : :) – int64int32nag_int array
Contains the grid structures for each grid level and must be passed unchanged to nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) in order to extract the grid information.
13:   lsol(nlev) – int64int32nag_int array
lsol(l)lsoll contains the pointer to the solution in sol at grid level ll and time t. (lsol(l)lsoll actually contains the array index immediately preceding the start of the solution in sol.)
14:   sol(lenrwk(6 × npde + 1)lenrwk-(6×npde+1)) – double array
Contains the solution uu at time t for each grid level ll in turn, positioned according to lsol. More precisely
u(i,j) = sol(lsol(l) + (j1) × nl + i)
uij=sol(lsoll+(j-1)×nl+i)
represents the jjth component of the solution at the iith grid point in the llth level, for i = 1,2,,nli=1,2,,nl, j = 1,2,,npdej=1,2,,npde and l = 1,2,,nlevl=1,2,,nlev, where nlnl is the number of grid points at level ll (obtainable by a call to nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz)).
15:   ierr – int64int32nag_int scalar
Will be initialized by nag_pde_2d_gen_order2_rectilinear (d03rb) to some value prior to internal calls to ierr.

Output Parameters

1:     ierr – int64int32nag_int scalar
Should be set to 11 to force a termination of the integration and an immediate return to the calling program with ifail = 4ifail=4. ierr should remain unchanged otherwise.
11:   opti(44) – int64int32nag_int array
May be set to control various options available in the integrator.
opti(1) = 0opti1=0
All the default options are employed.
opti(1) > 0opti1>0
The default value of opti(i)optii, for i = 2,3,4i=2,3,4, can be obtained by setting opti(i) = 0optii=0.
opti(1)opti1
Specifies the maximum number of grid levels allowed (including the base grid). opti(1)0opti10. The default value is opti(1) = 3opti1=3.
opti(2)opti2
Specifies the maximum number of Jacobian evaluations allowed during each nonlinear equations solution. opti(2)0opti20. The default value is opti(2) = 2opti2=2.
opti(3)opti3
Specifies the maximum number of Newton iterations in each nonlinear equations solution. opti(3)0opti30. The default value is opti(3) = 10opti3=10.
opti(4)opti4
Specifies the maximum number of iterations in each linear equations solution. opti(4)0opti40. The default value is opti(4) = 100opti4=100.
Constraint: opti(1)0opti10 and if opti(1) > 0opti1>0, opti(i)0optii0, for i = 2,3,4i=2,3,4.
12:   optr(33,npde) – double array
May be used to specify the optional vectors umaxumax, wsws and wtwt in the space and time monitors (see Section [Further Comments]).
If an optional vector is not required then all its components should be set to 1.01.0.
optr(1,j)optr1j, for j = 1,2,,npdej=1,2,,npde, specifies ujmaxujmax, the approximate maximum absolute value of the jjth component of uu, as used in (4) and (7). optr(1,j) > 0.0optr1j>0.0, for j = 1,2,,npdej=1,2,,npde.
optr(2,j)optr2j, for j = 1,2,,npdej=1,2,,npde, specifies wjswjs, the weighting factors used in the space monitor (see (4)) to indicate the relative importance of the jjth component of uu on the space monitor. optr(2,j)0.0optr2j0.0, for j = 1,2,,npdej=1,2,,npde.
optr(3,j)optr3j, for j = 1,2,,npdej=1,2,,npde, specifies wjtwjt, the weighting factors used in the time monitor (see (6)) to indicate the relative importance of the jjth component of uu on the time monitor. optr(3,j)0.0optr3j0.0, for j = 1,2,,npdej=1,2,,npde.
Constraints:
  • optr(1,j) > 0.0optr(1,j)>0.0, for j = 1,2,,npdej=1,2,,npde;
  • optr(i,j)0.0optrij0.0, for i = 2,3i=2,3 and j = 1,2,,npdej=1,2,,npde.
13:   rwk(lenrwk) – double array
The required value of lenrwk cannot be determined exactly in advance, but a suggested value is
lenrwk = maxpts × npde × (5 × l + 18 × npde + 9) + 2 × maxpts,
lenrwk=maxpts×npde×(5×l+18×npde+9)+2×maxpts,
where l = opti(1)l=opti1 if opti(1)0opti10 and l = 3l=3 otherwise, and maxptsmaxpts is the expected maximum number of grid points at any one level. If during the execution the supplied value is found to be too small then the function returns with ifail = 3ifail=3 and an estimated required size is printed on the current error message unit (see nag_file_set_unit_error (x04aa)).
Note:  the size of lenrwk cannot be checked upon initial entry to nag_pde_2d_gen_order2_rectilinear (d03rb) since the number of grid points on the base grid is not known.
14:   iwk(leniwk) – int64int32nag_int array
If ind = 0ind=0, iwk need not be set. Otherwise iwk must remain unchanged from a previous call to nag_pde_2d_gen_order2_rectilinear (d03rb).
15:   lenlwk – int64int32nag_int scalar
The dimension of the array lwk as declared in the (sub)program from which nag_pde_2d_gen_order2_rectilinear (d03rb) is called.
The required value of lenlwk cannot be determined exactly in advance, but a suggested value is
lenlwk = maxpts + 1,
lenlwk=maxpts+1,
where maxptsmaxpts is the expected maximum number of grid points at any one level. If during the execution the supplied value is found to be too small then the function returns with ifail = 3ifail=3 and an estimated required size is printed on the current error message unit (see nag_file_set_unit_error (x04aa)).
Note:  the size of lenlwk cannot be checked upon initial entry to nag_pde_2d_gen_order2_rectilinear (d03rb) since the number of grid points on the base grid is not known.
16:   itrace – int64int32nag_int scalar
The level of trace information required from nag_pde_2d_gen_order2_rectilinear (d03rb). itrace may take the value 1-1, 00, 11, 22 or 33.
itrace = 1itrace=-1
No output is generated.
itrace = 0itrace=0
Only warning messages are printed.
itrace > 0itrace>0
Output from the underlying solver is printed on the current advisory message unit (see nag_file_set_unit_advisory (x04ab)). This output contains details of the time integration, the nonlinear iteration and the linear solver.
If itrace < 1itrace<-1, then 1-1 is assumed and similarly if itrace > 3itrace>3, then 33 is assumed.
The advisory messages are given in greater detail as itrace increases. Setting itrace = 1itrace=1 allows you to monitor the progress of the integration without possibly excessive information.
17:   ind – int64int32nag_int scalar
Must be set to 00 or 11.
ind = 0ind=0
Starts the integration in time.
ind = 1ind=1
Continues the integration after an earlier exit from the function. In this case, only the following parameters may be reset between calls to nag_pde_2d_gen_order2_rectilinear (d03rb): tout, dt(2)dt2, dt(3)dt3, tols, tolt, opti, optr, itrace and ifail.
Constraint: 0ind10ind1.
Note:  for users of the NAG Toolbox only, additional values of ind are allowed. See Section [PDEs ()] in the (smpin) for further details.

Optional Input Parameters

1:     npde – int64int32nag_int scalar
Default: The dimension of the array optr.
The number of PDEs in the system.
Constraint: npde1npde1.
2:     lenrwk – int64int32nag_int scalar
Default: The dimension of the array rwk.
The required value of lenrwk cannot be determined exactly in advance, but a suggested value is
lenrwk = maxpts × npde × (5 × l + 18 × npde + 9) + 2 × maxpts,
lenrwk=maxpts×npde×(5×l+18×npde+9)+2×maxpts,
where l = opti(1)l=opti1 if opti(1)0opti10 and l = 3l=3 otherwise, and maxptsmaxpts is the expected maximum number of grid points at any one level. If during the execution the supplied value is found to be too small then the function returns with ifail = 3ifail=3 and an estimated required size is printed on the current error message unit (see nag_file_set_unit_error (x04aa)).
Note:  the size of lenrwk cannot be checked upon initial entry to nag_pde_2d_gen_order2_rectilinear (d03rb) since the number of grid points on the base grid is not known.
3:     leniwk – int64int32nag_int scalar
Default: The dimension of the array iwk.
The dimension of the array iwk as declared in the (sub)program from which nag_pde_2d_gen_order2_rectilinear (d03rb) is called.
The required value of leniwk cannot be determined exactly in advance, but a suggested value is
leniwk = maxpts × (14 + 5 × m) + 7 × m + 2,
leniwk=maxpts×(14+5×m)+7×m+2,
where maxptsmaxpts is the expected maximum number of grid points at any one level and m = opti(1)m=opti1 if opti(1) > 0opti1>0 and m = 3m=3 otherwise. If during the execution the supplied value is found to be too small then the function returns with ifail = 3ifail=3 and an estimated required size is printed on the current error message unit (see nag_file_set_unit_error (x04aa)).
Note:  the size of leniwk cannot be checked upon initial entry to nag_pde_2d_gen_order2_rectilinear (d03rb) since the number of grid points on the base grid is not known.

Input Parameters Omitted from the MATLAB Interface

lwk

Output Parameters

1:     ts – double scalar
The value of tt which has been reached. Normally ts = toutts=tout.
2:     dt(33) – double array
dt(1)dt1 contains the time step size for the next time step. dt(2)dt2 and dt(3)dt3 are unchanged or set to their default values if zero on entry.
3:     rwk(lenrwk) – double array
Communication array, used to store information between calls to nag_pde_2d_gen_order2_rectilinear (d03rb).
4:     iwk(leniwk) – int64int32nag_int array
The following components of the array iwk concern the efficiency of the integration. Here, mm is the maximum number of grid levels allowed (m = opti(1)m=opti1 if opti(1) > 1opti1>1 and m = 3m=3 otherwise), and ll is a grid level taking the values l = 1,2,,nll=1,2,,nl, where nlnl is the number of levels used.
iwk(1)iwk1
Contains the number of steps taken in time.
iwk(2)iwk2
Contains the number of rejected time steps.
iwk(2 + l)iwk2+l
Contains the total number of residual evaluations performed (i.e., the number of times pdedef was called) at grid level ll.
iwk(2 + m + l)iwk2+m+l
Contains the total number of Jacobian evaluations performed at grid level ll.
iwk(2 + 2 × m + l)iwk2+2×m+l
Contains the total number of Newton iterations performed at grid level ll.
iwk(2 + 3 × m + l)iwk2+3×m+l
Contains the total number of linear solver iterations performed at grid level ll.
iwk(2 + 4 × m + l)iwk2+4×m+l
Contains the maximum number of Newton iterations performed at any one time step at grid level ll.
iwk(2 + 5 × m + l)iwk2+5×m+l
Contains the maximum number of linear solver iterations performed at any one time step at grid level ll.
Note:  the total and maximum numbers are cumulative over all calls to nag_pde_2d_gen_order2_rectilinear (d03rb). If the specified maximum number of Newton or linear solver iterations is exceeded at any stage, then the maximums above are set to the specified maximum plus one.
5:     ind – int64int32nag_int scalar
ind = 1ind=1.
Note:  for users of the NAG Toolbox only, additional values of ind are allowed. See Section [PDEs ()] in the (smpin) for further details.
6:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

  ifail = 1ifail=1
On entry,npde < 1npde<1,
ortouttstoutts,
ortout is too close to ts,
orind = 0ind=0 and dt(1) < 0.0dt1<0.0,
or dt(i) < 0.0dti<0.0, for i = 2​ or ​3i=2​ or ​3,
ordt(2) > dt(3)dt2>dt3,
orind = 0ind=0 and 0.0 < dt(1) < 10 × machine precision × max (|ts|,|tout|)0.0<dt1<10×machine precision×max(|ts|,|tout|),
orind = 0ind=0 and dt(1) > touttsdt1>tout-ts,
orind = 0ind=0 and dt(1) < dt(2)​ or ​dt(1) > dt(3)dt1<dt2​ or ​dt1>dt3,
ortols or tolt0.0tolt0.0,
oropti(1) < 0opti1<0,
or opti(1) > 0opti1>0 and opti(j) < 0optij<0, for j = 2j=2, 33 or 44,
oroptr(1,j)0.0optr1j0.0, for some j = 1,2,,npdej=1,2,,npde,
oroptr(2,j) < 0.0optr2j<0.0, for some j = 1,2,,npdej=1,2,,npde,
oroptr(3,j) < 0.0optr3j<0.0, for some j = 1,2,,npdej=1,2,,npde,
orind0ind0 or 11,
orind = 1ind=1 on initial entry to nag_pde_2d_gen_order2_rectilinear (d03rb).
  ifail = 2ifail=2
The time step size to be attempted is less than the specified minimum size. This may occur following time step failures and subsequent step size reductions caused by one or more of the following:
  • the requested accuracy could not be achieved, i.e., tolt is too small,
  • the maximum number of linear solver iterations, Newton iterations or Jacobian evaluations is too small,
  • ILU decomposition of the Jacobian matrix could not be performed, possibly due to singularity of the Jacobian.
Setting itrace to a higher value may provide further information.
In the latter two cases you are advised to check their problem formulation in pdedef and/or bndary, and the initial values in pdeiv if appropriate.
  ifail = 3ifail=3
One or more of the workspace arrays is too small for the required number of grid points. At the initial time step this error may result because you set ierr to 1-1 in inidom or the internal check on the number of grid points following the call to inidom. An estimate of the required sizes for the current stage is output, but more space may be required at a later stage.
W ifail = 4ifail=4
ierr was set to 11 in monitr, forcing control to be passed back to calling program. Integration was successful as far as t = tst=ts.
W ifail = 5ifail=5
The integration has been completed but the maximum number of levels specified in opti(1)opti1 was insufficient at one or more time steps, meaning that the requested space accuracy could not be achieved. To avoid this warning either increase the value of opti(1)opti1 or decrease the value of tols.
  ifail = 6ifail=6
One or more of the output arguments of inidom was incorrectly specified, i.e.,
xminxmaxxminxmax,
orxmax too close to xmin,
oryminymaxyminymax,
orymax too close to ymin,
ornx or ny < 4ny<4,
ornrows < 4nrows<4,
ornrows > nynrows>ny,
ornpts > nx × nynpts>nx×ny,
ornbnds < 8nbnds<8,
ornbpts < 12nbpts<12,
ornbptsnptsnbptsnpts,
orlrow(i) < 1lrowi<1 or lrow(i) > nptslrowi>npts, for some i = 1,2,,nrowsi=1,2,,nrows,
orlrow(i)lrow(i1)lrowilrowi-1, for some i = 2,3,,nrowsi=2,3,,nrows,
orirow(i) < 0irowi<0 or irow(i) > nyirowi>ny, for some i = 1,2,,nrowsi=1,2,,nrows,
orirow(i)irow(i1)irowiirowi-1, for some i = 2,3,,nrowsi=2,3,,nrows,
oricol(i) < 0icoli<0 or icol(i) > nxicoli>nx, for some i = 1,2,,nptsi=1,2,,npts,
orllbnd(i) < 1llbndi<1 or llbnd(i) > nbptsllbndi>nbpts, for some i = 1,2,,nbndsi=1,2,,nbnds,
orllbnd(i)llbnd(i1)llbndillbndi-1, for some i = 2,3,,nbndsi=2,3,,nbnds,
orilbnd(i)1ilbndi1, 22, 33, 44, 1212, 2323, 3434, 4141, 2121, 3232, 4343 or 1414, for some i = 1,2,,nbndsi=1,2,,nbnds,
orlbnd(i) < 1lbndi<1 or lbnd(i) > nptslbndi>npts, for some i = 1,2,,nbptsi=1,2,,nbpts.

Accuracy

There are three sources of error in the algorithm: space and time discretization, and interpolation (linear) between grid levels. The space and time discretization errors are controlled separately using the parameters tols and tolt described in Section [Further Comments], and you should test the effects of varying these parameters. Interpolation errors are generally implicitly controlled by the refinement criterion since in areas where interpolation errors are potentially large, the space monitor will also be large. It can be shown that the global spatial accuracy is comparable to that which would be obtained on a uniform grid of the finest grid size. A full error analysis can be found in Trompert and Verwer (1993).

Further Comments

Algorithm Outline

The local uniform grid refinement method is summarised as follows.
  1. Initialize the course base grid, an initial solution and an initial time step.
  2. Solve the system of PDEs on the current grid with the current time step.
  3. If the required accuracy in space and the maximum number of grid levels have not yet been reached:
    (a) Determine new finer grid at forward time level.
    (b) Get solution values at previous time level(s) on new grid.
    (c) Interpolate internal boundary values from old grid at forward time.
    (d) Get initial values for the Newton process at forward time.
    (e) Go to 22.
  4. Update the coarser grid solution using the finer grid values.
  5. Estimate error in time integration. If time error is acceptable advance time level.
  6. Determine new step size then go to 22 with coarse base as current grid.

Refinement Strategy

For each grid point ii a space monitor μisμis is determined by
μis = max {γj(|Δx2(2)/( x2 )uj(xi,yi,t)| + |Δy2(2)/( y2 )uj(xi,yi,t)|)},
j = 1 , npde
μ i s = max j = 1 , npde { γ j ( | Δ x 2 2 x 2 u j ( x i , y i ,t) | + | Δ y 2 2 y 2 uj ( x i , y i ,t) | ) } ,
(3)
where ΔxΔx and ΔyΔy are the grid widths in the xx and yy directions; and xixi, yiyi are the (x,y)(x,y) coordinates at grid point ii. The parameter γjγj is obtained from
γj = (wjs)/(ujmax σ),
γj=wjs ujmax σ ,
(4)
where σσ is the user-supplied space tolerance; wjswjs is a weighting factor for the relative importance of the jjth PDE component on the space monitor; and ujmaxujmax is the approximate maximum absolute value of the jjth component. A value for σσ must be supplied by you. Values for wjswjs and ujmaxujmax must also be supplied but may be set to the values 1.01.0 if little information about the solution is known.
A new level of refinement is created if
max {μis} > 0.9  or  1.0,
i
maxi{μis}>0.9  or   1.0,
(5)
depending on the grid level at the previous step in order to avoid fluctuations in the number of grid levels between time steps. If (5) is satisfied then all grid points for which μis > 0.25μis>0.25 are flagged and surrounding cells are quartered in size.
No derefinement takes place as such, since at each time step the solution on the base grid is computed first and new finer grids are then created based on the new solution. Hence derefinement occurs implicitly. See Section [Algorithm Outline].

Time Integration

The time integration is controlled using a time monitor calculated at each level ll up to the maximum level used, given by
μlt = sqrt( 1/N j = 1 npde wjt i = 1 ngpts (l) (( Δ t )/(α i j )ut(xi,yi,t))2 )
μ l t = 1 N j = 1 npde w j t i = 1 ngpts ( l ) ( Δ t α i j u t ( x i , y i ,t) ) 2
(6)
where ngpts(l)ngpts(l) is the total number of points on grid level ll; N = ngpts(l) × npdeN=ngpts(l)×npde; ΔtΔt is the current time step; utut is the time derivative of uu which is approximated by first-order finite differences; wjtwjt is the time equivalent of the space weighting factor wjswjs; and αijαij is given by
α i j = τ ((ujmax)/100 + |u(xi,yi,t)|)
α i j = τ ( u j max 100 + | u ( x i , y i ,t) | )
(7)
where ujmaxujmax is as before, and ττ is the user-specified time tolerance.
An integration step is rejected and retried at all levels if
max {μlt} > 1.0.
l
maxl{μlt}>1.0.
(8)

Example

function nag_pde_2d_gen_order2_rectilinear_example
ts = 0;
tout = 0.25;
dt = [0.001;
     1e-07;
     0];
tols = 0.1;
tolt = 0.05;
opti = [int64(5);0;0;0];
optr = [1, 1;
     1, 1;
     1, 1];
rwk = zeros(426000, 1);
iwk = zeros(117037, 1, 'int64');
lenlwk = int64(6000);
itrace = int64(0);
ind = int64(0);
[tsOut, dtOut, ~, ~, indOut, ifail] = ...
    nag_pde_2d_gen_order2_rectilinear(ts, tout, dt, tols, tolt, @inidom, @pdedef, ...
    @bndary, @pdeiv, @monitr, opti, optr, ...
    rwk, iwk, lenlwk, itrace, ind)

function [xmin, xmax, ymin, ymax, nx, ny, npts, nrows, nbnds, nbpts, ...
          lrow, irow, icol, llbnd, ilbnd, lbnd, ierr] = inidom(maxpts, ierr)

  nrows = int64(11);
  npts = int64(105);
  nbnds = int64(28);
  nbpts = int64(72);

  lrow = zeros(nrows, 1, 'int64');
  irow = zeros(nrows, 1, 'int64');
  icol = zeros(npts, 1, 'int64');
  llbnd = zeros(nbnds, 1, 'int64');
  ilbnd = zeros(nbnds, 1, 'int64');
  lbnd = zeros(nbpts, 1, 'int64');



  icold = [int64(0),1,2,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9, ...
           10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,8,9,10,0,1,2,3,4,5, ...
           6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10, ...
           0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8];

  ilbndd = [int64(1),2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,12, ...
           23,34,41,43,14,21,32];

  irowd = [int64(0),1,2,3,4,5,6,7,8,9,10];

  lbndd = [int64(2),4,15,26,37,46,57,68,79,88,98,99,100,101,102,103, ...
           104,96,86,85,84,83,82,70,59,48,39,28,17,6,8,9,10,11,12,13, ...
           18,29,40,49,60,72,73,74,75,76,77,67,56,45,36,25,33,32,42,52, ...
           53,43,1,97,105,87,81,3,7,71,78,14,31,51,54,34];

  llbndd = [int64(1),2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,61,62, ...
           63,64,65,66,67,68,69,70,71,72];

  lrowd = [int64(1),4,15,26,37,46,57,68,79,88,97];

  nx = int64(11);
  ny = int64(11);

% check maxpts against rough estimate of npts
  if (maxpts < double(nx)*double(ny))
     ierr = -1;
     return;
  end

  xmin = 0.0d0;
  ymin = 0.0d0;
  xmax = 1.0d0;
  ymax = 1.0d0;

  for i = 1:double(nrows)
    lrow(i) = lrowd(i);
    irow(i) = irowd(i);
  end

  for i = 1:double(nbnds)
    llbnd(i) = llbndd(i);
    ilbnd(i) = ilbndd(i);
  end

  for i = 1:double(nbpts)
    lbnd(i) = lbndd(i);
  end

  for i = 1:double(npts)
    icol(i) = icold(i);
  end

function [res] = pdedef(npts, npde, t, x, y, u, ut, ux, uy, uxx, uxy, uyy)
  res = zeros(npts, npde);

  epsilon = 1e-3;
  for i = 1:double(npts)
    res(i,1) = ut(i,1)-(-u(i,1)*ux(i,1)-u(i,2)*uy(i,1)+epsilon*(uxx(i,1)+uyy(i,1)));
    res(i,2) = ut(i,2)-(-u(i,1)*ux(i,2)-u(i,2)*uy(i,2)+epsilon*(uxx(i,2)+uyy(i,2)));
  end

function [res] = ...
    bndary(npts, npde, t, x, y, u, ut, ux, uy, nbnds, nbpts, llbnd, ilbnd, lbnd, res)

      epsilon = 1e-3;

      for k = double(llbnd(1)):double(nbpts)
         i = lbnd(k);
         a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
         if (a <= 0.0d0)
            res(i,1) = u(i,1) - (0.75d0-0.25d0/(1.0d0+exp(a)));
            res(i,2) = u(i,2) - (0.75d0+0.25d0/(1.0d0+exp(a)));
         else
            res(i,1) = u(i,1) - (0.75d0-0.25d0*exp(-a)/(exp(-a)+1.0d0));
            res(i,2) = u(i,2) - (0.75d0+0.25d0*exp(-a)/(exp(-a)+1.0d0));
         end
      end

function [u] = pdeiv(npts, npde, t, x, y)
  u = zeros(npts, npde);

  epsilon = 1e-3;

  for i = 1:double(npts)
    a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
    if (a <= 0.0d0)
       u(i,1) = 0.75d0 - 0.25d0/(1.0d0+exp(a));
       u(i,2) = 0.75d0 + 0.25d0/(1.0d0+exp(a));
    else
       u(i,1) = 0.75d0 - 0.25d0*exp(-a)/(exp(-a)+1.0d0);
       u(i,2) = 0.75d0 + 0.25d0*exp(-a)/(exp(-a)+1.0d0);
    end
  end

function [ierr] = ...
    monitr(npde, t, dt, dtnew, tlast, nlev, xmin, ymin, dxb, dyb, lgrid, ...
                 istruc, lsol, sol, ierr)
 

tsOut =

    0.2500


dtOut =

    0.0249
    0.0000
    0.2500


indOut =

                    1


ifail =

                    0


function d03rb_example
ts = 0;
tout = 0.25;
dt = [0.001;
     1e-07;
     0];
tols = 0.1;
tolt = 0.05;
opti = [int64(5);0;0;0];
optr = [1, 1;
     1, 1;
     1, 1];
rwk = zeros(426000, 1);
iwk = zeros(117037, 1, 'int64');
lenlwk = int64(6000);
itrace = int64(0);
ind = int64(0);
[tsOut, dtOut, ~, ~, indOut, ifail] = ...
    d03rb(ts, tout, dt, tols, tolt, @inidom, @pdedef, ...
    @bndary, @pdeiv, @monitr, opti, optr, ...
    rwk, iwk, lenlwk, itrace, ind)

function [xmin, xmax, ymin, ymax, nx, ny, npts, nrows, nbnds, nbpts, ...
          lrow, irow, icol, llbnd, ilbnd, lbnd, ierr] = inidom(maxpts, ierr)

  nrows = int64(11);
  npts = int64(105);
  nbnds = int64(28);
  nbpts = int64(72);

  lrow = zeros(nrows, 1, 'int64');
  irow = zeros(nrows, 1, 'int64');
  icol = zeros(npts, 1, 'int64');
  llbnd = zeros(nbnds, 1, 'int64');
  ilbnd = zeros(nbnds, 1, 'int64');
  lbnd = zeros(nbpts, 1, 'int64');



  icold = [int64(0),1,2,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9, ...
           10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,8,9,10,0,1,2,3,4,5, ...
           6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10, ...
           0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8];

  ilbndd = [int64(1),2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,12, ...
           23,34,41,43,14,21,32];

  irowd = [int64(0),1,2,3,4,5,6,7,8,9,10];

  lbndd = [int64(2),4,15,26,37,46,57,68,79,88,98,99,100,101,102,103, ...
           104,96,86,85,84,83,82,70,59,48,39,28,17,6,8,9,10,11,12,13, ...
           18,29,40,49,60,72,73,74,75,76,77,67,56,45,36,25,33,32,42,52, ...
           53,43,1,97,105,87,81,3,7,71,78,14,31,51,54,34];

  llbndd = [int64(1),2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,61,62, ...
           63,64,65,66,67,68,69,70,71,72];

  lrowd = [int64(1),4,15,26,37,46,57,68,79,88,97];

  nx = int64(11);
  ny = int64(11);

% check maxpts against rough estimate of npts
  if (maxpts < double(nx)*double(ny))
     ierr = -1;
     return;
  end

  xmin = 0.0d0;
  ymin = 0.0d0;
  xmax = 1.0d0;
  ymax = 1.0d0;

  for i = 1:double(nrows)
    lrow(i) = lrowd(i);
    irow(i) = irowd(i);
  end

  for i = 1:double(nbnds)
    llbnd(i) = llbndd(i);
    ilbnd(i) = ilbndd(i);
  end

  for i = 1:double(nbpts)
    lbnd(i) = lbndd(i);
  end

  for i = 1:double(npts)
    icol(i) = icold(i);
  end

function [res] = pdedef(npts, npde, t, x, y, u, ut, ux, uy, uxx, uxy, uyy)
  res = zeros(npts, npde);

  epsilon = 1e-3;
  for i = 1:double(npts)
    res(i,1) = ut(i,1)-(-u(i,1)*ux(i,1)-u(i,2)*uy(i,1)+epsilon*(uxx(i,1)+uyy(i,1)));
    res(i,2) = ut(i,2)-(-u(i,1)*ux(i,2)-u(i,2)*uy(i,2)+epsilon*(uxx(i,2)+uyy(i,2)));
  end

function [res] = ...
    bndary(npts, npde, t, x, y, u, ut, ux, uy, nbnds, nbpts, llbnd, ilbnd, lbnd, res)

      epsilon = 1e-3;

      for k = double(llbnd(1)):double(nbpts)
         i = lbnd(k);
         a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
         if (a <= 0.0d0)
            res(i,1) = u(i,1) - (0.75d0-0.25d0/(1.0d0+exp(a)));
            res(i,2) = u(i,2) - (0.75d0+0.25d0/(1.0d0+exp(a)));
         else
            res(i,1) = u(i,1) - (0.75d0-0.25d0*exp(-a)/(exp(-a)+1.0d0));
            res(i,2) = u(i,2) - (0.75d0+0.25d0*exp(-a)/(exp(-a)+1.0d0));
         end
      end

function [u] = pdeiv(npts, npde, t, x, y)
  u = zeros(npts, npde);

  epsilon = 1e-3;

  for i = 1:double(npts)
    a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
    if (a <= 0.0d0)
       u(i,1) = 0.75d0 - 0.25d0/(1.0d0+exp(a));
       u(i,2) = 0.75d0 + 0.25d0/(1.0d0+exp(a));
    else
       u(i,1) = 0.75d0 - 0.25d0*exp(-a)/(exp(-a)+1.0d0);
       u(i,2) = 0.75d0 + 0.25d0*exp(-a)/(exp(-a)+1.0d0);
    end
  end

function [ierr] = ...
    monitr(npde, t, dt, dtnew, tlast, nlev, xmin, ymin, dxb, dyb, lgrid, ...
                 istruc, lsol, sol, ierr)
 

tsOut =

    0.2500


dtOut =

    0.0249
    0.0000
    0.2500


indOut =

                    1


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013