Skip to content

Wrong coefficient for 1D-QuinticSpline #421

Description

@Jasmine969
Image

The normalization coefficient for 1D-QuinticSpline is $1/(60h)$ rather than $1/(120h)$.

Validation:

def quintic_spline(r, cut, dimension=3):
    h = cut / 3
    q = r / h
    w = np.where(q < 3, (3 - q) ** 5, 0)
    w -= np.where(q < 2, 6 * (2 - q) ** 5, 0)
    w += np.where(q < 1, 15 * (1 - q) ** 5, 0)
    if dimension == 1:
        alphaD = 1 / (60 * h)
    elif dimension == 2:
        alphaD = 7 / (478 * np.pi * h ** 2)
    else:  # dimension == 3
        alphaD = 1 / (120 * np.pi * h ** 3)
    return w * alphaD

kernels = {'Quintic-Spline': quintic_spline}

def test_kernel(name, cutoff):
    from scipy.integrate import quad, dblquad, tplquad
    res1d, _ = quad(lambda r: kernels[name](r, cut=cutoff, dimension=1),
                    0, cutoff)
    res2d, _ = dblquad(lambda r, theta: kernels[name](r, cut=cutoff, dimension=2) * r,
                       0, 2 * np.pi, 0, cutoff)
    res3d, _ = tplquad(lambda r, phi, theta: kernels[name](r, cut=cutoff, dimension=3) * r * r * np.sin(phi),
                       0, 2 * np.pi, 0, np.pi, 0, cutoff)
    flag_pass = True
    for i, res in enumerate([res1d, res2d, res3d]):
        print(f'Dimension {i+1}: result={res}.')
        if abs(res - 1) > 1e-6:
            print(f'Dimension {i + 1} does not satisfy the normalization condition!')
            flag_pass = False
    if flag_pass:
        print(f'Kernel {name} has passed the normalization condition test!')

if __name__ == '__main__':
    test_kernel('Quintic-Spline', 3)

Output:

Dimension 1: result=0.9999999999999204.
Dimension 2: result=0.9999999999972096.
Dimension 3: result=0.9999999999634721.
Kernel Quintic-Spline has passed the normalization condition test!

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions