Chebfun `eigs` 在做什么

计算科学 matlab pde 数值分析 特征值 切比雪夫
2021-12-25 13:33:42

这是在做什么?看起来像原来的特征值问题被转换为具有不同维度搭配点的广义特征值问题。有人可以对此进行更多解释吗?

eigs.m

% Linear combination coefficients for convergence test. The convergence of the
% combination is the same as the worst constituent function. The nontrivial
% coefficents are to make accidental cancellations extremely unlikely.
coeff = 1./(2*(1:k)');
for dim = [dimVals NaN]
    [V, D, P] = getEigenvalues(discA, discB, k, sigma);

    % Combine the eigenfunctions into a composite.
    v = V*coeff(1:size(V,2));

    % Convert the different components into cells
    u = partition(discA, P*v);

    % Test the happiness of the function pieces:
    vscale = zeros(1, sum(isFun));   % intrinsic scaling only
    [isDone, cutoff] = testConvergence(discA, u(isFun), vscale, prefs);

    if ( all(isDone) )
        break
    elseif ( ~isnan(dim) )
        % Update the discretiztion dimension on unhappy pieces:
        discA.dimension(~isDone) = dim;
    end
end

eigs.m/getEigenvalues

function [V, D, P] = getEigenvalues(discA, discB, k, sigma)
% Formulate the discrete problem and solve for the eigenvalues

    % Discretize the LHS operator (incl. constraints/continuity):
    [PA, P, C, ignored, PS] = matrix(discA);

    % Discretize the RHS operator, or use identity.
    if ( ~isempty(discB) )
        % TODO: This is untidy. Can we make a method to do this? NH Apr 2014.
        discB.dimension = discA.dimension;
        PB = matrix(discB);
        % Project RHS matrix and prepend rows for the LHS constraints.
        PB = [ zeros(size(C)) ; PB ];
    else
        PB = [ zeros(size(C)) ; PS ];
    end

    % Compute eigenvalues.
    if ( length(PA) <= 2000 )
        [V, D] = eig(full(PA), full(PB));
        lam = diag(D);

        % Find the ones we're looking for.
        lam = deflate(lam, size(C,1));
        idx = nearest(lam, sigma);
        idx = filter(idx, P*V, k, discA);

        % Extract them:
        V = V(:,idx);
        D = D(idx,idx);

    else
        % TODO: Experimental.
        [V, D] = eigs(PA, PB, k, sigma);
    end

end

有关完整代码,请参见https://github.com/chebfun/chebfun/blob/development/%40linop/eigs.m

0个回答
没有发现任何回复~