Fractional derivatives (FDs) or derivatives of arbitrary order have been used in many applications, and it is envisioned that in the future they will appear in many functional minimization problems of practical interest. Since fractional derivatives have such properties as being non-local, it can be extremely challenging to find analytical solutions for fractional parametric optimization problems, and in many cases, analytical solutions may not exist. Therefore, it is of great importance to develop numerical methods for such problems. This paper presents a numerical scheme for a linear functional minimization problem that involves FD terms. The FD is defined in terms of the Riemann-Liouville definition; however, the scheme will also apply to Caputo derivatives, as well as other definitions of fractional derivatives. In this scheme, the spatial domain is discretized into several subdomains and 2-node one-dimensional linear elements are adopted to approximate the solution and its fractional derivative at point within the domain. The fractional optimization problem is converted to an eigenvalue problem, the solution of which leads to fractional orthogonal functions. Convergence study of the number of elements and error analysis of the results ensure that the algorithm yields stable results. Various fractional orders of derivative are considered, and as the order approaches the integer value of 1, the solution recovers the analytical result for the corresponding integer order problem.