Remes算法matlab代码注释

以下是Remes算法的MATLAB代码注释：

function [p, q, err] = remes(x, y, m, n)
% REMES Remes algorithm for rational function approximation
%   [P, Q, ERR] = REMES(X, Y, M, N) computes the Pade approximant P/Q
%   of degree M/N to the function f(x) = y(x) on the nodes X, where X is
%   a vector of distinct real numbers. The output ERR is the maximum
%   absolute error over the nodes X.
%
%   See also RESIDUE, PADEN, POLYVAL, POLYFIT

% Initialization
v = ones(size(x));
for j = 1:(2*n)
    v = [v, v(:, j).*x'];
end
A = zeros(2*n+1, m+n+1);
A(:,1:n+1) = v(:,1:n+1);
for j = 1:n
    A(:,j+1) = -x'.*A(:,j);
end
B = zeros(2*n+1,1);
B(1:n+1) = y';
% Iteration
for k = 1:m
    r = B-A(:,n+k+1:m+n+k)*[B(n+k+1:m+n+k);zeros(n,1)];
    h = A(:,1:n+k)\r;
    % Compute error and check convergence
    p = h(1:k+1);
    q = [1, h(k+2:n+k+1)'];
    err = max(abs(polyval(p, x)./polyval(q, x)-y));
    if err < 1.0e-16
        break
    end
    % Update A and B
    for j = 1:n
        A(:,k+n+1+j) = -x'.*A(:,k+n+j);
    end
    A(:,k+n+2) = v(:,k+n+2);
    for j = k+2:n+k+1
        A(:,j) = A(:,j)-p(j-k)*A(:,k+n+2-j+k);
    end
    B(n+k+2:m+n+1) = [zeros(k-n-1,1);p(end:-1:2)];
    B(1:n+k+1) = y'-polyval(q(end:-1:2),x).*B(n+k+2:m+n+1);
end
end

代码主要分为两个部分，初始化和迭代。

初始化部分首先将 x 向量的幂次作为矩阵 v 的列，然后构造矩阵 A 和向量 B，将其用于最小二乘求解。

迭代部分首先求解出当前的有理函数 P/Q 并计算其在节点上的最大误差 err。如果 err 小于某个阈值，算法就终止；否则，更新 A 和 B，并继续迭代。在更新 A 和 B 的过程中，需要将新的向量和矩阵插入到 A 和 B 的末尾，然后用新的 P 更新 A 和 B 中的一部分元素。