In this paper, a Nonlinear Program (NLP) matrix generation method based on the dual number theory for solving pseudospectral trajectory optimization problems is proposed. Dual numbers have the capacity to compute differentials with precision, which is utilized to determine the relationship between variables and constraints. In comparison to the traditional method using Not-a-Number (NaN), the proposed method has the advantage of simultaneously determining whether the relationship is linear or not. Using the proposed method, the relationship is classified into three categories: independent, linearly-correlated, and nonlinearly-correlated. Subsequently, the independent or linearly-correlated relationships lead to zero or constant elements for the Jacobian matrix and zero elements for the Hessian matrix. Only the nonlinearly-correlated relationships result in non-constant elements that need to be calculated iteratively during the solving process. Compared to the NaN method, the proposed method reduces the number of non-constant elements, thereby reducing the computational burden. Simulations verify that the proposed method is superior in determining zero and constant elements in NLP matrices, thus improving the generation efficiency.