What is this doing? Looks like the original eigenvalue problem is converted into generalized eigenvalue problems with different dimensions of collocation points. Can someone explain more about this?


% 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) )
    elseif ( ~isnan(dim) )
        % Update the discretiztion dimension on unhappy pieces:
        discA.dimension(~isDone) = dim;


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 ];
        PB = [ zeros(size(C)) ; PS ];

    % 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);

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


For the full code see https://github.com/chebfun/chebfun/blob/development/%40linop/eigs.m

  • 5
    $\begingroup$ Have you read the docs and the paper they reference? mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/… $\endgroup$ – David Ketcheson Dec 28 '17 at 4:49
  • 2
    $\begingroup$ Could you be a little more specific? These kinds of questions are really tricky if you don't have access to the paper that directly describes the code, so it would really help to say what you already know about it, as well as what exactly you want to know. It sets up a matrix eigenvalue problem corresponding to a linear operator eigenvalue problem, so do you mean to ask how it does that, or why there are different numbers of points in a loop, or something else? $\endgroup$ – Kirill Dec 29 '17 at 3:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.