NAG Toolbox: nag_pde_2d_gen_order2_rectangle (d03ra)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{ts}$ – double scalar

The initial value of the independent variable $t$.

Constraint: $\mathbf{ts}<\mathbf{tout}$.

2:     $\mathrm{tout}$ – double scalar

The final value of $t$ to which the integration is to be carried out.

3:     $\mathrm{dt}\left(3\right)$ – double array

The initial, minimum and maximum time step sizes respectively.

$\mathbf{dt}\left(1\right)$

Specifies the initial time step size to be used on the first entry, i.e., when $\mathbf{ind}=0$. If $\mathbf{dt}\left(1\right)=0.0$ then the default value $\mathbf{dt}\left(1\right)=0.01\times \left(\mathbf{tout}-\mathbf{ts}\right)$ is used. On subsequent entries ($\mathbf{ind}=1$), the value of $\mathbf{dt}\left(1\right)$ is not referenced.

$\mathbf{dt}\left(2\right)$

Specifies the minimum time step size to be attempted by the integrator. If $\mathbf{dt}\left(2\right)=0.0$ the default value $\mathbf{dt}\left(2\right)=10.0\times \mathit{machine precision}$ is used.

$\mathbf{dt}\left(3\right)$

Specifies the maximum time step size to be attempted by the integrator. If $\mathbf{dt}\left(3\right)=0.0$ the default value $\mathbf{dt}\left(3\right)=\mathbf{tout}-\mathbf{ts}$ is used.

Constraints:

• if $\mathbf{ind}=0$, $\mathbf{dt}\left(1\right)\ge 0.0$;
• if $\mathbf{ind}=0$ and $\mathbf{dt}\left(1\right)>0.0$,
  $10.0\times \mathit{machine precision}\times \max \left(\left|\mathbf{ts}\right|,\left|\mathbf{tout}\right|\right)\le \mathbf{dt}\left(1\right)\le \mathbf{tout}-\mathbf{ts}$ and
  $\mathbf{dt}\left(2\right)\le \mathbf{dt}\left(1\right)\le \mathbf{dt}\left(3\right)$, where the values of $\mathbf{dt}\left(2\right)$ and $\mathbf{dt}\left(3\right)$ will have been reset to their default values if zero on entry;
• $0\le \mathbf{dt}\left(2\right)\le \mathbf{dt}\left(3\right)$.

4:     $\mathrm{xmin}$ – double scalar

5:     $\mathrm{xmax}$ – double scalar

The extents of the rectangular domain in the $x$-direction, i.e., the $x$ coordinates of the left and right boundaries respectively.

Constraint: $\mathbf{xmin}<\mathbf{xmax}$ and xmax must be sufficiently distinguishable from xmin for the precision of the machine being used.

6:     $\mathrm{ymin}$ – double scalar

7:     $\mathrm{ymax}$ – double scalar

The extents of the rectangular domain in the $y$-direction, i.e., the $y$ coordinates of the lower and upper boundaries respectively.

Constraint: $\mathbf{ymin}<\mathbf{ymax}$ and ymax must be sufficiently distinguishable from ymin for the precision of the machine being used.

8:     $\mathrm{nx}$ – int64int32nag_int scalar

The number of grid points in the $x$-direction (including the boundary points).

Constraint: $\mathbf{nx}\ge 4$.

9:     $\mathrm{ny}$ – int64int32nag_int scalar

The number of grid points in the $y$-direction (including the boundary points).

Constraint: $\mathbf{ny}\ge 4$.

10:   $\mathrm{tols}$ – double scalar

The space tolerance used in the grid refinement strategy ($\sigma$ in equation (4)). See Refinement Strategy.

Constraint: $\mathbf{tols}>0.0$.

11:   $\mathrm{tolt}$ – double scalar

The time tolerance used to determine the time step size ($\tau$ in equation (7)). See Time Integration .

Constraint: $\mathbf{tolt}>0.0$.

12:   $\mathrm{pdedef}$ – function handle or string containing name of m-file

pdedef must evaluate the functions $F_j$, for $\mathit{j}=1,2,\dots ,\mathbf{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:     $\mathrm{npts}$ – int64int32nag_int scalar

The number of grid points in the current grid.

2:     $\mathrm{npde}$ – int64int32nag_int scalar

The number of PDEs in the system.

3:     $\mathrm{t}$ – double scalar

The current value of the independent variable $t$.

4:     $\mathrm{x}\left(\mathbf{npts}\right)$ – double array

$\mathbf{x}\left(\mathit{i}\right)$ contains the $x$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$.

5:     $\mathrm{y}\left(\mathbf{npts}\right)$ – double array

$\mathbf{y}\left(\mathit{i}\right)$ contains the $y$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$.

6:     $\mathrm{u}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{u}\left(\mathit{i},\mathit{j}\right)$ contains the value of the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

7:     $\mathrm{ut}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{ut}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial t}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

8:     $\mathrm{ux}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{ux}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial x}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

9:     $\mathrm{uy}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{uy}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial y}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

10:   $\mathrm{uxx}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{uxx}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{{\partial }^{2}u}{\partial x^2}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

11:   $\mathrm{uxy}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{uxy}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{{\partial }^{2}u}{\partial x\partial y}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

12:   $\mathrm{uyy}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{uyy}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{{\partial }^{2}u}{\partial y^2}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

Output Parameters

1:     $\mathrm{res}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{res}\left(\mathit{i},\mathit{j}\right)$ must contain the value of $F_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$, at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$, although the residuals at boundary points will be ignored (and overwritten later on) and so they need not be specified here.

13:   $\mathrm{bndary}$ – function handle or string containing name of m-file

bndary must evaluate the functions $G_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,\mathbf{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, nbpts, lbnd, res)

Input Parameters

1:     $\mathrm{npts}$ – int64int32nag_int scalar

The number of grid points in the current grid.

2:     $\mathrm{npde}$ – int64int32nag_int scalar

The number of PDEs in the system.

3:     $\mathrm{t}$ – double scalar

The current value of the independent variable $t$.

4:     $\mathrm{x}\left(\mathbf{npts}\right)$ – double array

$\mathbf{x}\left(\mathit{i}\right)$ contains the $x$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$.

5:     $\mathrm{y}\left(\mathbf{npts}\right)$ – double array

$\mathbf{y}\left(\mathit{i}\right)$ contains the $y$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$.

6:     $\mathrm{u}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{u}\left(\mathit{i},\mathit{j}\right)$ contains the value of the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

7:     $\mathrm{ut}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{ut}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial t}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

8:     $\mathrm{ux}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{ux}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial x}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

9:     $\mathrm{uy}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{uy}\left(\mathit{i},\mathit{j}\right)$ contains the value of $\frac{\partial u}{\partial y}$ for the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

10:   $\mathrm{nbpts}$ – int64int32nag_int scalar

The number of boundary points in the grid.

11:   $\mathrm{lbnd}\left(\mathbf{nbpts}\right)$ – int64int32nag_int array

$\mathbf{lbnd}\left(\mathit{i}\right)$ contains the grid index for the $\mathit{i}$th boundary point, for $\mathit{i}=1,2,\dots ,\mathbf{nbpts}$. Hence the $\mathit{i}$th boundary point has coordinates $\mathbf{x}\left(\mathbf{lbnd}\left(\mathit{i}\right)\right)$ and $\mathbf{y}\left(\mathbf{lbnd}\left(\mathit{i}\right)\right)$, and the corresponding solution values are $\mathbf{u}\left(\mathbf{lbnd}\left(\mathit{i}\right),\mathbf{npde}\right)$, etc.

12:   $\mathrm{res}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{res}\left(\mathit{i},\mathit{j}\right)$ contains the value of $F_{\mathit{j}}$, for $\mathit{i}=1,2,\dots ,\mathbf{npde}$, at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$, as returned by pdedef. The residuals at the boundary points will be overwritten and so need not have been set by pdedef.

Output Parameters

1:     $\mathrm{res}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{res}\left(\mathbf{lbnd}\left(\mathit{i}\right),\mathit{j}\right)$ must contain the value of $G_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$, at the $\mathit{i}$th boundary point, for $\mathit{i}=1,2,\dots ,\mathbf{nbpts}$.

Note: elements of res corresponding to interior points must not be altered.

14:   $\mathrm{pdeiv}$ – function handle or string containing name of m-file

pdeiv must specify the initial values of the PDE components $u$ at all points in the grid. pdeiv is not referenced if, on entry, $\mathbf{ind}=1$.

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

Input Parameters

1:     $\mathrm{npts}$ – int64int32nag_int scalar

The number of grid points in the grid.

2:     $\mathrm{npde}$ – int64int32nag_int scalar

The number of PDEs in the system.

3:     $\mathrm{t}$ – double scalar

The (initial) value of the independent variable $t$.

4:     $\mathrm{x}\left(\mathbf{npts}\right)$ – double array

$\mathbf{x}\left(\mathit{i}\right)$ contains the $x$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$.

5:     $\mathrm{y}\left(\mathbf{npts}\right)$ – double array

$\mathbf{y}\left(\mathit{i}\right)$ contains the $y$ coordinate of the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$.

Output Parameters

1:     $\mathrm{u}\left(\mathbf{npts},\mathbf{npde}\right)$ – double array

$\mathbf{u}\left(\mathit{i},\mathit{j}\right)$ must contain the value of the $\mathit{j}$th PDE component at the $\mathit{i}$th grid point, for $\mathit{i}=1,2,\dots ,\mathbf{npts}$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

15:   $\mathrm{monitr}$ – function handle or string containing name of m-file

monitr is called by nag_pde_2d_gen_order2_rectangle (d03ra) 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 argument tlast. The input arguments contain information about the grid and solution at all grid levels used.

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, ngpts, xpts, ypts, lsol, sol, ierr)

Input Parameters

1:     $\mathrm{npde}$ – int64int32nag_int scalar

The number of PDEs in the system.

2:     $\mathrm{t}$ – double scalar

The current value of the independent variable $t$, i.e., the time at the end of the integration step just completed.

3:     $\mathrm{dt}$ – double scalar

The current time step size $\Delta t$, i.e., the time step size used for the integration step just completed.

4:     $\mathrm{dtnew}$ – double scalar

The step size that will be used for the next time step.

5:     $\mathrm{tlast}$ – logical scalar

Indicates if intermediate or final time step. $\mathbf{tlast}=\mathit{false}$ for an intermediate step, $\mathbf{tlast}=\mathit{true}$ for the last call to monitr before returning to your program.

6:     $\mathrm{nlev}$ – int64int32nag_int scalar

The number of grid levels used at time t.

7:     $\mathrm{ngpts}\left(\mathbf{nlev}\right)$ – int64int32nag_int array

$\mathbf{ngpts}\left(\mathit{l}\right)$ contains the number of grid points at level $\mathit{l}$, for $\mathit{l}=1,2,\dots ,\mathbf{nlev}$.

8:     $\mathrm{xpts}\left(\mathit{lpts}\right)$ – double array

Contains the $x$ coordinates of the grid points in each level in turn, i.e., $\mathbf{x}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{ngpts}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,\mathbf{nlev}$.

So for level $\mathit{l}$, $\mathbf{x}\left(\mathit{i}\right)=\mathbf{xpts}\left(\mathit{k}+\mathit{i}\right)$, where $\mathit{k}=\mathbf{ngpts}\left(1\right)+\mathbf{ngpts}\left(2\right)+\cdots +\mathbf{ngpts}\left(\mathit{l}-1\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{ngpts}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,\mathbf{nlev}$.

9:     $\mathrm{ypts}\left(\mathit{lpts}\right)$ – double array

Contains the $y$ coordinates of the grid points in each level in turn, i.e., $\mathbf{y}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{ngpts}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,\mathbf{nlev}$.

So for level $\mathit{l}$, $\mathbf{y}\left(\mathit{i}\right)=\mathbf{ypts}\left(\mathit{k}+\mathit{i}\right)$, where $\mathit{k}=\mathbf{ngpts}\left(1\right)+\mathbf{ngpts}\left(2\right)+\cdots +\mathbf{ngpts}\left(\mathit{l}-1\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{ngpts}\left(\mathit{l}\right)$ and $\mathit{l}=1,2,\dots ,\mathbf{nlev}$.

10:   $\mathrm{lsol}\left(\mathbf{nlev}\right)$ – int64int32nag_int array

$\mathbf{lsol}\left(\mathit{l}\right)$ contains the pointer to the solution in sol at grid level $\mathit{l}$ and time t. ($\mathbf{lsol}\left(\mathit{l}\right)$ actually contains the array index immediately preceding the start of the solution in sol.)

11:   $\mathrm{sol}\left(:\right)$ – double array

Default: $\mathit{lenrwk}-6\times \mathbf{npde}$

Contains the solution $\mathbf{u}\left(\mathbf{ngpts}\left(\mathit{l}\right),\mathbf{npde}\right)$ at time t for each grid level $\mathit{l}$ in turn, positioned according to lsol, i.e., for level $\mathit{l}$, $\mathbf{u}\left(\mathit{i},\mathit{j}\right)=\mathbf{sol}\left(\mathbf{lsol}\left(\mathit{l}\right)+\left(\mathit{j}-1\right)\times \mathbf{ngpts}\left(\mathit{l}\right)+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{ngpts}\left(\mathit{l}\right)$, $\mathit{j}=1,2,\dots ,\mathbf{npde}$ and $\mathit{l}=1,2,\dots ,\mathbf{nlev}$.

12:   $\mathrm{ierr}$ – int64int32nag_int scalar

Will be set to $0$.

Output Parameters

1:     $\mathrm{ierr}$ – int64int32nag_int scalar

Should be set to $1$ to force a tidy termination and an immediate return to the calling program with $\mathbf{ifail}=\mathbf{4}$. ierr should remain unchanged otherwise.

16:   $\mathrm{opti}\left(4\right)$ – int64int32nag_int array

May be set to control various options available in the integrator.

$\mathbf{opti}\left(1\right)=0$

All the default options are employed.

$\mathbf{opti}\left(1\right)>0$

The default value of $\mathbf{opti}\left(\mathit{i}\right)$, for $\mathit{i}=2,3,4$, can be obtained by setting $\mathbf{opti}\left(\mathit{i}\right)=0$.

$\mathbf{opti}\left(1\right)$

Specifies the maximum number of grid levels allowed (including the base grid). $\mathbf{opti}\left(1\right)\ge 0$. The default value is $\mathbf{opti}\left(1\right)=3$.

$\mathbf{opti}\left(2\right)$

Specifies the maximum number of Jacobian evaluations allowed during each nonlinear equations solution. $\mathbf{opti}\left(2\right)\ge 0$. The default value is $\mathbf{opti}\left(2\right)=2$.

$\mathbf{opti}\left(3\right)$

Specifies the maximum number of Newton iterations in each nonlinear equations solution. $\mathbf{opti}\left(3\right)\ge 0$. The default value is $\mathbf{opti}\left(3\right)=10$.

$\mathbf{opti}\left(4\right)$

Specifies the maximum number of iterations in each linear equations solution. $\mathbf{opti}\left(4\right)\ge 0$. The default value is $\mathbf{opti}\left(4\right)=100$.

Constraint: $\mathbf{opti}\left(1\right)\ge 0$ and if $\mathbf{opti}\left(1\right)>0$, $\mathbf{opti}\left(\mathit{i}\right)\ge 0$, for $\mathit{i}=2,3,4$.

17:   $\mathrm{optr}\left(3,\mathbf{npde}\right)$ – double array

May be used to specify the optional vectors $u^{\max }$, $w^{\mathrm{s}}$ and $w^{\mathrm{t}}$ in the space and time monitors (see Further Comments).

If an optional vector is not required then all its components should be set to $1.0$.

$\mathbf{optr}\left(1,\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$, specifies $u_{\mathit{j}}^{\max }$, the approximate maximum absolute value of the $\mathit{j}$th component of $u$, as used in (4) and (7). $\mathbf{optr}\left(1,\mathit{j}\right)>0.0$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

$\mathbf{optr}\left(2,\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$, specifies $w_{\mathit{j}}^{\mathrm{s}}$, the weighting factors used in the space monitor (see (4)) to indicate the relative importance of the $\mathit{j}$th component of $u$ on the space monitor. $\mathbf{optr}\left(2,\mathit{j}\right)\ge 0.0$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

$\mathbf{optr}\left(3,\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$, specifies $w_{\mathit{j}}^{\mathrm{t}}$, the weighting factors used in the time monitor (see (6)) to indicate the relative importance of the $\mathit{j}$th component of $u$ on the time monitor. $\mathbf{optr}\left(3,\mathit{j}\right)\ge 0.0$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

Constraints:

• $\mathbf{optr}\left(1,\mathit{j}\right)>0.0$, for $\mathit{j}=1,2,\dots ,\mathbf{npde}$;
• $\mathbf{optr}\left(\mathit{i},\mathit{j}\right)\ge 0.0$, for $\mathit{i}=2,3$ and $\mathit{j}=1,2,\dots ,\mathbf{npde}$.

18:   $\mathrm{rwk}\left(\mathit{lenrwk}\right)$ – double array

lenrwk, the dimension of the array, must satisfy the constraint $\mathit{lenrwk}\ge \mathbf{nx}\times \mathbf{ny}\times \mathbf{npde}\times \left(14+18\times \mathbf{npde}\right)+2\times \mathbf{nx}\times \mathbf{ny}$ (the required size for the initial grid).

The required value of lenrwk cannot be determined exactly in advance, but a suggested value is

$$\mathit{lenrwk}=\mathit{maxpts}\times \mathbf{npde}\times \left(5\times \mathit{l}+18\times \mathbf{npde}+9\right)+2\times \mathit{maxpts}\text{,}$$

where $\mathit{l}=\mathbf{opti}\left(1\right)$ if $\mathbf{opti}\left(1\right)\ne 0$ and $\mathit{l}=3$ otherwise, and $\mathit{maxpts}$ 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 $\mathbf{ifail}=\mathbf{3}$ and an estimated required size is printed on the current error message unit (see nag_file_set_unit_error (x04aa)).

Constraint: $\mathit{lenrwk}\ge \mathbf{nx}\times \mathbf{ny}\times \mathbf{npde}\times \left(14+18\times \mathbf{npde}\right)+2\times \mathbf{nx}\times \mathbf{ny}$ (the required size for the initial grid).

19:   $\mathrm{iwk}\left(\mathbf{leniwk}\right)$ – int64int32nag_int array

If $\mathbf{ind}=0$, iwk need not be set. Otherwise iwk must remain unchanged from a previous call to nag_pde_2d_gen_order2_rectangle (d03ra).

20:   $\mathrm{itrace}$ – int64int32nag_int scalar

The level of trace information required from nag_pde_2d_gen_order2_rectangle (d03ra). itrace may take the value $-1$, $0$, $1$, $2$ or $3$.

$\mathbf{itrace}=-1$

No output is generated.

$\mathbf{itrace}=0$

Only warning messages are printed.

$\mathbf{itrace}>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 $\mathbf{itrace}<-1$, then $-1$ is assumed and similarly if $\mathbf{itrace}>3$, then $3$ is assumed.

The advisory messages are given in greater detail as itrace increases. Setting $\mathbf{itrace}=1$ allows you to monitor the progress of the integration without possibly excessive information.

21:   $\mathrm{ind}$ – int64int32nag_int scalar

Must be set to $0$ or $1$, alternatively $10$ or $11$.

$\mathbf{ind}=0$

Starts the integration in time. pdedef is assumed to be serial.

$\mathbf{ind}=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_rectangle (d03ra): tout, dt, tols, tolt, opti, optr, itrace and ifail. pdedef is assumed to be serial.

$\mathbf{ind}=10$

Equivalent to $\mathbf{ind}=0$. This option is included only for compatibility with other NAG library products.

$\mathbf{ind}=11$

Equivalent to $\mathbf{ind}=1$. This option is included only for compatibility with other NAG library products.

Constraint: $0\le \mathbf{ind}\le 1$ or $10\le \mathbf{ind}\le 11$.

Optional Input Parameters

1:     $\mathrm{npde}$ – int64int32nag_int scalar

Default: the dimension of the array optr.

The number of PDEs in the system.

Constraint: $\mathbf{npde}\ge 1$.

2:     $\mathrm{leniwk}$ – int64int32nag_int scalar

Default: the dimension of the array iwk.

The dimension of the array iwk.

The required value of leniwk cannot be determined exactly in advance, but a suggested value is

$$\mathbf{leniwk}=\mathit{maxpts}\times \left(14+5\times m\right)+7\times m+2\text{,}$$

where $\mathit{maxpts}$ is the expected maximum number of grid points at any one level and $m=\mathbf{opti}\left(1\right)$ if $\mathbf{opti}\left(1\right)>0$ and $m=3$ otherwise. If during the execution the supplied value is found to be too small then the function returns with $\mathbf{ifail}=\mathbf{3}$ and an estimated required size is printed on the current error message unit (see nag_file_set_unit_error (x04aa)).

Constraint: $\mathbf{leniwk}\ge 19\times \mathbf{nx}\times \mathbf{ny}+9$ (the required size for the initial grid).

3:     $\mathrm{lenlwk}$ – int64int32nag_int scalar

Default: $\mathit{maxpts}+1$ 

The dimension of the array lwk.

The required value of lenlwk cannot be determined exactly in advanced, but a suggested value is

$$\mathbf{lenlwk}=\mathit{maxpts}+1\text{,}$$

where $\mathit{maxpts}$ 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 $\mathbf{ifail}=\mathbf{3}$ and an estimated required size is printed on the current error message unit (see nag_file_set_unit_error (x04aa)).

Constraint: $\mathbf{lenlwk}\ge \mathbf{nx}\times \mathbf{ny}+1$ (the required size for the initial grid).

Output Parameters

1:     $\mathrm{ts}$ – double scalar

The value of $t$ which has been reached. Normally $\mathbf{ts}=\mathbf{tout}$.

2:     $\mathrm{dt}\left(3\right)$ – double array

$\mathbf{dt}\left(1\right)$ contains the time step size for the next time step. $\mathbf{dt}\left(2\right)$ and $\mathbf{dt}\left(3\right)$ are unchanged or set to their default values if zero on entry.

3:     $\mathrm{rwk}\left(\mathit{lenrwk}\right)$ – double array

$\mathit{lenrwk}=\mathbf{nx}\times \mathbf{ny}\times \mathbf{npde}\times \left(14+18\times \mathbf{npde}\right)+2\times \mathbf{nx}\times \mathbf{ny}$ (the required size for the initial grid).

Communication array, used to store information between calls to nag_pde_2d_gen_order2_rectangle (d03ra).

4:     $\mathrm{iwk}\left(\mathbf{leniwk}\right)$ – int64int32nag_int array

The following components of the array iwk concern the efficiency of the integration. Here, $m$ is the maximum number of grid levels allowed ($m=\mathbf{opti}\left(1\right)$ if $\mathbf{opti}\left(1\right)>1$ and $m=3$ otherwise), and $\mathit{l}$ is a grid level taking the values $\mathit{l}=1,2,\dots ,\mathit{nl}$, where $\mathit{nl}$ is the number of levels used.

$\mathbf{iwk}\left(1\right)$

Contains the number of steps taken in time.

$\mathbf{iwk}\left(2\right)$

Contains the number of rejected time steps.

$\mathbf{iwk}\left(2+\mathit{l}\right)$

Contains the total number of residual evaluations performed (i.e., the number of times pdedef was called) at grid level $\mathit{l}$.

$\mathbf{iwk}\left(2+m+\mathit{l}\right)$

Contains the total number of Jacobian evaluations performed at grid level $\mathit{l}$.

$\mathbf{iwk}\left(2+2\times m+\mathit{l}\right)$

Contains the total number of Newton iterations performed at grid level $\mathit{l}$.

$\mathbf{iwk}\left(2+3\times m+\mathit{l}\right)$

Contains the total number of linear solver iterations performed at grid level $\mathit{l}$.

$\mathbf{iwk}\left(2+4\times m+\mathit{l}\right)$

Contains the maximum number of Newton iterations performed at any one time step at grid level $\mathit{l}$.

$\mathbf{iwk}\left(2+5\times m+\mathit{l}\right)$

Contains the maximum number of linear solver iterations performed at any one time step at grid level $\mathit{l}$.

Note: the total and maximum numbers are cumulative over all calls to nag_pde_2d_gen_order2_rectangle (d03ra). 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:     $\mathrm{ind}$ – int64int32nag_int scalar

$\mathbf{ind}=1$, if ind on input was $0$ or $1$, or $\mathbf{ind}=11$, if ind on input was $10$ or $11$.

6:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

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

   $\mathbf{ifail}=1$

On entry,	$\mathbf{npde}<1$,
or	$\mathbf{tout}\le \mathbf{ts}$,
or	tout is too close to ts,
or	$\mathbf{ind}=0$ and $\mathbf{dt}\left(1\right)<0.0$,
or	$\mathbf{dt}\left(i\right)<0.0$, for $\mathit{i}=2\text{ or }3$,
or	$\mathbf{dt}\left(2\right)>\mathbf{dt}\left(3\right)$,
or	$\mathbf{ind}=0$ and $0.0<\mathbf{dt}\left(1\right)<10\times \mathit{machine precision}\times \max \left(\left|\mathbf{ts}\right|,\left|\mathbf{tout}\right|\right)$,
or	$\mathbf{ind}=0$ and $\mathbf{dt}\left(1\right)>\mathbf{tout}-\mathbf{ts}$,
or	$\mathbf{ind}=0$ and $\mathbf{dt}\left(1\right)<\mathbf{dt}\left(2\right)$ or $\mathbf{dt}\left(1\right)>\mathbf{dt}\left(3\right)$,
or	$\mathbf{xmin}\ge \mathbf{xmax}$,
or	xmax too close to xmin,
or	$\mathbf{ymin}\ge \mathbf{ymax}$,
or	ymax too close to ymin,
or	nx or $\mathbf{ny}<4$,
or	tols or $\mathbf{tolt}\le 0.0$,
or	$\mathbf{opti}\left(1\right)<0$,
or	$\mathbf{opti}\left(1\right)>0$ and $\mathbf{opti}\left(j\right)<0$, for $\mathit{j}=2$, $3$ or $4$,
or	$\mathbf{optr}\left(1,j\right)\le 0.0$, for some $j=1,2,\dots ,\mathbf{npde}$,
or	$\mathbf{optr}\left(2,j\right)<0.0$, for some $j=1,2,\dots ,\mathbf{npde}$,
or	$\mathbf{optr}\left(3,j\right)<0.0$, for some $j=1,2,\dots ,\mathbf{npde}$,
or	lenrwk, leniwk or lenlwk too small for initial grid level,
or	$\mathbf{ind}\ne 0$ or $1$,
or	$\mathbf{ind}=1$ on initial entry to nag_pde_2d_gen_order2_rectangle (d03ra).

   $\mathbf{ifail}=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.

   $\mathbf{ifail}=3$

One or more of the workspace arrays is too small for the required number of grid points. An estimate of the required sizes for the current stage is output, but more space may be required at a later stage.

W  $\mathbf{ifail}=4$

ierr was set to $1$ in monitr, forcing control to be passed back to calling program. Integration was successful as far as $\mathbf{t}=\mathbf{ts}$.

W  $\mathbf{ifail}=5$

The integration has been completed but the maximum number of levels specified in $\mathbf{opti}\left(1\right)$ 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 $\mathbf{opti}\left(1\right)$ or decrease the value of tols.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Algorithm Outline

Refinement Strategy

Time Integration

Example

Open in the MATLAB editor: d03ra_example

function d03ra_example


fprintf('d03ra example results\n\n');

global ssav1 ssav2;
global xsav ysav tsav;

ex1;
fig1 = figure;
ex1_plot;


ex2;
fig2 = figure;
ex2_plot(1);
fig3 = figure;
ex2_plot(2);



function ex1
% First d03ra example.  Temperature distribution for a single, one-step
% reaction of a mixture of two chemicals.
global heat_release activ_energy diffusion damkohler;
global iout isav;

% Set values for problem parameters.
mxlev = 3;
npde = 1;
npts = 2000;
xmin = 0;    xmax = 1;
ymin = 0;    ymax = 1;
tols = 0.5;
tolt = 0.01;
nx   = int64(21);
ny   = nx;

reaction_rate =  5;
activ_energy  = 20;
heat_release  =  1;
diffusion     =  0.1;
damkohler = reaction_rate*exp(activ_energy)/(heat_release*activ_energy);

leniwk = npts*(14 + 5*mxlev) + 7*mxlev + 2;
lenlwk = npts + 1;
lenrwk = npts*npde*(5*mxlev + 18*npde + 9) + 2*npts;

% The above are the recommended values from the documentation, but d03ra
% returns with ifail = 3 and the minumum workspace size requirements if
% these are used.  So we beef them up here with some fudge factors.
leniwk = 2*leniwk;
lenlwk = int64(2*lenlwk);
lenrwk = 2*lenrwk;

% Create some work arrays - see NAG documentation.
rwk = zeros(lenrwk, 1);
iwk = zeros(leniwk, 1, 'int64');

% Initialize some problem variables.
ind    = int64(0);
itrace = int64(0);
ts     = 0;
dt     = [0.1e-2; 0.0; 0.0];
twant  = [0.24; 0.25];
opti = zeros(4, 1, 'int64');
optr = ones(3, npde);

% index of next set of saved results
isav = 1;

% We run through the problem twice - once from t=0 up to t=twant(1), and
% then onwards to t=twant(2).
warning('Off');
fprintf('Example 1\n==========\n\n');
for iout = 1:2
  tout = twant(iout);
  [ts, dt, rwk, iwk, ind, ifail] = ...
  d03ra( ...
         ts, tout, dt, xmin, xmax, ymin, ymax, nx, ny, tols, tolt, ...
         @ex1_pdedef, @ex1_bndary, @ex1_pdeiv, @ex1_monitr, opti, ...
         optr, rwk, iwk, itrace, ind, 'lenlwk', lenlwk);

  fprintf('\nTime =%8.4f \n', ts);
  fprintf('Total number of accepted timesteps %5d \n', iwk(1));
  fprintf('Total number of rejected timesteps %5d \n', iwk(2));
  fprintf('\n              Total (max) number of\n');
  fprintf('          Residual  Jacobian    Newton   Lin sys\n');
  fprintf('             evals     evals     iters     iters\n');
  fprintf('At level\n');
  maxlev = 3;
  for j = 1:maxlev
    if iwk(j+2) ~= 0
      evals = iwk(j+2:maxlev:j+2+5*maxlev);
      fprintf('%8d%10d%10d%6d(%2d)%6d(%2d)\n', j, evals);
    end;
  end;
end


function [res] = ex1_pdedef(npts, npde, t, x, y, u, ut, ux, uy, uxx, uxy, uyy)
  % Evaluate the system of PDEs for this problem.

  global heat_release activ_energy diffusion damkohler;

  res = zeros(npts, npde);
  for i = 1:npts
    e = exp(-activ_energy/u(i,1));
    res(i,1) = ut(i,1) - diffusion*(uxx(i,1) + uyy(i,1)) - ...
               damkohler*(1.0 + heat_release - u(i,1))*e;
  end


function [u] = ex1_pdeiv(npts, npde, t, x, y)
  % Evaluate initial conditions for PDEs for this system.

  u = ones(npts, npde);


function [res] = ex1_bndary(npts, npde, t, x, y, u, ut, ux, uy, ...
                            nbpts, lbnd, res)
  % Implement boundary conditions for the domain.

  tol = 10*x02aj();
  for i = 1:nbpts
    j = lbnd(i);
    if (abs(x(j)) <= tol)
      res(j,1) = ux(j,1);
    elseif (abs(x(j)-1) <= tol)
      res(j,1) = u(j,1) - 1;
    elseif (abs(y(j)) <= tol)
      res(j,1) = uy(j,1);
    elseif (abs(y(j)-1) <= tol)
      res(j,1) = u(j,1) - 1;
    end
end


function [ierr] = ex1_monitr(npde, t, dt, dtnew, tlast, nlev, ...
                             ngpts, xpts, ypts, lsol, sol, ierr)
  % Save the results at specified times.
  times = [0.001; 0.228; 0.240; 0.25];

  global iout isav;
  global xsav ysav ssav1 tsav;

  % Save this time step only if the current time is equal to
  % or just past the next output time
  if (t < times(isav))
    return;
  end
  fprintf('Saving set %d at time = %f\n', isav, t);

  % Specify the grid level for extracting the solution.
  level = 1;
  npts = ngpts(level);
  ipsol = lsol(level);

  % Allocate space for saves the first time through.
  nside = round(sqrt(double(npts)));
  if isav == 1
    nsav = length(times);
    xsav = zeros(nside, 1);
    ysav = zeros(nside, 1);
    ssav1 = zeros(nside, nside, nsav);
    tsav = zeros(nsav, 1);
  end

  % Save this time step.
  tsav(isav) = t;
  for i = 1:nside
    for j = 1:nside
      k = (i-1)*nside + j;
      xsav(j) = xpts(k);
      ysav(i) = ypts(k);
      ssav1(j, i, isav) = sol(ipsol+k);
    end
  end

  % Look for the next time step, unless called for very last time.
  if ~(tlast && iout == 2)
    isav = isav+1;
  end


function ex1_plot
  % Plot the results.
  global xsav ysav ssav1 tsav isav;
  % Set the colours for each surface, and the array of handles for the plots.
  colours = [[1 0 0]; [0 1 0]; [1 0 1]; [0 0 1]; [0 1 1]; [2 1 0]; [2 0 1];
             [1 0 2]; [1 2 0]; [0 1 2]; [0 2 1]; [3 2 1]; [2 3 1]; [2 1 3];
             [3 1 2]; [1 3 2]; [1 2 3]; [0.3 2 1]; [2 0.3 1]; [2 1 0.3]; 
             [0.3 1 2]; [1 0.3 2]; [1 2 0.3]; [0.3 0.2 1]; [0.2 0.3 1];
             [0.2 1 0.3]; [0.3 1 0.2]; [1 0.3 0.2]; [1 0.2 0.3]];
             
  h = zeros(isav,1);
  % Check the allocation of colours against the number to be plotted.
  if (length(colours(:,1)) < isav)
    fprintf(['Not enough colours allocated in example1_plot: ', ...
             '%d needed\n'], isav);
  end
  % Plot all the surfaces, and colour each one.
  for i = 1:isav
    h(i) = mesh(xsav, ysav, ssav1(:,:,i));
    hold on;
    set(h(i), 'EdgeColor', colours(i,:));
  end
  hold off;
  % Label the axes, and set the title.
  xlabel('x');
  ylabel('y');
  zlabel('U(x,y,t)');
  title({'Model for a Single, One-step Reaction', ...
         'of a Mixture of Two Chemicals'});
  % Add a legend, using the plot handles and the array of time values.
  hleg = legend(h, num2str(tsav(1:isav)));
  % Set the view to something nice (determined empirically).
  view(15, 30);
  
  
function ex2
  % Second d03ra example.  Concentration distribution in a multispecies
  % (predator-prey) food web model.
  global alpha beta;
  global iout icount;

  % Set values for problem parameters.
  npde = int64(2);
  npts = int64(4000);
  xmin = 0;  xmax = 1;
  ymin = 0;  ymax = 1;
  tols = 0.075;
  tolt = 0.1;
  nx = int64(21);
  ny = int64(21);
  alpha = 50.0;
  beta = 300.0;

  maxlev = 4;
  leniwk = npts*(14 + 5*maxlev) + 7*maxlev + 2;
  lenlwk = npts + 1;
  lenrwk = npts*npde*(5*maxlev + 18*npde + 9) + 2*npts;

  % The above are the recommended values from the documentation, but d03ra
  % returns with ifail = 3 and the minumum workspace size requirements if
  % these are used.  So we beef them up here with some fudge factors.
  leniwk = 2*leniwk;
  lenlwk = int64(3*lenlwk);
  lenrwk = 3*lenrwk;

  % Create some work arrays - see NAG documentation.
  rwk = zeros(lenrwk, 1);
  iwk = zeros(leniwk, 1, 'int64');

% Initialize some problem variables.
  ind    = int64(0);
  itrace = int64(0);
  ts     = 0;
  dt     = [0.5e-3; 1.0e-6; 0.0];
  twant  = [0.01; 0.025];
  opti    = zeros(4, 1, 'int64');
  opti(1) = 4;
  optr      = ones(3, npde);
  optr(1,1) = 250;
  optr(1,2) = 1.5e+6;

  % Counts how many times monitor routine is called.
  icount = 0;
  fprintf('\n\nExample 2\n=========\n\n');

  % We run through the problem twice - once from t=0 up to t=twant(1), and
  % then onwards to t=twant(2).
  for iout = 1:2
    tout = twant(iout);
    [ts, dt, rwk, iwk, ind, ifail] = ...
    d03ra( ...
           ts, tout, dt, xmin, xmax, ymin, ymax, nx, ny, tols, tolt, ...
           @ex2_pdedef, @ex2_bndary, @ex2_pdeiv, @ex2_monitr, opti, optr, ...
           rwk, iwk, itrace, ind, 'lenlwk', lenlwk);

    % Output some statistics.
    fprintf('\nTime =%8.4f \n', ts);
    fprintf('Total number of accepted timesteps %5d \n', iwk(1));
    fprintf('Total number of rejected timesteps %5d \n', iwk(2));
    fprintf('\n              Total (max) number of\n');
    fprintf('          Residual  Jacobian    Newton   Lin sys\n');
    fprintf('             evals     evals     iters     iters\n');
    fprintf('At level\n');
    for j = 1:maxlev
      if iwk(j+2) ~= 0
        evals = iwk(j+2:maxlev:j+2+5*maxlev);
        fprintf('%8d%10d%10d%6d(%2d)%6d(%2d)\n', j, evals);
      end;
    end;
end


function [res] = ex2_pdedef(npts, npde, t, x, y, u, ut, ux, uy, ...
                            uxx, uxy, uyy)
  % Evaluate the system of PDEs for this problem.

  global alpha beta;

  fp = 4*pi;
  res = zeros(npts, npde);
  for i = 1:npts
    b1 = 1 + alpha*x(i)*y(i) + beta*sin(fp*x(i))*sin(fp*y(i));
    res(i,1) = ut(i,1) - uxx(i,1) - uyy(i,1) - ...
               u(i,1)*(b1 - u(i,1) - 0.5e-6*u(i,2));
    res(i,2) = -0.05*(uxx(i,2) + uyy(i,2)) - ...
               u(i,2)*(-b1 + 1.0e4*u(i,1) - u(i,2));
  end


function [u] = ex2_pdeiv(npts, npde, t, x, y)
  % Evaluate initial conditions for PDEs for this system.

  global alpha beta;

  fp = 4*pi;
  u = zeros(npts, npde);
  for i = 1:npts
    u(i,1) = 10 + (16*x(i)*(1 - x(i))*y(i)*(1 - y(i)))^2;
    u(i,2) = -1 - alpha*x(i)*y(i) - beta*sin(fp*x(i))*sin(fp*y(i)) + ...
             1.0e4*u(i,1);
  end


function [res] = ex2_bndary(npts, npde, t, x, y, u, ut, ux, uy, ...
                            nbpts, lbnd, res)
  % Implement boundary conditions for the domain.

  tol = 10*x02aj;
  for i = 1:nbpts
    j = lbnd(i);
    if (abs(x(j)) <= tol || abs(x(j)-1) <= tol)
      res(j,1) = ux(j,1);
      res(j,2) = ux(j,2);
    elseif (abs(y(j)) <= tol || abs(y(j)-1) <= tol)
      res(j,1) = uy(j,1);
      res(j,2) = uy(j,2);
    end
  end


function [ierr] = ex2_monitr(npde, t, dt, dtnew, tlast, nlev, ...
                             ngpts, xpts, ypts, lsol, sol, ierr)
  % Monitor the results at specified intervals.
           
  global iout icount;
  global xsav ysav ssav2 tsav isav;

  % Specify the maximum number of times the results get saved, and the
  % interval for saves.
%  nsav = 5;
%  nint = 10;
 size(xpts);
 size(ypts);

  % Do we want to save this time?
  icount = icount+1;
 % if mod(icount,nint) ~= 1 && ~tlast
 %   return;
 % end

 % isav = 1 + round(icount/nint);
 % if isav > nsav
 %   fprintf('Not enough save space allocated in monitr2: %d\n', ...
 %           nsav);
 %   return;
 % end
  isav = icount;

  % Specify the grid level for extracting solution.
  level = 1;
  npts = ngpts(level);
  ipsol = lsol(level);

  % Allocate space for saves the first time through.
  nside = round(sqrt(double(npts)));

  % Save this time step.  For this problem, calculating two
  % concentrations as a function of x and y (and time).
  tsav(isav) = t;
  if (icount==1)
    xsav(1:nside) = xpts(1:nside);
    ysav(1:nside) = ypts(1:nside:nside^2);
  end

  s1 = reshape(sol(ipsol+1:ipsol+nside^2),[nside,nside]);
  s2 = reshape(sol(ipsol+npts+1:ipsol+npts+nside^2),[nside,nside]);
  ssav2(:,:,isav,1) = s1';
  ssav2(:,:,isav,2) = s2';

function ex2_plot(id)
  % Plot the results.
  global xsav ysav ssav2 tsav isav;
  % Check the value of the array identifier.
  if (id ~= 1 && id ~= 2)
    fprintf('Illegal value for array identifier in plot: %d\n', id);
    return;
  end
  % Label the axes, and set the title according to the array identifier.
  v = ssav2(:,:,:,id);
  slice(xsav, ysav, tsav, v ,[0,0.4,0.6,0.8],[1],[0.01,0.02]);
  colormap hot;
  colorbar;
  xlabel('x');
  ylabel('y');
  zlabel('time');
  if (id == 1)
    title({'Multispecies Food Web Model:', ...
           'Time-dependent Predator Concentration'});
  else
    title({'Multispecies Food Web Model:', ...
           'Time-dependent Prey Concentration'});
  end
  % Add a legend, using the plot handles and the array of time values.
  % Set the view to something nice (determined empirically).
  view(24,42);

d03ra example results

Example 1
==========

Saving set 1 at time = 0.001000
Saving set 2 at time = 0.228996
Saving set 3 at time = 0.240000

Time =  0.2400 
Total number of accepted timesteps    75 
Total number of rejected timesteps     0 

              Total (max) number of
          Residual  Jacobian    Newton   Lin sys
             evals     evals     iters     iters
At level
       1       600        75   150(159)     2( 2)
Saving set 4 at time = 0.250000

Time =  0.2500 
Total number of accepted timesteps   180 
Total number of rejected timesteps     1 

              Total (max) number of
          Residual  Jacobian    Newton   Lin sys
             evals     evals     iters     iters
At level
       1      1468       181   382(391)     4( 2)
       2       662        82   170(170)     4( 1)
       3       177        22    45(45)     3( 1)


Example 2
=========


Time =  0.0100 
Total number of accepted timesteps    14 
Total number of rejected timesteps     0 

              Total (max) number of
          Residual  Jacobian    Newton   Lin sys
             evals     evals     iters     iters
At level
       1       196        14    28(39)     2( 2)
       2        84         6    12(19)     2( 3)

Time =  0.0250 
Total number of accepted timesteps    29 
Total number of rejected timesteps     0 

              Total (max) number of
          Residual  Jacobian    Newton   Lin sys
             evals     evals     iters     iters
At level
       1       406        29    58(84)     2( 2)
       2       294        21    42(64)     2( 3)
       3        98         7    14(28)     2( 3)