Convert geodetic coordinates to geocentric (cartesian)
I can't explain the discrepancy with the online conversion tool without seeing the source code for the latter.

1. Review


The comment for the ellipsoid parameters should give the units for the semi-major axis. Also, I think "reciprocal flattening" would be a better description of the second parameter, given that you use the name rf later on in the code.
The docstring for geodetic_to_geocentric should give the units (degrees) for the latitude and longitude parameters.
"Height above ellipsoid" would be a clearer description for the h parameter.
The parameters to geodetic_to_geocentric could be spelled out in full, which would help understand their meaning. Parameter names are documentation too.
Comments would help understand the meaning of values like N.
The expressions (N + h) * cos(lat), sin(lat), and (1 - 1 / rf) ** 2 are evaluated twice. These could be stored in local variables to avoid duplicated work.


2. Revised code

from math import cos, radians, sin, sqrt

# Ellipsoid parameters: semi major axis in metres, reciprocal flattening.
GRS80 = 6378137, 298.257222100882711
WGS84 = 6378137, 298.257223563

def geodetic_to_geocentric(ellipsoid, latitude, longitude, height):
    """Return geocentric (Cartesian) Coordinates x, y, z corresponding to
    the geodetic coordinates given by latitude and longitude (in
    degrees) and height above ellipsoid. The ellipsoid must be
    specified by a pair (semi-major axis, reciprocal flattening).

    """
    φ = radians(latitude)
    λ = radians(longitude)
    sin_φ = sin(φ)
    a, rf = ellipsoid           # semi-major axis, reciprocal flattening
    e2 = 1 - (1 - 1 / rf) ** 2  # eccentricity squared
    n = a / sqrt(1 - e2 * sin_φ ** 2) # prime vertical radius
    r = (n + height) * cos(φ)   # perpendicular distance from z axis
    x = r * cos(λ)
    y = r * sin(λ)
    z = (n * (1 - e2) + height) * sin_φ
    return x, y, z