这个问题我之前在谱方法计算中也踩过坑——roots函数在处理高次、系数量级差异极大的多项式时，数值稳定性会直接崩盘，尤其是次数超过20左右时，浮点数精度的限制会让伴随矩阵特征值求解过程引入虚假的虚部。本质原因是：roots依赖伴随矩阵的特征值来求根，而伴随矩阵的构造会放大系数间的量级差异，导致数值误差急剧积累。

针对GLL节点的求解，完全不需要通过构造多项式系数再调用roots的方式，因为GLL节点有明确的数学性质，我们可以用更稳定的迭代方法直接求解。

核心思路 #

GLL节点（共$N$个）包含端点$x=-1$和$x=1$，内部节点是$N-2$次勒让德多项式导数$P_{N-1}'(x)$的根。我们可以：

1. 用Gauss节点作为内部节点的初始迭代值（已经足够接近真实根）；
2. 通过牛顿迭代法，直接针对$P_{N-1}'(x)=0$进行求解，迭代过程中用勒让德多项式的递推公式计算函数值和导数值，避免构造系数。

实现代码（Matlab） #

下面是一个完整的GLL节点生成函数，完全避开roots，数值稳定性拉满：

function [x_nodes, weights] = generate_gll_nodes(N)
    % 生成N个Gauss-Legendre-Lobatto节点及对应权重
    if N == 1
        x_nodes = 0;
        weights = 2;
        return;
    elseif N == 2
        x_nodes = [-1, 1];
        weights = [1, 1];
        return;
    end

    % 1. 生成初始迭代值：用N-2个Gauss-Legendre节点作为内部节点起点
    [gauss_x, ~] = get_gauss_legendre_nodes(N-2);
    x_guess = [-1; gauss_x; 1];

    % 2. 牛顿迭代求解内部节点（端点保持不动）
    tol = 1e-14;
    max_iter = 10;
    for iter = 1:max_iter
        % 计算P_{N-1}(x)和其一阶导数
        [P, P_prime] = legendre_poly_and_deriv(N-1, x_guess);
        % 内部节点满足P_{N-1}'(x)=0，计算迭代增量
        f = P_prime(2:end-1);
        f_prime = legendre_poly_second_deriv(N-1, x_guess(2:end-1));
        
        delta = f ./ f_prime;
        x_guess(2:end-1) = x_guess(2:end-1) - delta;

        if max(abs(delta)) < tol
            break;
        end
    end

    x_nodes = x_guess;
    % 计算GLL权重（可选，根据需求保留）
    [P, ~] = legendre_poly_and_deriv(N-1, x_nodes);
    weights = 2 ./ ((N*(N-1)) .* P.^2);
end

function [P, P_prime] = legendre_poly_and_deriv(n, x)
    % 用递推公式计算n次勒让德多项式及其一阶导数
    if n == 0
        P = ones(size(x));
        P_prime = zeros(size(x));
        return;
    elseif n == 1
        P = x;
        P_prime = ones(size(x));
        return;
    end

    P_prev2 = ones(size(x)); % P_0(x)
    P_prev1 = x;             % P_1(x)
    P_p_prev2 = zeros(size(x)); % P_0'(x)
    P_p_prev1 = ones(size(x));  % P_1'(x)

    for k = 2:n
        P = ((2*k-1)*x.*P_prev1 - (k-1)*P_prev2)/k;
        P_prime = P_p_prev2 + (2*k-1)*P_prev1;
        % 更新递推变量
        P_prev2 = P_prev1;
        P_prev1 = P;
        P_p_prev2 = P_p_prev1;
        P_p_prev1 = P_prime;
    end
end

function P_double_prime = legendre_poly_second_deriv(n, x)
    % 计算n次勒让德多项式的二阶导数
    [P, P_prime] = legendre_poly_and_deriv(n, x);
    mask = abs(x) ~= 1;
    P_double_prime = zeros(size(x));
    
    % 非端点用递推公式：(1-x²)P'' = -2xP' + n(n+1)P
    P_double_prime(mask) = (-2*x(mask).*P_prime(mask) + n*(n+1)*P(mask))./(1 - x(mask).^2);
    
    % 端点处的二阶导数用解析公式
    idx_pos = x == 1;
    if any(idx_pos)
        P_double_prime(idx_pos) = n*(n+1)*(n+2)*(n-1)/6;
    end
    idx_neg = x == -1;
    if any(idx_neg)
        P_double_prime(idx_neg) = (-1)^(n-1)*n*(n+1)*(n+2)*(n-1)/6;
    end
end

function [x, w] = get_gauss_legendre_nodes(n)
    % 生成n个Gauss-Legendre节点（用于初始迭代值）
    if n == 0
        x = []; w = []; return;
    end
    beta = 0.5./sqrt(1-(2*(1:n-1)).^(-2));
    T = diag(beta, 1) + diag(beta, -1);
    [V, D] = eig(T);
    x = diag(D);
    [x, idx] = sort(x);
    w = 2*V(1, idx).^2;
end

为什么这个方法更可靠？ #

1. 避免了系数量级爆炸：递推计算勒让德多项式和导数，完全不需要构造高次多项式的系数，从根源上消除了系数极差带来的数值误差；
2. 迭代收敛快：Gauss节点作为初始值已经非常接近GLL的内部节点，牛顿迭代通常3-5次就能达到1e-14以上的精度；
3. 稳定性高：即使N取到100甚至更大，也能输出完全实的、准确的GLL节点。

使用示例 #

% 生成22个GLL节点
[x, w] = generate_gll_nodes(22);
disp(x); % 全部为实根，无虚部

内容的提问来源于stack exchange，提问作者David Lu