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_ode_bvp_ps_lin_solve (d02ue)

## Purpose

nag_ode_bvp_ps_lin_solve (d02ue) finds the solution of a linear constant coefficient boundary value problem by using the Chebyshev integration formulation on a Chebyshev Gauss–Lobatto grid.

## Syntax

[bmat, f, uc, resid, ifail] = d02ue(n, a, b, c, bmat, y, bvec, f, 'm', m)
[bmat, f, uc, resid, ifail] = nag_ode_bvp_ps_lin_solve(n, a, b, c, bmat, y, bvec, f, 'm', m)

## Description

nag_ode_bvp_ps_lin_solve (d02ue) solves the constant linear coefficient ordinary differential problem
 m ∑ fj + 1(dju)/(dxj) = f(x), x ∈ [a,b] j = 0
$∑ j=0 m fj+1 dju dxj = f(x) , x ∈ [a,b]$
subject to a set of m$m$ linear constraints at points yi[a,b] ${y}_{\mathit{i}}\in \left[a,b\right]$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$:
 m ∑ Bi,j + 1 ((dju)/(dxj))(x = yi) = βi, j = 0
$∑ j=0 m B i,j+1 ( dju dxj ) (x=yi) = βi ,$
where 1m4$1\le m\le 4$, B$B$ is an m × (m + 1)$m×\left(m+1\right)$ matrix of constant coefficients and βi${\beta }_{i}$ are constants. The points yi${y}_{i}$ are usually either a$a$ or b$b$.
The function f(x)$f\left(x\right)$ is supplied as an array of Chebyshev coefficients cj${c}_{j}$, j = 0,1,,n$j=0,1,\dots ,n$ for the function discretized on n + 1$n+1$ Chebyshev Gauss–Lobatto points (as returned by nag_ode_bvp_ps_lin_cgl_grid (d02uc)); the coefficients are normally obtained by a previous call to nag_ode_bvp_ps_lin_coeffs (d02ua). The solution and its derivatives (up to order m$m$) are returned, in the form of their Chebyshev series representation, as arrays of Chebyshev coefficients; subsequent calls to nag_ode_bvp_ps_lin_cgl_vals (d02ub) will return the corresponding function and derivative values at the Chebyshev Gauss–Lobatto discretization points on [a,b]$\left[a,b\right]$. Function and derivative values can be obtained on any uniform grid over the same range [a, b ]$\left[a,b\right]$ by calling the interpolation function nag_ode_bvp_ps_lin_grid_vals (d02uw).

## References

Clenshaw C W (1957) The numerical solution of linear differential equations in Chebyshev series Proc. Camb. Phil. Soc. 53 134–149
Coutsias E A, Hagstrom T and Torres D (1996) An efficient spectral method for ordinary differential equations with rational function coefficients Mathematics of Computation 65(214) 611–635
Greengard L (1991) Spectral integration and two-point boundary value problems SIAM J. Numer. Anal. 28(4) 1071–80
Lundbladh A, Hennigson D S and Johannson A V (1992) An efficient spectral integration method for the solution of the Navier–Stokes equations Technical report FFA–TN 1992–28 Aeronautical Research Institute of Sweden
Muite B K (2010) A numerical comparison of Chebyshev methods for solving fourth-order semilinear initial boundary value problems Journal of Computational and Applied Mathematics 234(2) 317–342

## Parameters

### Compulsory Input Parameters

1:     n – int64int32nag_int scalar
n$n$, where the number of grid points is n + 1$n+1$.
Constraint: n8${\mathbf{n}}\ge 8$ and n is even.
2:     a – double scalar
a$a$, the lower bound of domain [a,b]$\left[a,b\right]$.
Constraint: a < b${\mathbf{a}}<{\mathbf{b}}$.
3:     b – double scalar
b$b$, the upper bound of domain [a,b]$\left[a,b\right]$.
Constraint: b > a${\mathbf{b}}>{\mathbf{a}}$.
4:     c(n + 1${\mathbf{n}}+1$) – double array
The Chebyshev coefficients cj${c}_{j}$, j = 0,1,,n$j=0,1,\dots ,n$, for the right hand side of the boundary value problem. Usually these are obtained by a previous call of nag_ode_bvp_ps_lin_coeffs (d02ua).
5:     bmat(m,m + 1${\mathbf{m}}+1$) – double array
m, the first dimension of the array, must satisfy the constraint 1m4$1\le {\mathbf{m}}\le 4$.
bmat(i,j + 1)${\mathbf{bmat}}\left(\mathit{i},\mathit{j}+1\right)$ must contain the coefficients Bi,j + 1${B}_{\mathit{i},\mathit{j}+1}$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$ and j = 0,1,,m$\mathit{j}=0,1,\dots ,m$, in the problem formulation of Section [Description].
6:     y(m) – double array
m, the dimension of the array, must satisfy the constraint 1m4$1\le {\mathbf{m}}\le 4$.
The points, yi${y}_{\mathit{i}}$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$, where the boundary conditions are discretized.
7:     bvec(m) – double array
m, the dimension of the array, must satisfy the constraint 1m4$1\le {\mathbf{m}}\le 4$.
The values, βi${\beta }_{\mathit{i}}$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$, in the formulation of the boundary conditions given in Section [Description].
8:     f(m + 1${\mathbf{m}}+1$) – double array
The coefficients, fj${f}_{\mathit{j}}$, for j = 1,2,,m + 1$\mathit{j}=1,2,\dots ,m+1$, in the formulation of the linear boundary value problem given in Section [Description]. The highest order term, f(m + 1)${\mathbf{f}}\left({\mathbf{m}}+1\right)$, needs to be nonzero to have a well posed problem.

### Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the array bmat and the dimension of the arrays y, bvec. (An error is raised if these dimensions are not equal.)
The order, m$m$, of the boundary value problem to be solved.
Constraint: 1m4$1\le {\mathbf{m}}\le 4$.

None.

### Output Parameters

1:     bmat(m,m + 1${\mathbf{m}}+1$) – double array
The coefficients have been scaled to form an equivalent problem defined on the domain [1,1]$\left[-1,1\right]$.
2:     f(m + 1${\mathbf{m}}+1$) – double array
The coefficients have been scaled to form an equivalent problem defined on the domain [1,1]$\left[-1,1\right]$.
3:     uc(n + 1${\mathbf{n}}+1$,m + 1${\mathbf{m}}+1$) – double array
The Chebyshev coefficients in the Chebyshev series representations of the solution and derivatives of the solution to the boundary value problem. The n + 1$n+1$ elements uc(1 : n + 1,1)${\mathbf{uc}}\left(1:{\mathbf{n}}+1,1\right)$ contain the coefficients representing the solution U(xi)$U\left({x}_{\mathit{i}}\right)$, for i = 0,1,,n$\mathit{i}=0,1,\dots ,n$. uc(1 : n + 1,j + 1)${\mathbf{uc}}\left(1:{\mathbf{n}}+1,\mathit{j}+1\right)$ contains the coefficients representing the j$\mathit{j}$th derivative of U$U$, for j = 1,2,,m$\mathit{j}=1,2,\dots ,m$.
4:     resid – double scalar
The maximum residual resulting from substituting the solution vectors returned in uc into both linear equations of Section [Description] representing the linear boundary value problem and associated boundary conditions. That is
max { maxi = 1,m
 (m ) ∑ Bi,j + 1 ((dju)/(dxj))(x = y_i) − βi j = 0
, maxi = 1, n + 1
 (m ) ∑ fj + 1 ((dju)/(dxj))(x = x_i) − f(x) j = 0
} .
$max { max i=1,m ( | ∑ j=0 m B i,j+1 ( dju dxj ) ( x=yi ) - βi | ) , max i=1, n+1 ( | ∑ j=0 m f j+1 ( dju dxj ) (x=xi) - f(x) | ) } .$
5:     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:

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

ifail = 1${\mathbf{ifail}}=1$
Constraint: n is even.
Constraint: n8${\mathbf{n}}\ge 8$.
ifail = 2${\mathbf{ifail}}=2$
Constraint: a < b${\mathbf{a}}<{\mathbf{b}}$.
ifail = 3${\mathbf{ifail}}=3$
On entry, f(m + 1) = 0.0${\mathbf{f}}\left({\mathbf{m}}+1\right)=0.0$.
ifail = 6${\mathbf{ifail}}=6$
Constraint: 1m4$1\le {\mathbf{m}}\le 4$. Constraint: 1m4$1\le {\mathbf{m}}\le 4$.
ifail = 7${\mathbf{ifail}}=7$
Internal error while unpacking matrix during iterative refinement.
ifail = 8${\mathbf{ifail}}=8$
Singular matrix encountered during iterative refinement.
W ifail = 9${\mathbf{ifail}}=9$
During iterative refinement, the maximum number of iterations was reached.
W ifail = 10${\mathbf{ifail}}=10$
During iterative refinement, convergence was achieved, but the residual is more than 100 × machine precision.
ifail = 999${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

The accuracy should be close to machine precision for well conditioned boundary value problems.

The number of operations is of the order n log(n) $n\mathrm{log}\left(n\right)$ and the memory requirements are O(n)$\mathit{O}\left(n\right)$; thus the computation remains efficient and practical for very fine discretizations (very large values of n$n$). Collocation methods will be faster for small problems, but the method of nag_ode_bvp_ps_lin_solve (d02ue) should be faster for larger discretizations.

## Example

```function nag_ode_bvp_ps_lin_solve_example
n = int64(16);
a = -pi/2;
b =  pi/2;

% Set up boundary condition on left side of domain
y = [a, b];
% Set up Dirichlet condition using exact solution at x=a.
bmat = zeros(2, 3);
bmat(1, 1:2) = [1, 1];
bmat(2, 1:2) = [1, 1];
bvec = [cos(a) - sin(a), cos(b) - sin(b)];

% Set up problem definition
f = [1, 2, 3];

% Set up solution grid
[x, ifail] = nag_ode_bvp_ps_lin_cgl_grid(n, a, b);

% Set up problem right hand sides for grid and transform
f0 = -2*sin(x) - 2*cos(x);
[c, ifail] = nag_ode_bvp_ps_lin_coeffs(n, f0);

% Solve in coefficient space
[bmat, f, uc, resid, ifail] = nag_ode_bvp_ps_lin_solve(n, a, b, c, bmat, y, bvec, f);

% Evaluate solution and derivatives on Chebyshev grid
u = zeros(17, 3);
for q=0:2
[u(:, q+1),  ifail] = nag_ode_bvp_ps_lin_cgl_vals(n, a, b, int64(q), uc(:, q+1));
end

% Print Solution
fprintf('\nNumerical solution U and its first two derivatives\n');
fprintf('      x          U          Ux         Uxx\n');
for i=1:17
fprintf('%10.4f %10.4f %10.4f %10.4f\n', x(i), u(i, :));
end
```
```

Numerical solution U and its first two derivatives
x          U          Ux         Uxx
-1.5708    -0.0000     1.0000     0.0000
-1.5406     0.0302     0.9995    -0.0302
-1.4512     0.1193     0.9929    -0.1193
-1.3061     0.2616     0.9652    -0.2616
-1.1107     0.4440     0.8960    -0.4440
-0.8727     0.6428     0.7661    -0.6428
-0.6011     0.8247     0.5656    -0.8247
-0.3064     0.9534     0.3017    -0.9534
-0.0000     1.0000    -0.0000    -1.0000
0.3064     0.9534    -0.3017    -0.9534
0.6011     0.8247    -0.5656    -0.8247
0.8727     0.6428    -0.7661    -0.6428
1.1107     0.4440    -0.8960    -0.4440
1.3061     0.2616    -0.9652    -0.2616
1.4512     0.1193    -0.9929    -0.1193
1.5406     0.0302    -0.9995    -0.0302
1.5708    -0.0000    -1.0000     0.0000

```
```function d02ue_example
n = int64(16);
a = -pi/2;
b =  pi/2;

% Set up boundary condition on left side of domain
y = [a, b];
% Set up Dirichlet condition using exact solution at x=a.
bmat = zeros(2, 3);
bmat(1, 1:2) = [1, 1];
bmat(2, 1:2) = [1, 1];
bvec = [cos(a) - sin(a), cos(b) - sin(b)];

% Set up problem definition
f = [1, 2, 3];

% Set up solution grid
[x, ifail] = d02uc(n, a, b);

% Set up problem right hand sides for grid and transform
f0 = -2*sin(x) - 2*cos(x);
[c, ifail] = d02ua(n, f0);

% Solve in coefficient space
[bmat, f, uc, resid, ifail] = d02ue(n, a, b, c, bmat, y, bvec, f);

% Evaluate solution and derivatives on Chebyshev grid
u = zeros(17, 3);
for q=0:2
[u(:, q+1),  ifail] = d02ub(n, a, b, int64(q), uc(:, q+1));
end

% Print Solution
fprintf('\nNumerical solution U and its first two derivatives\n');
fprintf('      x          U          Ux         Uxx\n');
for i=1:17
fprintf('%10.4f %10.4f %10.4f %10.4f\n', x(i), u(i, :));
end
```
```

Numerical solution U and its first two derivatives
x          U          Ux         Uxx
-1.5708    -0.0000     1.0000     0.0000
-1.5406     0.0302     0.9995    -0.0302
-1.4512     0.1193     0.9929    -0.1193
-1.3061     0.2616     0.9652    -0.2616
-1.1107     0.4440     0.8960    -0.4440
-0.8727     0.6428     0.7661    -0.6428
-0.6011     0.8247     0.5656    -0.8247
-0.3064     0.9534     0.3017    -0.9534
-0.0000     1.0000    -0.0000    -1.0000
0.3064     0.9534    -0.3017    -0.9534
0.6011     0.8247    -0.5656    -0.8247
0.8727     0.6428    -0.7661    -0.6428
1.1107     0.4440    -0.8960    -0.4440
1.3061     0.2616    -0.9652    -0.2616
1.4512     0.1193    -0.9929    -0.1193
1.5406     0.0302    -0.9995    -0.0302
1.5708    -0.0000    -1.0000     0.0000

```