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_extractgrid (d03rz)

Purpose

nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) is designed to be used in conjunction with nag_pde_2d_gen_order2_rectilinear (d03rb). It can be called from the monitr to obtain the number of grid points and their (x,y)(x,y) coordinates on a solution grid.

Syntax

[npts, x, y, ifail] = d03rz(level, nlev, xmin, ymin, dxb, dyb, lgrid, istruc, lenxy)
[npts, x, y, ifail] = nag_pde_2d_gen_order2_rectilinear_extractgrid(level, nlev, xmin, ymin, dxb, dyb, lgrid, istruc, lenxy)

Description

nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) extracts the number of grid points and their (x,y)(x,y) coordinates on a specific solution grid produced by nag_pde_2d_gen_order2_rectilinear (d03rb). It must be called only from within the monitr. The parameters nlev, xmin, ymin, dxb, dyb, lgrid and istruc to monitr must be passed unchanged to nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz).

References

None.

Parameters

Compulsory Input Parameters

1:     level – int64int32nag_int scalar
The grid level at which the coordinates are required.
Constraint: 1levelnlev1levelnlev.
2:     nlev – int64int32nag_int scalar
3:     xmin – double scalar
4:     ymin – double scalar
5:     dxb – double scalar
6:     dyb – double scalar
7:     lgrid( : :) – int64int32nag_int array
8:     istruc( : :) – int64int32nag_int array
Note: the dimension of the array lgrid must be at least nlevnlev.
The dimension of the array istruc must be at least lgrid(nlev) + 2 × nrows + npts + 1lgridnlev+2×nrows+npts+1 where nrowsnrows is stored in istruc(lgrid(nlev))istruc(lgrid(nlev)) and is the number of rows in the grid at level nlev.
nlev, xmin, ymin, dxb, dyb, lgrid and istruc as supplied to monitr must be passed unchanged to nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz).
9:     lenxy – int64int32nag_int scalar
The dimension of the arrays x and y as declared in the (sub)program from which nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) is called.
Constraint: lenxynptslenxynpts.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     npts – int64int32nag_int scalar
The number of grid points in the grid level level.
2:     x(lenxy) – double array
3:     y(lenxy) – double array
x(i)xi and y(i)yi contain the (x,y)(x,y) coordinates respectively of the iith grid point, for i = 1,2,,nptsi=1,2,,npts.
4:     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:
  ifail = 1ifail=1
On entry,level < 1level<1,
orlevel > nlevlevel>nlev.
  ifail = 2ifail=2
The dimension of the arrays x and y is too small for the requested grid level, i.e., lenxy < nptslenxy<npts.

Accuracy

Not applicable.

Further Comments

None.

Example

function nag_pde_2d_gen_order2_rectilinear_extractgrid_example
ts = 0;
twant = [0.25; 1];
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);
global iout;
for iout = 1:2
  tout = twant(iout);
  [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);
  fprintf('\nStatistics\n');
  fprintf('Time = %8.4f\n', ts);
  fprintf('Total number of accepted timesteps = %d\n', iwk(1));
  fprintf('Total number of rejected timesteps = %d\n', iwk(2));
  fprintf('\n             T o t a l  n u m b e r  o f    \n');
  fprintf('          Residual  Jacobian    Newton  Lin sys\n');
  fprintf('             evals     evals     iters    iters\n');
  fprintf(' At level \n');
  maxlev = double(opti(1));
  for j = 1:maxlev
    if (iwk(j+2) ~= 0)
      fprintf('%6d %10d %10d %10d %10d\n', j, iwk(j+2), iwk(j+2+maxlev), ...
                                 iwk(j+2+2*maxlev), iwk(j+2+3*maxlev) );
    end
  end
  fprintf('\n             M a x i m u m  n u m b e r  o f\n');
  fprintf('              Newton iters    Lin sys iters \n');
  fprintf(' At level \n');
  for j = 1:maxlev
    if (iwk(j+2) ~= 0)
      fprintf('%6d%14d%14d\n', j, iwk(j+2+4*maxlev), iwk(j+2+5*maxlev) );
    end
  end

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 [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 [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)

  lenxy=int64(2500);

  global iout;

  if (tlast && iout == 2)
    for level = 1:double(nlev)
      ipsol = lsol(level);

      % Get grid information
      [npts, x, y, ifail] = ...
    nag_pde_2d_gen_order2_rectilinear_extractgrid(int64(level), nlev,  ...
                                 xmin, ymin, dxb, dyb, lgrid, istruc, lenxy);
      if (ifail ~= 0)
        ierr = 1;
      elseif (level == 1)
        % Get exact solution
        [uex] = pdeiv(npts, npde, t, x, y);

        fprintf('\n Solution at every 2nd grid point in level 1 at time %8.4f\n', t);
        fprintf('       x          y        approx u     exact u    approx v     exact v\n');
        for ipt = 1:2:double(npts)
          fprintf('%11.4E %11.4e %11.4e %11.4e %11.4e %11.4e\n', x(ipt), ...
           y(ipt), sol(double(ipsol)+double(ipt)), uex(ipt,1), ...
           sol(double(ipsol)+double(npts)+double(ipt)), uex(ipt,2));
        end
        fprintf('\n');
      end
    end
  end
 

Statistics
Time =   0.2500
Total number of accepted timesteps = 14
Total number of rejected timesteps = 0

             T o t a l  n u m b e r  o f    
          Residual  Jacobian    Newton  Lin sys
             evals     evals     iters    iters
 At level 
     1        196         14         28         14
     2        196         14         28         22
     3        196         14         28         25
     4        196         14         28         31
     5        141         10         21         29

             M a x i m u m  n u m b e r  o f
              Newton iters    Lin sys iters 
 At level 
     1             2             1
     2             2             1
     3             2             1
     4             2             2
     5             3             2

 Solution at every 2nd grid point in level 1 at time   1.0000
       x          y        approx u     exact u    approx v     exact v
 0.0000E+00  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 2.0000E-01  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  1.0000e-01  5.0022e-01  5.0000e-01  9.9978e-01  1.0000e+00
 3.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 5.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 7.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 0.0000E+00  2.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 2.0000E-01  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 4.0000E-01  2.0000e-01  5.0011e-01  5.0000e-01  9.9989e-01  1.0000e+00
 6.0000E-01  2.0000e-01  4.9994e-01  5.0000e-01  1.0001e+00  1.0000e+00
 8.0000E-01  2.0000e-01  4.9999e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E+00  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  3.0000e-01  4.9997e-01  5.0048e-01  1.0000e+00  9.9952e-01
 3.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 5.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 7.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  3.0000e-01  5.0003e-01  5.0000e-01  9.9997e-01  1.0000e+00
 0.0000E+00  4.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  4.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 4.0000E-01  4.0000e-01  5.0018e-01  5.0000e-01  9.9982e-01  1.0000e+00
 8.0000E-01  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E+00  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  5.0000e-01  7.4998e-01  7.5000e-01  7.5002e-01  7.5000e-01
 3.0000E-01  5.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 5.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 7.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  5.0000e-01  5.0005e-01  5.0000e-01  9.9995e-01  1.0000e+00
 0.0000E+00  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 4.0000E-01  6.0000e-01  4.9998e-01  5.0048e-01  1.0000e+00  9.9952e-01
 6.0000E-01  6.0000e-01  4.9991e-01  5.0000e-01  1.0001e+00  1.0000e+00
 8.0000E-01  6.0000e-01  4.9983e-01  5.0000e-01  1.0002e+00  1.0000e+00
 1.0000E+00  6.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  7.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
 3.0000E-01  7.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 5.0000E-01  7.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 7.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 0.0000E+00  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 4.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 6.0000E-01  8.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 8.0000E-01  8.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  9.0000e-01  7.5001e-01  7.5000e-01  7.4999e-01  7.5000e-01
 3.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
 5.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
 7.0000E-01  9.0000e-01  4.9994e-01  5.0048e-01  1.0001e+00  9.9952e-01
 0.0000E+00  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 4.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 6.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 8.0000E-01  1.0000e+00  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01


Statistics
Time =   1.0000
Total number of accepted timesteps = 45
Total number of rejected timesteps = 0

             T o t a l  n u m b e r  o f    
          Residual  Jacobian    Newton  Lin sys
             evals     evals     iters    iters
 At level 
     1        630         45         90         45
     2        630         45         90         78
     3        630         45         90         87
     4        630         45         90        124
     5        575         41         83        122

             M a x i m u m  n u m b e r  o f
              Newton iters    Lin sys iters 
 At level 
     1             2             1
     2             2             1
     3             2             1
     4             2             2
     5             3             2

function d03rz_example
ts = 0;
twant = [0.25; 1];
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);
global iout;
for iout = 1:2
  tout = twant(iout);
  [ts, dt, rwk, iwk, ind, ifail] = ...
    d03rb(ts, tout, dt, tols, tolt, @inidom, @pdedef, ...
    @bndary, @pdeiv, @monitr, opti, optr, ...
    rwk, iwk, lenlwk, itrace, ind);
  fprintf('\nStatistics\n');
  fprintf('Time = %8.4f\n', ts);
  fprintf('Total number of accepted timesteps = %d\n', iwk(1));
  fprintf('Total number of rejected timesteps = %d\n', iwk(2));
  fprintf('\n             T o t a l  n u m b e r  o f    \n');
  fprintf('          Residual  Jacobian    Newton  Lin sys\n');
  fprintf('             evals     evals     iters    iters\n');
  fprintf(' At level \n');
  maxlev = double(opti(1));
  for j = 1:maxlev
    if (iwk(j+2) ~= 0)
      fprintf('%6d %10d %10d %10d %10d\n', j, iwk(j+2), iwk(j+2+maxlev), ...
                                 iwk(j+2+2*maxlev), iwk(j+2+3*maxlev) );
    end
  end
  fprintf('\n             M a x i m u m  n u m b e r  o f\n');
  fprintf('              Newton iters    Lin sys iters \n');
  fprintf(' At level \n');
  for j = 1:maxlev
    if (iwk(j+2) ~= 0)
      fprintf('%6d%14d%14d\n', j, iwk(j+2+4*maxlev), iwk(j+2+5*maxlev) );
    end
  end

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 [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 [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)

  lenxy=int64(2500);

  global iout;

  if (tlast && iout == 2)
    for level = 1:double(nlev)
      ipsol = lsol(level);

      % Get grid information
      [npts, x, y, ifail] = ...
    d03rz(int64(level), nlev, xmin, ymin, dxb, dyb, lgrid, istruc, lenxy);
      if (ifail ~= 0)
        ierr = 1;
      elseif (level == 1)
        % Get exact solution
        [uex] = pdeiv(npts, npde, t, x, y);

        fprintf('\n Solution at every 2nd grid point in level 1 at time %8.4f\n', t);
        fprintf('       x          y        approx u     exact u    approx v     exact v\n');
        for ipt = 1:2:double(npts)
          fprintf('%11.4E %11.4e %11.4e %11.4e %11.4e %11.4e\n', x(ipt), ...
           y(ipt), sol(double(ipsol)+double(ipt)), uex(ipt,1), ...
           sol(double(ipsol)+double(npts)+double(ipt)), uex(ipt,2));
        end
        fprintf('\n');
      end
    end
  end
 

Statistics
Time =   0.2500
Total number of accepted timesteps = 14
Total number of rejected timesteps = 0

             T o t a l  n u m b e r  o f    
          Residual  Jacobian    Newton  Lin sys
             evals     evals     iters    iters
 At level 
     1        196         14         28         14
     2        196         14         28         22
     3        196         14         28         25
     4        196         14         28         31
     5        141         10         21         29

             M a x i m u m  n u m b e r  o f
              Newton iters    Lin sys iters 
 At level 
     1             2             1
     2             2             1
     3             2             1
     4             2             2
     5             3             2

 Solution at every 2nd grid point in level 1 at time   1.0000
       x          y        approx u     exact u    approx v     exact v
 0.0000E+00  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 2.0000E-01  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  1.0000e-01  5.0022e-01  5.0000e-01  9.9978e-01  1.0000e+00
 3.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 5.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 7.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 0.0000E+00  2.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 2.0000E-01  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 4.0000E-01  2.0000e-01  5.0011e-01  5.0000e-01  9.9989e-01  1.0000e+00
 6.0000E-01  2.0000e-01  4.9994e-01  5.0000e-01  1.0001e+00  1.0000e+00
 8.0000E-01  2.0000e-01  4.9999e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E+00  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  3.0000e-01  4.9997e-01  5.0048e-01  1.0000e+00  9.9952e-01
 3.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 5.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 7.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  3.0000e-01  5.0003e-01  5.0000e-01  9.9997e-01  1.0000e+00
 0.0000E+00  4.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  4.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 4.0000E-01  4.0000e-01  5.0018e-01  5.0000e-01  9.9982e-01  1.0000e+00
 8.0000E-01  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E+00  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  5.0000e-01  7.4998e-01  7.5000e-01  7.5002e-01  7.5000e-01
 3.0000E-01  5.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 5.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 7.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  5.0000e-01  5.0005e-01  5.0000e-01  9.9995e-01  1.0000e+00
 0.0000E+00  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 4.0000E-01  6.0000e-01  4.9998e-01  5.0048e-01  1.0000e+00  9.9952e-01
 6.0000E-01  6.0000e-01  4.9991e-01  5.0000e-01  1.0001e+00  1.0000e+00
 8.0000E-01  6.0000e-01  4.9983e-01  5.0000e-01  1.0002e+00  1.0000e+00
 1.0000E+00  6.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  7.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
 3.0000E-01  7.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 5.0000E-01  7.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 7.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 9.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 0.0000E+00  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 4.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 6.0000E-01  8.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
 8.0000E-01  8.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
 1.0000E-01  9.0000e-01  7.5001e-01  7.5000e-01  7.4999e-01  7.5000e-01
 3.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
 5.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
 7.0000E-01  9.0000e-01  4.9994e-01  5.0048e-01  1.0001e+00  9.9952e-01
 0.0000E+00  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 2.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 4.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 6.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
 8.0000E-01  1.0000e+00  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01


Statistics
Time =   1.0000
Total number of accepted timesteps = 45
Total number of rejected timesteps = 0

             T o t a l  n u m b e r  o f    
          Residual  Jacobian    Newton  Lin sys
             evals     evals     iters    iters
 At level 
     1        630         45         90         45
     2        630         45         90         78
     3        630         45         90         87
     4        630         45         90        124
     5        575         41         83        122

             M a x i m u m  n u m b e r  o f
              Newton iters    Lin sys iters 
 At level 
     1             2             1
     2             2             1
     3             2             1
     4             2             2
     5             3             2


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