Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_pde_1d_parab_coll_interp (d03py)

## Purpose

nag_pde_1d_parab_coll_interp (d03py) may be used in conjunction with either nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj). It computes the solution and its first derivative at user-specified points in the spatial coordinate.

## Syntax

[up, rsave, ifail] = d03py(u, xbkpts, npoly, xp, itype, rsave, 'npde', npde, 'nbkpts', nbkpts, 'npts', npts, 'intpts', intpts)
[up, rsave, ifail] = nag_pde_1d_parab_coll_interp(u, xbkpts, npoly, xp, itype, rsave, 'npde', npde, 'nbkpts', nbkpts, 'npts', npts, 'intpts', intpts)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: lrsave has been removed from the interface
.

## Description

nag_pde_1d_parab_coll_interp (d03py) is an interpolation function for evaluating the solution of a system of partial differential equations (PDEs), or the PDE components of a system of PDEs with coupled ordinary differential equations (ODEs), at a set of user-specified points. The solution of a system of equations can be computed using nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj) on a set of mesh points; nag_pde_1d_parab_coll_interp (d03py) can then be employed to compute the solution at a set of points other than those originally used in nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj). It can also evaluate the first derivative of the solution. Polynomial interpolation is used between each of the break points xbkpts(i)${\mathbf{xbkpts}}\left(\mathit{i}\right)$, for i = 1,2,,nbkpts$\mathit{i}=1,2,\dots ,{\mathbf{nbkpts}}$. When the derivative is needed (itype = 2${\mathbf{itype}}=2$), the array xp(intpts)${\mathbf{xp}}\left({\mathbf{intpts}}\right)$ must not contain any of the break points, as the method, and consequently the interpolation scheme, assumes that only the solution is continuous at these points.

None.

## Parameters

Note: the parameters u, npts, npde, xbkpts, nbkpts, rsave and lrsave must be supplied unchanged from either nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj).

### Compulsory Input Parameters

1:     u(npde,npts) – double array
npde, the first dimension of the array, must satisfy the constraint npde1${\mathbf{npde}}\ge 1$.
The PDE part of the original solution returned in the parameter u by the function nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj).
2:     xbkpts(nbkpts) – double array
nbkpts, the dimension of the array, must satisfy the constraint nbkpts2${\mathbf{nbkpts}}\ge 2$.
xbkpts(i)${\mathbf{xbkpts}}\left(\mathit{i}\right)$, for i = 1,2,,nbkpts$\mathit{i}=1,2,\dots ,{\mathbf{nbkpts}}$, must contain the break points as used by nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj).
Constraint: xbkpts(1) < xbkpts(2) < < ${\mathbf{xbkpts}}\left(1\right)<{\mathbf{xbkpts}}\left(2\right)<\cdots <{\mathbf{xbkpts}}\left({\mathbf{nbkpts}}\right)$.
3:     npoly – int64int32nag_int scalar
The degree of the Chebyshev polynomial used for approximation as used by nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj).
Constraint: 1npoly49$1\le {\mathbf{npoly}}\le 49$.
4:     xp(intpts) – double array
intpts, the dimension of the array, must satisfy the constraint intpts1${\mathbf{intpts}}\ge 1$.
xp(i)${\mathbf{xp}}\left(\mathit{i}\right)$, for i = 1,2,,intpts$\mathit{i}=1,2,\dots ,{\mathbf{intpts}}$, must contain the spatial interpolation points.
Constraints:
• xbkpts(1)xp(1) < xp(2) < < xp(intpts)${\mathbf{xbkpts}}\left(1\right)\le {\mathbf{xp}}\left(1\right)<{\mathbf{xp}}\left(2\right)<\cdots <{\mathbf{xp}}\left({\mathbf{intpts}}\right)\le {\mathbf{xbkpts}}\left({\mathbf{nbkpts}}\right)$;
• if itype = 2${\mathbf{itype}}=2$, xp(i)xbkpts(j)${\mathbf{xp}}\left(\mathit{i}\right)\ne {\mathbf{xbkpts}}\left(\mathit{j}\right)$, for i = 1,2,,intpts$\mathit{i}=1,2,\dots ,{\mathbf{intpts}}$ and j = 2,3,,nbkpts1$\mathit{j}=2,3,\dots ,{\mathbf{nbkpts}}-1$.
5:     itype – int64int32nag_int scalar
Specifies the interpolation to be performed.
itype = 1${\mathbf{itype}}=1$
The solution at the interpolation points are computed.
itype = 2${\mathbf{itype}}=2$
Both the solution and the first derivative at the interpolation points are computed.
Constraint: itype = 1${\mathbf{itype}}=1$ or 2$2$.
6:     rsave(lrsave) – double array
The array rsave contains information required by nag_pde_1d_parab_coll_interp (d03py) as returned by nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj). The contents of rsave must not be changed from the call to nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj). Some elements of this array are overwritten on exit.

### Optional Input Parameters

1:     npde – int64int32nag_int scalar
Default: The first dimension of the array u.
The number of PDEs.
Constraint: npde1${\mathbf{npde}}\ge 1$.
2:     nbkpts – int64int32nag_int scalar
Default: The dimension of the array xbkpts.
The number of break points.
Constraint: nbkpts2${\mathbf{nbkpts}}\ge 2$.
3:     npts – int64int32nag_int scalar
Default: The second dimension of the array u.
The number of mesh points as used by nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj).
Constraint: npts = (nbkpts1) × npoly + 1${\mathbf{npts}}=\left({\mathbf{nbkpts}}-1\right)×{\mathbf{npoly}}+1$.
4:     intpts – int64int32nag_int scalar
Default: The dimension of the array xp.
The number of interpolation points.
Constraint: intpts1${\mathbf{intpts}}\ge 1$.

lrsave

### Output Parameters

1:     up(npde,intpts,itype) – double array
If itype = 1${\mathbf{itype}}=1$, up(i,j,1)${\mathbf{up}}\left(\mathit{i},\mathit{j},1\right)$, contains the value of the solution Ui(xj,tout)${U}_{\mathit{i}}\left({x}_{\mathit{j}},{t}_{\mathrm{out}}\right)$, at the interpolation points xj = xp(j)${x}_{\mathit{j}}={\mathbf{xp}}\left(\mathit{j}\right)$, for j = 1,2,,intpts$\mathit{j}=1,2,\dots ,{\mathbf{intpts}}$ and i = 1,2,,npde$\mathit{i}=1,2,\dots ,{\mathbf{npde}}$.
If itype = 2${\mathbf{itype}}=2$, up(i,j,1)${\mathbf{up}}\left(\mathit{i},\mathit{j},1\right)$ contains Ui(xj,tout)${U}_{\mathit{i}}\left({x}_{\mathit{j}},{t}_{\mathrm{out}}\right)$ and up(i,j,2)${\mathbf{up}}\left(\mathit{i},\mathit{j},2\right)$ contains (Ui)/(x) $\frac{\partial {U}_{\mathit{i}}}{\partial x}$ at these points.
2:     rsave(lrsave) – double array
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, itype ≠ 1${\mathbf{itype}}\ne 1$ or 2$2$, or npoly < 1${\mathbf{npoly}}<1$, or npde < 1${\mathbf{npde}}<1$, or nbkpts < 2${\mathbf{nbkpts}}<2$, or intpts < 1${\mathbf{intpts}}<1$, or npts ≠ (nbkpts − 1) × npoly + 1${\mathbf{npts}}\ne \left({\mathbf{nbkpts}}-1\right)×{\mathbf{npoly}}+1$, or xbkpts(i)${\mathbf{xbkpts}}\left(\mathit{i}\right)$, for i = 1,2, … ,nbkpts$\mathit{i}=1,2,\dots ,{\mathbf{nbkpts}}$, are not ordered.
ifail = 2${\mathbf{ifail}}=2$
On entry, the interpolation points xp(i)${\mathbf{xp}}\left(\mathit{i}\right)$, for i = 1,2,,intpts$\mathit{i}=1,2,\dots ,{\mathbf{intpts}}$, are not in strictly increasing order, or when itype = 2${\mathbf{itype}}=2$, at least one of the interpolation points stored in xp is equal to one of the break points stored in xbkpts.
ifail = 3${\mathbf{ifail}}=3$
You are attempting extrapolation, that is, one of the interpolation points xp(i)${\mathbf{xp}}\left(i\right)$, for some i$i$, lies outside the interval [xbkpts(1),${\mathbf{xbkpts}}\left(1\right),{\mathbf{xbkpts}}\left({\mathbf{nbkpts}}\right)$]. Extrapolation is not permitted.

## Accuracy

See the documents for nag_pde_1d_parab_coll (d03pd) or nag_pde_1d_parab_dae_coll (d03pj).

None.

## Example

```function nag_pde_1d_parab_coll_interp_example
m = int64(0);
ts = 0;
tout = 0.0001;
u = zeros(2, 28);
xbkpts = [-1;
-0.7777777777777778;
-0.5555555555555556;
-0.3333333333333333;
-0.1111111111111111;
0.1111111111111111;
0.3333333333333333;
0.5555555555555556;
0.7777777777777778;
1];
npoly = int64(3);
acc = 0.0001;
rsave = zeros(2085, 1);
isave = zeros(80, 1, 'int64');
itrace = int64(0);
ind = int64(0);
cwsav = {''; ''; ''; ''; ''; ''; ''; ''; ''; ''};
lwsav = false(100, 1);
iwsav = zeros(505, 1, 'int64');
rwsav = zeros(1100, 1);
xp = [-1;
-0.6;
-0.2;
0.2;
0.6;
1];
itype = int64(1);
[ts, u, x, rsave, isave, ind, user, cwsav, lwsav, iwsav, rwsav, ifail] = ...
nag_pde_1d_parab_coll(m, ts, tout, @pdedef, @bndary, u, xbkpts, npoly, ...
@uinit, acc, rsave, isave, itask, itrace, ind, ...
cwsav, lwsav, iwsav, rwsav);
[up, rsaveOut, ifail] = nag_pde_1d_parab_coll_interp(u, xbkpts, npoly, xp, itype, rsave);
ifail

function [p, q, r, ires, user] = pdedef(npde, t, x, nptl, u, ux, ires, user)
p = zeros(npde, npde, nptl);
q = zeros(npde, nptl);
r = zeros(npde, nptl);

for i = 1:double(nptl)
q(1,i) = u(2,i);
q(2,i) = u(1,i)*ux(2,i) - ux(1,i)*u(2,i);
r(1,i) = ux(1,i);
r(2,i) = ux(2,i);
p(1,1,i) = 0;
p(1,2,i) = 0;
p(2,1,i) = 0;
p(2,2,i) = 1;
end;

function [beta, gamma, ires, user] = bndary(npde, t, u, ux, ibnd, ires, user)
beta = zeros(npde, 1);
gamma = zeros(npde, 1);

if (ibnd == 0)
beta(1) = 1;
gamma(1) = 0;
beta(2) = 0;
gamma(2) = u(1) - 1;
else
beta(1) = 1;
gamma(1) = 0;
beta(2) = 0;
gamma(2) = u(1) + 1;
end

function [u, user] = uinit(npde, npts, x, user)
u = zeros(npde, npts);

piby2 = pi/2;
for i = 1:double(npts)
u(1,i) = -sin(piby2*x(i));
u(2,i) = -piby2*piby2*u(1,i);
end
```
```

ifail =

0

```
```function d03py_example
m = int64(0);
ts = 0;
tout = 0.0001;
u = zeros(2, 28);
xbkpts = [-1;
-0.7777777777777778;
-0.5555555555555556;
-0.3333333333333333;
-0.1111111111111111;
0.1111111111111111;
0.3333333333333333;
0.5555555555555556;
0.7777777777777778;
1];
npoly = int64(3);
acc = 0.0001;
rsave = zeros(2085, 1);
isave = zeros(80, 1, 'int64');
itrace = int64(0);
ind = int64(0);
cwsav = {''; ''; ''; ''; ''; ''; ''; ''; ''; ''};
lwsav = false(100, 1);
iwsav = zeros(505, 1, 'int64');
rwsav = zeros(1100, 1);
xp = [-1;
-0.6;
-0.2;
0.2;
0.6;
1];
itype = int64(1);
[ts, u, x, rsave, isave, ind, user, cwsav, lwsav, iwsav, rwsav, ifail] = ...
d03pd(m, ts, tout, @pdedef, @bndary, u, xbkpts, npoly, ...
@uinit, acc, rsave, isave, itask, itrace, ind, ...
cwsav, lwsav, iwsav, rwsav);
[up, rsaveOut, ifail] = d03py(u, xbkpts, npoly, xp, itype, rsave);
ifail

function [p, q, r, ires, user] = pdedef(npde, t, x, nptl, u, ux, ires, user)
p = zeros(npde, npde, nptl);
q = zeros(npde, nptl);
r = zeros(npde, nptl);

for i = 1:double(nptl)
q(1,i) = u(2,i);
q(2,i) = u(1,i)*ux(2,i) - ux(1,i)*u(2,i);
r(1,i) = ux(1,i);
r(2,i) = ux(2,i);
p(1,1,i) = 0;
p(1,2,i) = 0;
p(2,1,i) = 0;
p(2,2,i) = 1;
end;

function [beta, gamma, ires, user] = bndary(npde, t, u, ux, ibnd, ires, user)
beta = zeros(npde, 1);
gamma = zeros(npde, 1);

if (ibnd == 0)
beta(1) = 1;
gamma(1) = 0;
beta(2) = 0;
gamma(2) = u(1) - 1;
else
beta(1) = 1;
gamma(1) = 0;
beta(2) = 0;
gamma(2) = u(1) + 1;
end

function [u, user] = uinit(npde, npts, x, user)
u = zeros(npde, npts);

piby2 = pi/2;
for i = 1:double(npts)
u(1,i) = -sin(piby2*x(i));
u(2,i) = -piby2*piby2*u(1,i);
end
```
```

ifail =

0

```