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

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

• According to Matlab if x is a vector of the roots of $P_N(x)+P_{N+1}(x)$ for N=5 then legendreP(5,x) + legendreP(6,x) is the zero vector and legendreP(5,-x) - legendreP(6,-x) is also the zero vector, so I think that answers my question about generating the flipped LGR points. – lamont Jun 19 '19 at 22:03
• I think also the obvious thing for me to do is to just take a function that is trivial to differentiate and check applying the derivative operator against it vs. the known analytic solution. Already did that with my known-working LGL code. Can then invert that to test integration as well.. – lamont Jun 20 '19 at 17:15
• The matlab code associate with the book SPECTRAL METHODS: Algorithms, Analysis and Applications is here and legsrddiff.m by changing dy=dy1+dy2; y=y1+y2; to dy=dy1-dy2; y=y1-y2; -- still working on adding the noncollocated point. – lamont Jun 25 '19 at 15:33

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
% x - collocation points
% 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.

• Actually that's incredibly timely. I did stumble through it myself, but went back and for simplicity was using LGL for a few months. Recently hit a problem where I have a suspicion that LGR might be more stable, so I'm looking at adding that to my solver, so I can immediately use this! – lamont Nov 15 '19 at 19:29
• And yes, I need the non-collated endpoints and need to use quadrature for those for the purposes of embedding in a BVP solver. – lamont Nov 15 '19 at 19:33