Accuracy of the eigenvalue solutions#

Every eigenpair returned by lb() and freq() passed the verification described in Eigenvalue solvers: a relative residual not larger than check_rtol and, for symmetric matrices, the inertia check, which tells that no eigenvalue below the highest one returned is missing. The parameters below control what is computed and how strictly it is verified.

Repeated and clustered eigenvalues#

Structures with symmetries may have repeated eigenvalues, e.g. the two equal critical load multipliers of the cylinder model of the example below. The eigenvalue is then unique, but its eigenvectors are any basis of a subspace: two solvers, or two runs with a different num_eigvalues, return the same eigenvalues and possibly different eigenvectors, all of them valid. Eigenvectors of repeated or very close eigenvalues are compared by their subspaces, e.g. with scipy.linalg.subspace_angles(), and not one by one.

Number of eigenvalues#

num_eigvalues eigenvalues are computed by the sparse solvers. The inertia check guarantees that none is missing below the highest one returned, but a group of repeated eigenvalues may be cut at the last returned one, e.g. num_eigvalues=1 returns one mode of a pair. A margin above the number of modes needed avoids cutting a group whose modes are all needed, e.g. for a modal superposition or a comparison of mode shapes.

Tolerances#

  • check_rtol: maximum relative residual of an eigenpair, default 1e-3 in lb() and 1e-6 in freq(). A smaller value rejects less accurate sparse solutions, at the cost of more dense fallbacks or a RuntimeError when ARPACK does not reach it.

  • tol: tolerance of ARPACK, 0 for the machine precision.

Dense solutions#

  • sparse_solver=False solves the full problem with a dense solver and returns all the eigenvalues. Its cost grows with the cube of the number of dofs; it is a reference for validations.

  • max_dense_size: largest problem solved with the dense fallback when the sparse solvers failed. 0 disables the fallback, such that a failure of the sparse solvers raises a RuntimeError, and warnings.simplefilter('error', DenseFallbackWarning) turns the fallback into an error, see Eigenvalue solvers.

Examples#

Two equal critical load multipliers of a cylinder model of panels, whose matrices are saved in tests/data:

def test_lb_repeated_load_multipliers():
    """Two equal critical load multipliers: the sparse and the dense solvers
    agree on the load multipliers and on the eigenspace of the pair, but
    not necessarily on the individual eigenvectors, any basis of the
    eigenspace being valid. With num_eigvalues=1 only one of them is
    returned"""
    from scipy.linalg import subspace_angles

    K, KG = load_saved_matrices('cylinder_5panels')
    eigvals, eigvecs = lb(K, KG, silent=True, num_eigvalues=4)
    eig_dense, eigvecs_dense = lb(K, KG, silent=True, sparse_solver=False,
                                  num_eigvalues=4)
    np.testing.assert_allclose(eigvals, eig_dense[:4], rtol=1e-7)
    np.testing.assert_allclose(eigvals[0], eigvals[1], rtol=1e-7)
    assert np.all(relative_residuals(K, KG, eigvals, eigvecs) < 1e-3)
    # same eigenspace of the pair
    assert subspace_angles(eigvecs[:, :2], eigvecs_dense[:, :2]).max() < 1e-4
    # the individual eigenvectors of both solvers are equally valid
    assert np.all(relative_residuals(K, KG, eig_dense[:2],
                                     eigvecs_dense[:, :2]) < 1e-8)

    eigvals1, eigvecs1 = lb(K, KG, silent=True, num_eigvalues=1)
    assert eigvals1.shape == (1,)
    assert eigvecs1.shape == (K.shape[0], 1)
    np.testing.assert_allclose(eigvals1[0], eigvals[0], rtol=1e-7)

The verification detects wrong eigenpairs and missing load multipliers:

def test_lb_check_eigenpairs():
    """The verification detects wrong eigenpairs and a missing lowest load
    multiplier"""
    from structsolve.linear_buckling import _check_eigenpairs, _eigh_condensed

    K, KG = load_saved_matrices('plate_ssss_Nxx')
    mu, eigvecs = _eigh_condensed(K, KG, 10)
    eigvals = -1./mu
    assert _check_eigenpairs(K, KG, eigvals, eigvecs, 1e-3) is None

    # eigenvectors of other eigenvalues
    error = _check_eigenpairs(K, KG, eigvals, eigvecs[:, ::-1], 1e-3)
    assert 'relative residual' in error

    # wrong eigenvalues
    error = _check_eigenpairs(K, KG, eigvals*1.01, eigvecs, 1e-3)
    assert 'relative residual' in error

    # lowest load multiplier missing
    error = _check_eigenpairs(K, KG, eigvals[1:], eigvecs[:, 1:], 1e-3)
    assert 'missing' in error

    # both equal critical load multipliers missing
    K, KG = load_saved_matrices('cylinder_5panels')
    mu, eigvecs = _eigh_condensed(K, KG, 10)
    eigvals = -1./mu
    assert _check_eigenpairs(K, KG, eigvals, eigvecs, 1e-3) is None
    error = _check_eigenpairs(K, KG, eigvals[2:], eigvecs[:, 2:], 1e-3)
    assert 'missing' in error

The sparse solver of freq() against the dense solution:

@pytest.mark.parametrize('num_eigvalues', [5, 25, 40])
@pytest.mark.parametrize('symmetric', [True, False])
def test_freq_lowest_frequencies(symmetric, num_eigvalues):
    """The sparse solver, without dense fallback, returns the lowest
    frequencies of the full dense solution, also for many modes, which
    eigsh() computes in several passes when ncv is capped"""
    K, M = freq_problem(symmetric)
    lambda2, eigvecs = freq(K, M, silent=True, symmetric=symmetric,
                            num_eigvalues=num_eigvalues, max_dense_size=0,
                            sort=False)
    ref, _ = freq(K, M, silent=True, symmetric=symmetric,
                  sparse_solver=False, sort=False)
    w_ref = np.sort((-ref).real)[:num_eigvalues]
    assert lambda2.shape == (num_eigvalues,)
    np.testing.assert_allclose(np.sort((-lambda2).real), w_ref, rtol=1e-8)
    if symmetric:
        assert np.isrealobj(lambda2)