Need an example Legendre-Gauss-Radau pseudospectral differentiation matrix or Matlab code

Question

I'm trying to implement various kinds of pseudospectral methods for direct optimization in Matlab using IPOPT. I've got some working Legendre-Gauss-Lobatto code, but would like to use the flipped Radau method. What greatly helped me out writing the LGL code was having an example 6x6 differentiation matrix to test against until my code could reproduce it. What I could use is a similar LGR differentiation matrix for e.g. the N=6 case (a 6x7 "unsparse" matrix due to the addition of the uncollocated -1 point).

I'd also like to know if the flipped Radau points are literally just inverted around 0, or if the algorithm to determine the LGR points needs to be modified to find $P_N(x)-P_{N+1}(x)$ instead of $P_N(x)+P_{N+1}(x)$?

Alternatively pointers to any Matlab or Python code which can reproduce these differentiation matrices would be appreciated.

For background on LGR/LGL/LG methods:

Divya Garg, Michael Patterson, William Hager, Anil Rao, David Benson, et al.. An overview of three pseudospectral methods for the numerical solution of optimal control problems. 2017. hal-01615132

Answer 1

This is probably too late to help you, but here is a compact code that will do points, weights and first derivatives for Gauss, Lobatto or either Radau.

function [x,w,A] = OCnonsymGLReig(n,meth)
   % code for nonsymmetric orthogonal collocation applications on 0 < x < 1
   % n - interior points
   % meth = 1,2,3,4 for Gauss, Lobatto, Radau (right), Radau (left)
   % x - collocation points
   % w - quadrature weights
   % A - 1st derivative
   na = [1 0 0 1];  nb = [1 0 1 0];  nt = n + 2;
   a = 1.0 - na(meth);  b = 1.0 - nb(meth);
   ab = r_jacobi(n,a,b);   ab(2:n,2) = sqrt(ab(2:n,2));
   T = diag(ab(2:n,2),-1) + diag(ab(:,1)) + diag(ab(2:n,2),+1);
   x = eig(T);  x=sort(x);
   x=0.5*(x+1.0);  x = [0.0;x;1.0];
   xdif = x-x'+eye(nt);     dpx = prod(xdif,2);
   w = (x.^nb(meth)).*((1.0 .- x).^na(meth))./(dpx.*dpx); w = w/sum(w); % quadrature
   A = dpx./(dpx'.*xdif);  A(1:nt+1:nt*nt) = 1.0 - sum(A,2);            % derivative
end

It relies on Gautschi's r_jacobi routine to calculate the recurrent coefficients. It is available at his site or mine. There are other routines you might find useful at my site - http://TildenTechnologies.com/Numerics. This is compact, but not very efficient. I am in the process of updating these routines for some that are more efficient, so stay tuned.

The code uses nodes at both ends, regardless of the quadrature selected. This is what you want for a BVP (see my article - Orthogonal Collocation Revisited (Mar. 2019), https://doi.org/10.1016/j.cma.2018.10.019). Many people get this wrong for Gauss points. Not sure what you need, but it would not be difficult to modify the derivative calculation.

It should work on Matlab, but I have tested it only with Octave.