我正在尝试使用IPOPT在Matlab中实现各种伪谱方法直接优化。我有一些有效的 Legendre-Gauss-Lobatto 代码，但想使用翻转的 Radau 方法。对我编写 LGL 代码有很大帮助的是有一个示例 6x6 微分矩阵进行测试，直到我的代码可以重现它。我可以使用类似的 LGR 微分矩阵，例如 N=6 的情况（由于添加了未并置的 -1 点，所以一个 6x7“非稀疏”矩阵）。
我还想知道翻转的 Radau 点是否实际上只是在 0 左右反转，或者是否需要修改确定 LGR 点的算法以找到而不是 ?
或者指向可以重现这些微分矩阵的任何 Matlab 或 Python 代码的指针将不胜感激。
有关 LGR/LGL/LG 方法的背景：
这可能为时已晚，无法帮助您，但这里有一个紧凑的代码，它将为 Gauss、Lobatto 或 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
它依赖 Gautschi 的 r_jacobi 例程来计算循环系数。它可以在他的网站或我的网站上找到。在我的网站http://TildenTechnologies.com/Numerics上您可能会发现其他一些有用的例程。这是紧凑的，但不是很有效。我正在为一些更高效的例程更新这些例程，敬请期待。
代码使用两端的节点，而不管选择的正交。这就是您想要的 BVP（请参阅我的文章 - 重新审视正交搭配（2019 年 3 月），https://doi.org/10.1016/j.cma.2018.10.019）。许多人对高斯点都误解了这一点。不确定您需要什么，但修改导数计算并不难。
它应该在 Matlab 上工作，但我只用 Octave 测试过它。