一个可能的解决方案（@gammatester 建议）是使用 Jacobi 多项式。这避免了通过“朴素”多项式评估添加大多项式系数时的灾难性抵消问题。

径向 Zernike 多项式可以用 Jacobi 多项式表示如下（见方程（6））

$$R^m_n(\rho) = (-1)^{(n-m)/2}\rho^m \cdot P^{(m,0)}_{(n-m)/2} \Big(1-2\rho^2 \Big)$$

然而，在 MATLAB 中，jacobiP(n,a,b,x)对于x=rho. 该jacobiP函数实际上是符号工具箱的一部分，多项式的计算被推迟到符号引擎，它以速度换取任意精度。因此需要手动执行 Jacobi 多项式。

由于 Jacobi 函数的参数都是非负的 ($\alpha=m$,$\beta=0$,$n^*=(n-m/2)$)，我们可以使用以下表达式（参见Wikipedia，注意我填写的值$s$)

$$\begin{multline} P_n^{(\alpha,\beta)}(\rho) = (n+\alpha)!(n+\beta)! \quad \cdot\\ \sum_{s={0}}^n \left[\frac{1}{s!(n+\alpha-s)!(\beta+s)!(n-s)!} \left(\frac{x-1}{2}\right)^{n-s}\left(\frac{x+1}{2}\right)^{s} \right]\end{multline}$$

在 MATLAB 中，这转换为 (Jacobi p olice d epartment P olynomial, ' D ouble' implementation)

function P = jacobiPD(n,a,b,x)
    P = 0;
    for  s  0:n
        P = P + ...
            1/(factorial(s)*factorial(n+a-s)*factorial(b+s)*factorial(n-s)) * ...
            ((x-1)/2).^(n-s).*((x+1)/2).^s;
    end
    P = P*factorial(n+a) * factorial(n+b);
end

因此，实际的径向 Zernike 多项式是 (对于m=abs(m))

Z = @(m,n,rho) (-1)^((n-m)/2) * rho.^m .* jacobiPD((n-m)/2,m,0,1-2*rho.^2);

---

注意：这个自我回答只是一个实际的解决方案；随意标记另一个解释为什么它有效的答案。

在本文中，Honarvar 和 Paramesran 推导出了一种有趣的方法，以非常好的递归方式计算径向 Zernike 多项式。递归公式非常简单，无需除以或乘以大整数：

$$R^m_n(\rho) = \rho \left(R^{|m-1|}_{n-1}(\rho)+R^{m+1}_{n-1}(\rho)\right) - R^{m}_{n-2}(\rho)$$ 我建议看一下 Honarvar 和 Paramesran 论文中的图 1，它清楚地说明了不同 Zernike 多项式之间的依赖关系。

这是在以下 Octave 脚本中实现的：

clear                                     % Tested with Octave instead of Matlab
N = 120;
n_r = 1000;
R = cell(N+1,N+1);
rho = [0:n_r]/n_r;
rho_x_2 = 2*[0:n_r]/n_r;

R{0+1,0+1} = ones(1,n_r+1);               % R^0_0  Unfortunately zero based cell indexing is not possible
R{1+1,1+1} = R{0+1,0+1}.*rho;             % R^1_1  ==>  R{...+1,...+1} etc.
for n = 2:N,
    if bitget(n,1) == 0,                  % n is even
        R{0+1,n+1} = -R{0+1,n-2+1}+rho_x_2.*R{1+1,n-1+1};                % R^0_n
        m_lo = 2;
        m_hi = n-2;
    else
        m_lo = 1;
        m_hi = n-1;
    end
    for m = m_lo:2:m_hi,
        R{m+1,n+1} = rho.*(R{m-1+1,n-1+1}+R{m+1+1,n-1+1})-R{m+1,n-2+1};  % R^m_n
    end
    R{n+1,n+1} = rho.*R{n-1+1,n-1+1};                                    % R^n_n
end;


Z = @(m,n,rho) (-1)^((n-m)/2) * rho.^m .* jacobiPD((n-m)/2,m,0,1-2*rho.^2);
m = 22;
n = 112;
figure
plot(rho,Z(m,n,rho))
hold on
plot(rho,R{m+1,n+1},'r');
xlabel("rho")
ylabel("R^{22}_{112}(rho)")
legend("via Jacobi","recursive");
%print -djpg plt.jpg

m = 0;
n = 46;
max_diff_m_0_n_46 = norm(Z(m,n,rho)-R{m+1,n+1},inf)

例如，此代码生成的图显示，$m = 22$， 和$n = 112$，灾难性的取消发生在附近$\rho = 0.7$，如果 Zernike 径向多项式是通过 Jacobi 多项式计算的。因此，人们还必须担心低次泽尼克多项式的准确性。

递归方法似乎更适合以稳定的方式计算这些高阶 Zernike 多项式。尽管如此，对于$m = 0$和$n = 46$，雅可比和递归方法之间的最大差异是（仅？）1.4e-10，这对于您的应用程序可能足够准确。