# basisChebyshev class

Defines a class to represent a univariate Chebyshev basis.

Objects of this class are of type 'handle', useful for passing the basis by reference.

Objects of class basisChebyshev have the following properties:

• n: scalar, number of nodes
• a: scalar, lower bound
• b: scalar, upper bound
• nodes: d.1 vector, basis nodes
• D: operators to compute derivative
• I: operators to compute integrals
• nodetype: node type, i.e. 'gaussian','lobatto', or 'endpoint'
• varname: string, variable name. Defaults to 'V0'.
• WarnOutOfBounds: boolean, warns if interpolating outside [a,b] if true. Defaults to 'false'.

Object of class basisChebyshev have the following methods:

• basisChebyshev: class constructor
• SetNodes: computes the nodes
• Interpolation: returns the interpolation matrix
• DiffOperator: computes operators for derivatives and integrals
• Interpolate: evaluates the interpolated function
• nearestNode: returns the basis node that is closest to input argument
• plot: plots the basis functions
• display: prints a summary of the basis
• copy: makes an independent copy
• CheckValidity: checks consistency of properties
• Reset: computes nodes and deletes D and I operators
• Setter methods for a, b, n, nodetype: change respective value and resets the basis.

Last updated: October 4, 2014.

## Contents

```classdef basisChebyshev < handle
properties (SetAccess = protected)
nodes       % nodes for own dimension
D           % operators to compute derivative
I           % operators to compute integrals
end

properties (SetAccess = public)
n  =  3     % number of nodes
a  = -inf     % lower limits
b  =  inf     % upper limits
nodetype = 'gaussian'   % type of nodes
varname = 'V0'             % optional name for variable
WarnOutOfBounds = true    % warns if interpolating outsite [a,b] interval
end

methods
```

## basisChebyshev

```        function B = basisChebyshev(...
n,...           scalar, number of nodes
a,...           scalar, lower bound
b,...           scalar, upper bound
nodetype,...    string, 'gaussian','endpoint', or 'lobatto'
varname)%       string, name for variable
```

Constructor for basisChebyshev

```B = basisChebyshev(n,a,b,nodetype,varname)
```

The following inputs are required:

• n: scalar, number of nodes
• a: scalar, lower bound
• b: scalar, upper bound

Optional inputs

• nodetype: node type, one of 'gaussian' (default), 'endpoint', or 'lobatto'
• varname: string, variable name. Defaults to 'V0'

If no inputs are provided (needed to preallocate array of basis)

```            if nargin ==0
B.n = 3;
B.a = -1;
B.b = 1;
return
end
```

SET PROPERTIES

```            B.n = n;
B.a = a;
B.b = b;
```

Set nodetype

```            if nargin > 3 && ~isempty(nodetype)
B.nodetype = lower(nodetype);
end
```

Set variable name

```            if nargin==5
B.varname = varname;
end
```

Check validity of inputs and set nodes

```            B.CheckValidity;
B.SetNodes;

end %basisChebyshev1
```

## copy

```        function Bcopy = copy(B)
```

Bcopy = B.copy

Makes a (value) copy of basis B. Since class basisChebyshev is a handle, in this example

```B2 = B;
B3 = B.copy;
```
```B2.a = 0;
B3.n = 9;
```

the line B2.a = 0 sets the lower limit a to zero in both B2 and B (B2 is just a pointer to B), whereas the line B3.n = 9 does not modify B (on assignment, B3 is simply a copy of B, not a pointer to B).

```            Bcopy = basisChebyshev(B.n,B.a,B.b,B.nodetype, B.varname);
end
```

## checkValidity

```        function CheckValidity(B)
```

B.CheckValidity

Checks the consistency of some basis B parameters. Called by constructor to ensure the following conditions :

• the basis is unidimensional,
• a < b,
• n >= 2,
• nodetype is either 'gaussian', 'endpoint', or 'lobatto'
```            if ~isscalar(B.a); error('Parameter ''a'' must be scalar'); end
if ~isscalar(B.b); error('Parameter ''b'' must be scalar'); end
if ~isscalar(B.n); error('Parameter ''n'' must be scalar'); end

if B.a > B.b;  error('Lower bound must be less than upper bound'), end
if B.n < 2;    error('n must be greater than one'), end

if ~ismember(lower(B.nodetype),{'gaussian' 'endpoint' 'lobatto'})
warning(strcat('Unknown nodetype ',B.nodetype,': using gaussian instead'))
B.nodetype = 'gaussian';
end

end %CheckValidity
```

## SetNodes

```        function SetNodes(B)
```

B.SetNodes

Computes the nodes of the basis B, depending on B.nodetype, which can be 'gaussian', 'endpoint', or 'lobatto'. This method is called by the constructor.

```            n = B.n;
a = B.a;
b = B.b;

s=(b-a)/2;
m=(b+a)/2;

switch B.nodetype
case 'gaussian'   % Gaussian nodes
k = pi*(0.5:(n-0.5))';
x = m-cos(k/n)*s;
case 'endpoint'     % Extend nodes to endpoints
k = pi*(0.5:(n-0.5))';
x = m-cos(k/n)*s;
aa = x(1);
bb = x(end);
x = (bb*a-aa*b)/(bb-aa)+(b-a)/(bb-aa)*x;
case 'lobatto'     % Lobatto nodes
k = pi*(0:(n-1))';
x = m - cos(k/(n-1))*s;
end

x(abs(x)<eps) = 0;
B.nodes = x;
end %SetNodes
```

## DiffOperator

```        function DiffOperator(B,...
orderD,...            order of derivative
integrate...          integral if true
)
```

B.DiffOperator(orderD,integrate)

Compute operators to differentiate and integrate. This method is needed by Interpolation when input order is not zero; it should not be called directly by the user.

It updates the operators in B.D (if integrate is false) or in B.I (otherwise). Input orderD is the maximum order to compute.

```            if nargin < 2, orderD = 1; end
if nargin < 3, integrate = false; else integrate = logical(integrate); end
if orderD < 0, error('In DiffOperator: order must be nonnegative'), end

derivative = ~integrate;

if orderD==0, return, end

% Use previously stored values for orders
if (derivative && length(B.D) >= orderD), return, end
if (integrate  && length(B.I) >= orderD), return, end

n = B.n;
a = B.a;
b = B.b;

if(n-orderD < 2 && derivative)
warnmsg = strcat('Insufficient nodes: truncating order to = ',num2str(n-2));
warning(warnmsg)
orderD = n-2;
end

if derivative
%===============TAKE DERIVATIVE=========
if isempty(B.D)
[j,i] = meshgrid(1:n,1:n);
rc = mod(i+j,2) & (j>i);
d = zeros(n,n);
d(rc) = (4/(b-a))*(j(rc)-1);
d(1,:) = d(1,:)/2;
d = sparse(d(1:end-1,:));
B.D{1} = d;
else
d = B.D{1};
end

h = length(B.D) + 1;  % index of 1st missing element in B.D
for ii = h:orderD
B.D{ii} = d(1:n-ii,1:n-ii+1) * B.D{ii-1};
end
else
%===============TAKE INTEGRAL=========
nn = n + orderD;
i = (0.25*(b-a))./(1:nn);
d = sparse(diag(i) + diag(-i(1:nn-2),2));  %good
d(1,1) = 2*d(1,1);

d0 = ((-1).^(0:nn-1)).*sum(d);

if isempty(B.I)
B.I{1} = [d0(1:n);d(1:n,1:n)];
end
h = length(B.I)+1;    % index of 1st missing element in B.I

for ii = h:orderD
B.I{ii}=[d0(1:n+ii-1);d(1:n+ii-1,1:n+ii-1)] * B.I{ii-1};
end

end
end
```

## Interpolation

```        function Phi = Interpolation(B,... ##CompEcon:  chebbas.m
x,...                      k-vector of the evaluation points
order,...                  order of derivative
integrate)%                integrate if true
```
```Phi = B.Interpolation(x,order,integrate)
```

Computes interpolation matrices Phi for basis B.

Its inputs are

• x, k.1 vector of evaluation points. Defaults to basis nodes.
• order, h.1 vector, for order of derivatives. Defaults to 0.
• integrate, boolean, integrate if true, derivative if false. Defaults to false

Output Phi returns the interpolation k.n.h array, where n is number of basis functions.

```            if nargin<2 || isempty(x)
x = B.nodes;  % use basis nodes if x is missing
hasArg_x = false;
nx = B.n;
else
x = x(:);   % make sure x is column vector;
hasArg_x = true;
nx = numel(x);
end

if nargin<3, order = 0; end
if nargin<4, integrate = false; end

maxorder = max(order);
nn = B.n + maxorder * integrate;

% check if out of bounds
if hasArg_x && B.WarnOutOfBounds
isBelow =  x < B.a - eps(B.a);
isAbove =  x > B.b + eps(B.b);
isOut = isBelow | isAbove;

if any(isOut)
xout = x(isOut);
error(['OutOfBounds!: Interpolating %s = %f failed because basis defined only in [%4.2f, %4.2f]\n',...
'To allow interpolation outside this interval, set basis field ''WarnOutOfBounds'' to false'],...
B.varname,xout(1),B.a,B.b)
end

end

% Compute order 0 interpolation matrix
if hasArg_x || ~strcmp(B.nodetype,'gaussian') % evaluate at arbitrary nodes
m = length(x);
z = (2/(B.b - B.a))* (x-(B.a + B.b)/2);
bas = zeros(m,nn);
bas(:,1) = 1;
bas(:,2) = z;
z = 2*z;
for i = 3:nn
bas(:,i) = z.*bas(:,i-1) - bas(:,i-2);
end
else  % evaluate at usual Gaussian nodes
temp=((B.n-0.5):-1:0.5)';
bas  = cos((pi/B.n) * temp *(0:nn-1));
end

% Compute Phi ////////////////////////
%Phi = cell(length(order),1);
Phi = zeros(nx,B.n,length(order));

done = nan(length(order),1);

for ii = 1:length(order)
%Phi{ii} = Phi{find(order(ii)==done,1)};
Phi(:,:,ii) = Phi(:,:,find(order(ii)==done,1));
else
```

The operators I and D are saved into the basis when created. The try/catch sequence avoids calling the DiffOperator method (to speed up execution)

```                    if order(ii)==0
%Phi{ii} = bas;
Phi(:,:,ii) = bas;
elseif integrate
%Phi{ii} = bas(:,1:B.n+order(ii)) * B.I{order(ii)};
try
Phi(:,:,ii) = bas(:,1:B.n+order(ii)) * B.I{order(ii)};
catch
DiffOperator(B,maxorder,integrate);
Phi(:,:,ii) = bas(:,1:B.n+order(ii)) * B.I{order(ii)};
end
else % derivative
%Phi{ii} = bas(:,1:B.n-order(ii)) * B.D{order(ii)};
try
Phi(:,:,ii) = bas(:,1:B.n-order(ii)) * B.D{order(ii)};
catch
DiffOperator(B,maxorder,integrate);
Phi(:,:,ii) = bas(:,1:B.n-order(ii)) * B.D{order(ii)};
end
end
Phi(abs(Phi)<eps(10)) = 0;

end
done(ii) = order(ii);

end

end %Interpolation
```

## Evaluate

```        function y = Interpolate(B,...      ##CompEcon:  chebbas.m
coef,...             coefficient matrix
varargin...          inputs for Interpolation method
)
```
```y = B.Evaluate(coef,x,order,integrate)
```

Evaluates the interpolated functions defined by basis B and coefficients coef. Inputs x, order and integrate are as required by Interpolation method, while coef is a n.m matrix (n basis functions and m interpolated functions)

Output y returns the interpolated functions as a k.m.h array, where k is the number of evaluation points, m the number of functions (number of columns in coef), and h is number of order derivatives.

```            if nargin<2, error('Missing ''coef'' input');end

Phi = B.Interpolation(varargin{:});

nx = size(Phi,1);  % number of evaluation points
ny = size(coef,2); % number of evaluated functions
no = size(Phi,3);  % number of order evaluations

y = zeros(nx,ny,no);

for h = 1:no
y(:,:,h) = Phi(:,:,h) * coef;
end

if ny==1  % only one function
y = squeeze(y);
end

end %Evaluate
```

## plot

```         function plot(B,ii)
```
```B.plot(ii)
```

Plots the basis functions. Input ii is optional, to plot only the specified basis functions.

```            if nargin<2
ii = 1:B.n;
else
ii(ii>B.n) = [];
end
x = linspace(B.a,B.b,1+min(120,10*B.n))';
Phi = B.Interpolation(x);
Phi = Phi{1};
hold off
plot(x,0*x,'k'); hold on
plot(x,Phi(:,ii),'LineWidth',1.75)
plot(B.nodes,0*B.nodes,'ro')
title(B.varname)
hold off

end
```

## Nearest node

```        function [xnode,idx] = nearestNode(B,x)
```
```[xnode,idx] = nearestNode(B,x)
```

For a given value x, it finds the nearest node in basis B. This method is useful in simulations, when an initial guess solution at x can be given by xnode. idx is the associated index, e.g.

```xnode = B.nodes(idx,:)'
```
```            gap = abs(B.nodes - x);
[~, idx] = min(gap);
xnode = B.nodes(idx);

end
```

## display

```        function display(B)
```

Prints a summary of the basis

```            for i=1:numel(B)
fprintf('\nChebyshev Basis of one dimension with %s nodes\n\n',B(i).nodetype)

fprintf('\t%-22s %-12s %-16s\n','    Variable','# of nodes','    Interval')
fprintf('\t%-22s     %-8.0f [%6.2f, %6.2f]\n',B(i).varname,B(i).n,B(i).a,B(i).b)
fprintf('\n\n')
end
end
```

## Reset

```        function Reset(B)
```

Following a change in the basis definition, the nodes must be recomputed and the fields D and I are no longer valid. This method is called by the setter methods, so the user does not need to call it explicitly.

```            B.CheckValidity;
B.SetNodes;
B.I = [];
B.D = [];
end
```

## Setter methods

The following setter methods allow changing the basis definition (parameters n,a,b,nodetype). Changes to these parameters require recomputing the nodes and deleting the operators D and I.

```        function set.a(B,value)
B.a = value;
B.Reset;
end
```

```        function set.b(B,value)
B.b = value;
B.Reset;
end
```

```        function set.n(B,value)
B.n = value;
B.Reset;
end
```

```        function set.nodetype(B,value)
```    end %methods