**7. Numerical Results and Discussions**

In this section, we investigate the accuracy of the method by testing it on some numerical examples. We apply the numerical algorithm to two test problems using different orthogonal polynomials as a basis. The results for the test problems are shown through the figures and tables.

**Example 1.** *Consider a non-linear fractional variational problem as in Equation (1) with g*(*x*) = *h*(*x*) = <sup>1</sup> <sup>1</sup>+*x<sup>β</sup> ; we then have the following non-linear fractional variational problem [19]:*

$$J(\mathbf{y}) = \int\_0^1 \left( \frac{1}{1+\mathbf{x}^\delta} D^\mathbf{z} y(\mathbf{x}) - \left( I^{1-\mathbf{z}} y(\mathbf{x}) + 1 \right) \frac{\beta \mathbf{x}^{\delta-1}}{\left( 1+\mathbf{x}^\delta \right)^2} \right)^2 d\mathbf{x} \tag{56}$$

*under the constraints*

$$y(0) = 0, \; I^{1-\alpha}y(1) = \epsilon.$$

*The exact solution of the above equation is given as*

$$y\_{\text{exact}}(\mathbf{x}) = \left(\frac{1}{2}(1+\epsilon) - 1\right) \left(\frac{\Gamma(\beta+2)}{\Gamma(\beta+\alpha+1)}\mathbf{x}^{\beta+\alpha} + \frac{1}{\Gamma(\alpha+1)}\mathbf{x}^{\alpha}\right) + \frac{\Gamma(\beta+1)}{\Gamma(\alpha+\beta)}\mathbf{x}^{\beta+\alpha-1}.$$

We discuss this example for different values of *α* = 0.5, 0.6, 0.7, 0.8, 0.9, or 1, *β* = 5, and = 1. In Figures 1–5, it is shown that the solutions for the two different values of *α* = 0.8 and *α* = 1 coincide with the exact solutions for different orthogonal polynomials at *n* = 5.

**Figure 1.** Comparison of exact and numerical solutions using S1 for *α* = 0.8 and *α* = 1, Example 1.

**Figure 2.** Comparison of exact and numerical solutions using S2 for *α* = 0.8 and *α* = 1, Example 1.

**Figure 3.** Comparison of exact and numerical solutions using S3 for *α* = 0.8 and *α* = 1, Example 1.

**Figure 4.** Comparison of exact and numerical solutions using S4 for *α* = 0.8 and *α* = 1, Example 1.

**Figure 5.** Comparison of exact and numerical solutions using S5 for *α* = 0.8 and *α* = 1, Example 1.

In Figures 6–10, it is shown that the solution varies continuously for Shifted Legendre polynomials, Shifted Chebyshev polynomials of the second kind, Shifted Chebyshev polynomials of the third kind, Shifted Chebyshev polynomials of the fourth kind, and Gegenbauer polynomials, respectively, with different values of fractional order.

**Figure 6.** The behavior of solutions using S1 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 1.

**Figure 7.** The behavior of solutions using S2 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 1.

**Figure 8.** The behavior of solutions using S3 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 1.

**Figure 9.** The behavior of solutions using S4 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 1.

**Figure 10.** The behavior of solutions using S5 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 1.

In Table 1, we have listed the maximum absolute errors (MAE) and root-mean-square errors (RMSE) for Example 1 for the two different *n* values of 2 and 6.


**Table 1.** Result comparison of Example 1 for different orthogonal polynomials at different values of *n*.

In Table 1, we have compared results for different polynomials, and it is observed that the results for Shifted Legendre polynomials and Gegenbauer polynomials are better than those for the other polynomials. It is also observed that the MAE and RMSE decrease with increasing *n*.

**Example 2.** *Consider a non-linear fractional variational problem as in Equation (1) with g*(*x*) = *h*(*x*) = *e*−*vx; we then have the following non-linear fractional variational problem [19]:*

$$J(\mathbf{y}) = \int\_0^1 \left( e^{-\upsilon \mathbf{x}} D^\mathbf{z} y(\mathbf{x}) - \upsilon \left( I^{1-a} y(\mathbf{x}) + 1 \right) e^{-\upsilon \mathbf{x}} \right)^2 d\mathbf{x} \tag{57}$$

*under the constraints*

*y*(0) = 0, *I* <sup>1</sup>−*αy*(1) = .

*The exact solution of the above equation is given as*

$$y\_{exact}(\mathbf{x}) = \left(\boldsymbol{\varepsilon}^{-1}(1+\boldsymbol{\varepsilon}) - 1\right)\boldsymbol{\upsilon}^{-a}\left(\sum\_{k=0}^{\infty} \frac{(k+1)}{\Gamma(k+a+1)}(\boldsymbol{\upsilon}\mathbf{x})^{k+a}\right) + \boldsymbol{\upsilon}^{a-1}E\_{1,a}(\boldsymbol{\upsilon}\mathbf{x}) - \frac{\mathbf{x}^{a-1}}{\Gamma(a)}$$

*where Ea*,*b*(*x*) *is the Mittag-Leffler function of order a and b and is defined as*

$$E\_{a,b}(\mathbf{x}) = \sum\_{k=0}^{\infty} \frac{\mathbf{x}^k}{\Gamma(ak+b)}.$$

We discuss Example 2 for different *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1 and = 2.

In Figures 11–15, it is shown that the solutions for the two different values of *α* = 0.8 and *α* = 1 coincide with the exact solutions for different orthogonal polynomials at *n* = 5.

**Figure 11.** Comparison of exact and numerical solutions using S1 for *α* = 0.8 and *α* = 1, Example 2.

**Figure 12.** Comparison of exact and numerical solutions using S2 for *α* = 0.8 and *α* = 1, Example 2.

**Figure 13.** Comparison of exact and numerical solutions using S3 for *α* = 0.8 and *α* = 1, Example 2.

**Figure 14.** Comparison of exact and numerical solutions using S4 for *α* = 0.8 and *α* = 1, Example 2.

**Figure 15.** Comparison of exact and numerical solutions using S5 for *α* = 0.8 and *α* = 1, Example 2.

Figures 16–20 reflect that the approximate solution varies continuously for Shifted Legendre polynomials, Shifted Chebyshev polynomials of the second kind, Shifted Chebyshev polynomials of the third kind, Shifted Chebyshev polynomials of the fourth kind, and Gegenbauer polynomials, respectively, with different values of fractional order.

**Figure 16.** The behavior of solutions using S1 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 2.

**Figure 17.** The behavior of solutions using S2 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 2.

**Figure 18.** The behavior of solutions using S3 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 2.

**Figure 19.** The behavior of solutions using S4 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 2.

**Figure 20.** The behavior of solutions using S5 for *α* values of 0.5, 0.6, 0.7, 0.8, 0.9, and 1, Example 2.

In Table 2, we have listed the maximum absolute errors (MAE) and root-mean-square errors (RMSE) for Example 2 for the two *n* values 2 and 6.


**Table 2.** Result comparison of Example 2 for different orthogonal polynomials at different values of *n*.

In Table 2, we have compared results for different polynomials, and it is observed that the results for the Shifted Legendre polynomial are better than those for the other polynomials. It is also observed that the MAE and RMSE decrease as *n* increases.

#### **8. Conclusions**

We extended the Ritz method [18,20–22,32,33] for solving a class of NLFVPs using different orthogonal polynomials such as shifted Legendre polynomials, shifted Chebyshev polynomials of the first kind, shifted Chebyshev polynomials of the third kind, shifted Chebyshev polynomials of the fourth kind, and Gegenbauer polynomials. These polynomials were used to approximate the unknown function in the NLFVP. The advantage of the method is that it converts the given NLFVPs into a set of non-linear algebraic equations which are then solved numerically. The error bound of the approximation method for NLFVP was established. It was also shown that the approximate numerical solution converges to the exact solution as we increase the number of basis functions in the approximation. At the end, numerical results were provided by applying the method to two test examples, and it was observed that the results showed good agreement with the exact solution. Numerical results obtained using different orthogonal polynomials were compared. A comparative study showed that the shifted Legendre polynomials were more accurate in approximating the numerical solution.

**Author Contributions:** All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are very grateful to the referees for their constructive comments and suggestions for the improvement of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A**

**Theorem A1.** *Let <sup>f</sup>* : [0, 1] <sup>→</sup> *<sup>R</sup> be a function such that <sup>f</sup>* <sup>∈</sup> *<sup>C</sup>*(*N*+1)[0, 1] *and let fN*(*x*) *be the Nth approximation of the function from P*(*a*, *<sup>b</sup>*) *<sup>N</sup>* (*x*) = *span* {Ψ0(*x*), Ψ1(*x*),..., Ψ*N*(*x*)}*; then [45]*

$$\left\|f(\mathbf{x}) - f\_N(\mathbf{x})\right\|\_{\mathbf{w}^{(a,b)}}^2 \le \frac{K}{(N+1)!} \sqrt{\frac{\Gamma(1+a)\Gamma(3+2N+b)}{\Gamma(4+2N+a+b)}}\lambda$$

*where K* = *max* 9:;< *x*∈[0,1] *<sup>f</sup>* (*N*+1)(*x*) .

**Proof .** Since *<sup>f</sup>* <sup>∈</sup> *<sup>C</sup>*(*N*+1)[0, 1], the Taylor polynomial of *<sup>f</sup>* at *<sup>x</sup>* <sup>=</sup> 0, is given as

$$g\_1(\mathbf{x}) = f(0) + f'(0)\mathbf{x} + \dots + f^N(0)\frac{\mathbf{x}^N}{N!}.$$

The upper bound of the error of the Taylor polynomial is given as

$$|f(x) - g\_1(x)| \le \frac{Kx^{N+1}}{(N+1)!}r$$

where *K* = *max* 9:;< *x*∈[0,1] *<sup>f</sup>* (*N*+1)(*x*) .

Since *fN*(*x*) and *<sup>g</sup>*1(*x*) <sup>∈</sup> *<sup>P</sup>*(*a*, *<sup>b</sup>*) *<sup>N</sup>* (*x*), we have

$$\begin{split} \left\| f(\mathbf{x}) - f\_{\mathcal{N}}(\mathbf{x}) \right\|\_{\mathfrak{w}^{(a,b)}}^2 &\leq \left\| f(\mathbf{x}) - g\_{\mathcal{I}}(\mathbf{x}) \right\|\_{\mathfrak{w}^{(a,b)}}^2 \leq \left( \frac{\mathsf{K}}{(\mathsf{N}+1)!} \right)^2 \int\_0^1 \mathbf{x}^{2\mathsf{N}+2+b} (1-\mathsf{x})^a d\mathsf{x} \\ &= \left( \frac{\mathsf{K}}{(\mathsf{N}+1)!} \right)^2 \frac{\Gamma(1+a)\Gamma(3+2\mathsf{N}+b)}{\Gamma(4+2\mathsf{N}+a+b)}, \\ \left\| f(\mathbf{x}) - f\_{\mathcal{N}}(\mathbf{x}) \right\|\_{\mathfrak{w}^{(a,b)}}^2 &\leq \frac{\mathsf{K}}{(\mathsf{N}+1)!} \sqrt{\frac{\Gamma(1+a)\Gamma(3+2\mathsf{N}+b)}{\Gamma(4+2\mathsf{N}+a+b)}}, \end{split}$$

which shows that lim*N*→<sup>∞</sup> *<sup>f</sup>*(*x*) <sup>−</sup> *fN*(*x*) <sup>2</sup> *<sup>w</sup>*(*a*, *<sup>b</sup>*) = 0.
