The variational determination of the two-fermion reduced density matrix is described for harmonically trapped, ultracold few-fermion systems in one dimension with equal spin populations. This is accomplished by formulating the problem as a semi-definite program, with the two-fermion reduced density matrix being subject to well-known N-representability conditions. The ground-state energies, as well as the density, pair-correlation function, and lower-order eigenvalues of the two fermion reduced density matrix of various fermionic systems are found by utilising an augmented Lagrangian method for semi-definite programming. The ground-state energies are found to match well to those determined by full-configuration interaction and coupled-cluster calculations and the density, pair-correlation function, and eigenvalue results demonstrate that the salient features of these systems are well-described by this method. These results collectively demonstrate the utility of the reduced density matrix method firstly in describing strong correlation arising from short-range interactions, suggesting that the well-known N-representability conditions are sufficient to model ultracold fermionic systems, and secondly in illustrating the prospect of treating larger systems currently out of the reach of established methods.