Abstract: The adoption of the first-order local approximation functions has improved the accuracy, but it has also led to rank deficiency of the global stiffness matrix for the higher-order numerical manifold method, meaning the existence of the linear dependence. The first-order Taylor’s expansions with regard to the interpolation point are adopted as the local displacement functions, which makes the degrees of freedom defined on the physical cover have definite physical meanings. Then the first-order partial differential derivatives are expressed by the strain components, leading to the PC-u-ε. The strain components are further replaced by the stress components in the local framework, creating the PC-u-σ. In this way, both the displacement and the stress boundary conditions are easily applied. Numerical examples show that the PC-u-ε alone significantly causes a drastic decrease in rank deficiency, while deploying the PC-u-σ along the stress boundary completely eliminates the rank deficiency and retains higher accuracy.