Salih Yalçinbaş
Department of Mathematics, Manisa Celal Bayar University, Manisa, Türkiye.
DOI: 10.4236/jamp.2026.1410198   PDF    HTML   XML   0 Downloads   25 Views   Citations

Abstract

This study introduces an adaptive sparse Fibonacci collocation method for the numerical solution of ordinary differential equations with quadratic nonlinearities. Unlike conventional Fibonacci polynomial approaches based on a predefined truncation order and a full set of basis functions, the proposed framework incorporates the selection of significant Fibonacci modes directly into the solution process. The approximate solution is represented by a reduced Fibonacci expansion whose active basis is iteratively determined according to coefficient significance and residual information. Derivative terms are expressed through an operational matrix representation, while quadratic nonlinearities involving the solution and its derivatives are transformed into a nonlinear algebraic system at the collocation points. A residual-controlled enrichment and pruning mechanism is employed to balance approximation accuracy and basis complexity, allowing unnecessary Fibonacci modes to be eliminated without substantially degrading the solution. The accuracy and computational behavior of the method are assessed through benchmark nonlinear differential equations with available exact solutions using maximum absolute errors, residual norms, active basis size, and convergence characteristics. The proposed formulation provides a compact and interpretable extension of conventional Fibonacci collocation and offers an adaptive computational framework for quadratically nonlinear differential equations.

Keywords

Fibonacci Polynomials, Nonlinear Differential Equations, Quadratic Nonlinearity, Adaptive Collocation, Sparse Approximation, Operational Matrix, Residual Analysis

Share and Cite:

1. Introduction

Nonlinear ordinary differential equations (ODEs) arise in many areas of applied mathematics, physics, and engineering. In particular, equations involving quadratic interactions between the solution and its derivatives are often difficult to solve analytically, which has motivated the development of various numerical and semi-analytical approaches [1]-[13]. Polynomial and spectral techniques, including Chebyshev, Bessel, and Jacobi-based collocation methods, have proved effective for transforming such problems into finite-dimensional algebraic systems [2] [7]-[9] [13]-[16].

Fibonacci polynomials have also emerged as useful approximation functions because of their simple recurrence structure and convenient algebraic properties [17]. Fibonacci collocation and operational matrix techniques have been successfully applied to integral equations, differential equations, Lane-Emden-type problems, and fractional models [18]-[23]. In particular, Fibonacci collocation has already been employed for nonlinear differential equations [24], demonstrating the potential of this polynomial family for nonlinear problems.

Polynomial matrix-collocation techniques have also been developed specifically for differential equations with quadratic nonlinearities. For example, Bicer and Sezer [1] represented the solution by a truncated Bernoulli expansion and transformed the governing equation into a nonlinear matrix system. Their numerical results demonstrate improved accuracy as the truncation parameter is increased.

However, most conventional polynomial collocation methods employ a fixed approximation space,

$y_N\left(x\right)={\displaystyle \sum _{n=0}^N{a}_n{\varphi }_n}\left(x\right),$

in which all basis functions up to a prescribed order $N$ are retained. Such formulations do not explicitly determine which modes are actually necessary for representing the solution. Sparse approximation concepts [25]-[28] provide a natural mechanism for addressing this limitation.

Motivated by this observation, this study proposes an adaptive sparse Fibonacci collocation method for quadratically nonlinear differential equations. Instead of retaining the complete Fibonacci basis, the approximate solution is expressed as

$y_A\left(x\right)={\displaystyle \sum _{n\in A}{a}_n{F}_n}\left(x\right),A\subseteq \left\{0,\cdots ,N_{\max }\right\},$

where the active set $A$ is determined adaptively. Residual-controlled enrichment and coefficient-based pruning are employed to identify significant Fibonacci modes while maintaining the required approximation accuracy.

The main contribution of the proposed approach is therefore the integration of Fibonacci polynomial collocation, adaptive basis discovery, and sparse approximation within a unified framework for quadratic nonlinear ODEs. The method is evaluated using benchmark problems with known exact solutions, considering maximum absolute error, residual error, convergence behavior, and the number of active Fibonacci modes.

The remainder of the paper presents the Fibonacci approximation framework, mathematical formulation, proposed adaptive sparse method, error analysis, numerical experiments, and conclusions.

2. Fibonacci Polynomials and Sparse Approximation

Fibonacci polynomials provide a convenient basis for constructing polynomial approximations because of their recursive structure and simple algebraic representation [17]-[24]. They can be defined by

$F_0\left(x\right)=0,F_1\left(x\right)=1,$

and the recurrence relation

$F_{n+1}\left(x\right)=x{F}_n\left(x\right)+F_{n-1}\left(x\right),n\ge 1.$

The first few Fibonacci polynomials are

$F_1\left(x\right)=1,F_2\left(x\right)=x,F_3\left(x\right)=x^2+1,F_4\left(x\right)=x^3+2x,F_5\left(x\right)=x^4+3x^2+1.$

These functions have been successfully incorporated into collocation and operational-matrix techniques for differential and integral equations [18]-[24].

For a sufficiently smooth function $y\left(x\right)$ , the conventional truncated Fibonacci approximation can be written as

$y_N\left(x\right)={\displaystyle \sum _{n=1}^N{a}_{n}F}\left(x\right)=F\left(x\right)A,$

where

$F\left(x\right)=\left[F_1\left(x\right),F_2\left(x\right),\cdots ,F_{N+1}\left(x\right)\right]$

and

$A={\left[a_1,a_2,\cdots ,a_{N+1}\right]}^{\text{T}}.$

The derivatives of the approximation can be expressed through a Fibonacci operational differentiation matrix $D_F$ [19] [20]:

$y_N^{\left(k\right)}\left(x\right)=F\left(x\right)D_N^{\left(k\right)}A.$

This representation enables differential operators to be transformed into algebraic matrix equations and provides the basis for the collocation formulation developed in this study.

In a conventional Fibonacci approximation, all modes up to the prescribed truncation order are retained. The present approach instead considers a candidate space

$V_{N_{\max }}=\text{span}\left\{F_1,F_2,\cdots ,F_{N_{\max }+1}\right\}$

and searches for an active subset

$A\subseteq \left\{1,2,\cdots ,N_{\max }+1\right\}.$

The resulting sparse representation is

$y_A\left(x\right)={\displaystyle \sum _{n\in A}{a}_n{F}_n}\left(x\right).$

The objective is to obtain

$\left|A\right|\ll N_{\max +1}$ $\left|A\right|$

whenever the solution admits a compact Fibonacci representation. Inspired by sparse approximation principles [25]-[29], insignificant modes can be removed while new modes are introduced when they provide a meaningful reduction in the differential-equation residual.

Therefore, the approximation space is not fixed before computation but is progressively determined according to the solution structure. This adaptive representation forms the basis of the sparse Fibonacci collocation method developed in the following sections.

3. Mathematical Formulation

Consider the $m\text{th}$ -order ordinary differential equation with quadratic nonlinearities

$L\left[y\right]\left(x\right)+N\left[y\right]\left(x\right)=g\left(x\right),a\le x\le b,$ (1)

where the linear differential operator is defined by

$L\left[y\right]\left(x\right)={\displaystyle \sum _{k=0}^m{P}_k\left(x\right)y^{\left(k\right)}}\left(x\right),$ (2)

and the quadratic nonlinear operator is

$N\left[y\right]\left(x\right)={\displaystyle \sum _{p=0}^r{\displaystyle \sum _{q=0}^r{Q}_{pq}\left(x\right)y^{\left(p\right)}\left(x\right)y^{\left(q\right)}}}\left(x\right),$ (3)

Here, $P_k\left(x\right),Q_{pq}\left(x\right)$ and $g\left(x\right)$ are assumed to be sufficiently smooth functions on $\left[a,b\right]$ . The formulation includes nonlinear terms such as

$y^2,y{y}^{\prime },{\left(y^{\prime }\right)}^2,y{y}^{″},y^{\prime }{y}^{″},$

depending on the values of $p$ and $q$ .

The problem is supplemented by $m$ initial or boundary conditions, which may be written in the general form

${\displaystyle \sum _{k=0}^{m-1}\left[{\alpha }_{jk}{y}^{\left(k\right)}\left(a\right)+{\beta }_{jk}{y}^{\left(k\right)}\left(b\right)+{\gamma }_{jk}{y}^{\left(k\right)}\left(c_j\right)\right]}={\lambda }_j,j=0,1,2,\cdots ,m-1.$ (4)

This general structure permits initial, boundary, and mixed conditions of the type considered in nonlinear collocation formulations [1].

The approximate solution is sought in the Fibonacci polynomial space introduced in Section 2:

$y_A\left(x\right)={\displaystyle \sum _{n\in A}{a}_n{F}_n\left(x\right),}$ (5)

where

$A\subseteq \left\{1,2,\cdots ,N_{\max }+1\right\}$

denotes the active Fibonacci index set. In matrix form,

$y_A\left(x\right)=F_A\left(x\right)A_A,$ (6)

where

$F_A\left(x\right)=\left[F_{n_1}\left(x\right),F_{n_2}\left(x\right),\cdots ,F_{n_M}\left(x\right)\right]$

and

$A_A={\left[a_{n_1},a_{n_2},\cdots ,a_{n_M}\right]}^{\text{T}},M=\left|A\right|.$

Using the Fibonacci differentiation operator, the $k\text{th}$ derivative is represented as

$y_A^{\left(k\right)}\left(x\right)=F_A^{\left(k\right)}\left(x\right)A_A,$ (7)

Substitution of Equations (5) - (7) into Equation (1) gives the residual

$R\left(x;A_A\right)={\displaystyle \sum _{k=0}^m{P}_k\left(x\right)y_A^{\left(k\right)}\left(x\right)}+{\displaystyle \sum _{p=0}^r{\displaystyle \sum _{q=0}^r{Q}_{pq}\left(x\right)y_A^{\left(p\right)}\left(x\right)y_A^{\left(q\right)}\left(x\right)}}-g\left(x\right).$ (8)

For a set of collocation points

$x_i\in \left[a,b\right],i=0,1,\cdots ,N_c,$

the governing equation is enforced through

$R\left(x;A_A\right)=0$ (9)

Collocation Nodes and Algebraic System

For each adaptive approximation, the collocation points are selected on the computational interval using a uniform nodal rule. Let the computational interval be $t\in \left[a,b\right]$ , and let $N_c$ denote the number of collocation points. The nodes are defined by

$t_i=a+\frac{i-1}{N_c-1}\left(b-a\right),i=1,2,\cdots ,N_c.$ (10)

At every adaptive iteration, the governing differential equation is enforced at the same set of $N_c$ collocation points. Thus, for an $m\text{-th}$ order differential equation, the differential-equation residual equations are evaluated at $N_c$ points, while the prescribed initial or boundary conditions are imposed separately.

If $N_A=\left|\mathcal{A}\right|$ denotes the number of active Fibonacci modes, the resulting nonlinear system is assembled from the collocation residual equations and the required initial/boundary constraint equations. In particular, for an $m\text{-th}$ order problem with $m$ independent conditions, the system contains $N_c+m$ equations before any equation replacement or row selection is applied.

To obtain a square nonlinear algebraic system when required, $m$ collocation equations are replaced by the corresponding initial or boundary conditions. The replacement is performed at the rows associated with the selected constraint locations. Consequently, the final system contains $N_c$ equations for $N_A$ unknown Fibonacci coefficients when $N_A=N_c$ .

For first-order initial-value problems, one collocation equation is therefore replaced by the initial condition. For higher-order initial-value or boundary-value problems, the number of replaced equations is equal to the number of independent prescribed conditions.

The same node construction and equation-assembly rule is retained throughout the enrichment and pruning stages so that changes in the numerical solution can be attributed to the active Fibonacci basis rather than to a changing collocation grid.

4. Adaptive Sparse Fibonacci Collocation Method

The proposed method combines Fibonacci polynomial collocation with adaptive basis selection. Unlike conventional Fibonacci and other polynomial collocation approaches using a fixed truncation order [18] [22] [30], only significant Fibonacci modes are retained during the approximation process.

Let

$F\left(x\right)=\left[F_1\left(x\right),F_2\left(x\right),\cdots ,F_{N_{\max }+1}\left(x\right)\right]$

be the candidate Fibonacci basis. For an active index set $A$ , the approximate solution is written as

$y_A\left(x\right)={\displaystyle \sum _{j\in A}{a}_j}{F}_j\left(x\right)=F_A\left(x\right)A_A.$ (11)

Using the Fibonacci operational differentiation matrix $D_F$ [19] [20], the derivatives are represented by

$y_A^{\left(k\right)}\left(x\right)=F_A^{\left(k\right)}\left(x\right)A_A.$ (12)

4.1. Differentiation Operator for a Nonconsecutive Active Fibonacci Set

Because the adaptive procedure may produce a nonconsecutive active Fibonacci index set, the differentiation operator cannot in general be defined by simply restricting the full Fibonacci differentiation matrix to the active indices. Let the candidate basis be

$\Phi \left(t\right)={\left[F_0\left(t\right),F_1\left(t\right),\cdots ,F_M\left(t\right)\right]}^{\text{T}},$

and let the corresponding full differentiation matrix be $D$ , defined by

$\frac{\text{d}}{\text{d}t}\Phi \left(t\right)=D\Phi \left(t\right).$ (13)

For an active index set

$\mathcal{A}=\left\{j_1,j_2,\cdots ,j_{N_A}\right\},j_1<j_2<\cdots <j_{N_A},$

the approximation is written as

$y_N\left(t\right)={\Phi }_{\mathcal{A}}^{\text{T}}\left(t\right)c_{\mathcal{A}},$

where ${\Phi }_{\mathcal{A}}$ contains only the active Fibonacci modes. The derivative is nevertheless evaluated through the full candidate-space differentiation operator:

${y^{\prime }}_N\left(t\right)=D_{\mathcal{C},\mathcal{A}}{c}_{\mathcal{A}}\cdot \Phi \left(t\right),$ (14)

where $D_{\mathcal{C},\mathcal{A}}$ denotes the columns of the full differentiation matrix corresponding to the active modes.

Thus, if the derivative of an active mode $F_{jk}\left(t\right)$ contains Fibonacci modes whose indices are not members of $\mathcal{A}$ , those modes are not discarded. They remain in the derivative representation within the candidate space and are included when the differential-equation residual is evaluated.

More explicitly,

${F^{\prime }}_{jk}\left(t\right)={\displaystyle \sum _{\ell =0}^M{D}_{\ell j_k}}{F}_{\ell }\left(t\right),j_k\in \mathcal{A}.$ (15)

Consequently,

${y^{\prime }}_N\left(t\right)={\displaystyle \sum _{j\in \mathcal{A}}{c}_j}{\displaystyle \sum _{\ell =0}^M{D}_{\ell j}}{F}_{\ell }\left(t\right).$ (16)

The same construction is applied recursively for higher-order derivatives. For example,

$y_N^{\left(r\right)}\left(t\right)={\Phi }^{\text{T}}\left(t\right)D_{\mathcal{C},\mathcal{A}}^r{c}_{\mathcal{A}},$ (17)

where the active columns are selected only after the full candidate-space differentiation operator has been constructed.

This formulation separates basis selection from differentiation. The adaptive procedure determines which Fibonacci modes are used in the approximation itself, whereas differentiation is performed in the full candidate space. Therefore, nonconsecutive active sets do not introduce an artificial truncation of derivative contributions. A mode is removed from the approximation only through the prescribed pruning criterion, and not merely because its derivative representation is outside the active set.

4.2. Nonlinear Solution Procedure

At each adaptive iteration, the nonlinear algebraic system associated with the current active Fibonacci set is solved for the unknown expansion coefficients using the MATLAB fsolve solver. For the initial active set, the coefficient vector is initialized consistently with the prescribed initial or boundary conditions. At subsequent stages, the converged coefficient vector obtained from the preceding iteration serves as an effective initial guess for the current solve. The iterative solution process is terminated when both the nonlinear-system residual and the coefficient update satisfy a numerical tolerance of ${\epsilon }_{\text{solver}}={10}^{-10}$ .

The adaptive residual tolerance and the coefficient-pruning tolerance are denoted by ${\tau }_{\text{res}}$ and ${\tau }_{\text{prune}}$ , respectively. Following a successful nonlinear solve, the differential residual norm is evaluated to govern the adaptive enrichment and pruning criteria. If a proposed basis modification leads to a failure in fsolve convergence within the prescribed iteration limit, the tentative modification is rejected, and the previous converged active basis and coefficient vector are strictly retained. The procedure then either investigates the next admissible basis adjustment or terminates if convergence criteria are met, ensuring a consistent and robust numerical implementation across all benchmark problems.

Algorithm 1. Adaptive Sparse Fibonacci Collocation Method

Input: Candidate Fibonacci index set $\mathcal{C}$ , initial active set ${\mathcal{A}}^{\left(0\right)}$ , collocation nodes, initial/boundary conditions, nonlinear-solver tolerance ${\epsilon }_{NL}$ , residual tolerance ${\epsilon }_R$ , coefficient-pruning tolerance ${\epsilon }_c$ , and maximum number of adaptive iterations $K_{\max }$ .

Output: Final active index set ${\mathcal{A}}^{\ast }$ and coefficient vector $c^{\ast }$ .

1. Initialize $\mathcal{A}\leftarrow {\mathcal{A}}^{\left(0\right)}$ and construct the corresponding Fibonacci approximation.

2. Assemble the nonlinear collocation system using the current active set and the prescribed initial/boundary-conditions.

3. Initialize the coefficient vector. For the first iteration, use the prescribed conditions to determine the constrained coefficients and initialize the remaining coefficients. At subsequent iterations, use the previously converged coefficient vector as the initial guess.

4. Solve the nonlinear algebraic system using the selected nonlinear solver until the prescribed nonlinear-solver tolerance ${\epsilon }_{NL}$ is satisfied.

5. If the nonlinear solve is unsuccessful, reject the current tentative basis modification and retain the previous converged active set and coefficient vector.

6. Evaluate the collocation residual

$ℛ\left(\mathcal{A}\right)=\underset{i}{\max }\left|R\left(t_i\right)\right|.$

7. Apply the coefficient-based pruning criterion and remove a mode only when its removal satisfies the prescribed pruning condition.

8. For each inactive candidate mode, evaluate the residual reduction obtained by tentative enrichment.

9. If several candidate modes satisfy the enrichment criterion, select the mode giving the largest admissible residual reduction.

10. Add the selected mode, reassemble the nonlinear system, and re-solve for all active coefficients.

11. If the nonlinear solve associated with the enrichment is unsuccessful, reject that enrichment and retain the previous converged representation.

12. Repeat the enrichment-pruning cycle until

$ℛ\left(\mathcal{A}\right)\le {\epsilon }_R$ ,

no admissible enrichment produces a meaningful residual reduction, and no active coefficient satisfies the pruning criterion.

13. If $K_{\max }$ is reached, terminate the adaptive procedure and retain the latest converged representation.

14. Return ${\mathcal{A}}^{\ast }$ and $c^{\ast }$ .

5. Convergence and Error Analysis

The accuracy of the proposed adaptive sparse Fibonacci approximation is evaluated using the residual and, whenever available, the exact solution. For the active basis $A$ , the residual is defined as

$R_A\left(x\right)=L\left[y_A\right]\left(x\right)+N\left[y_A\right]\left(x\right)-g\left(x\right).$

The corresponding residual norm is

$E_R={\Vert R_A\Vert }_{\infty }=\underset{x\in \left[a,b\right]}{\max }\left|R_A\left(x\right)\right|.$

When the exact solution $y\left(x\right)$ is known, the maximum absolute error is computed by

$E_{\infty }={\Vert y-y_A\Vert }_{\infty }=\underset{x\in \left[a,b\right]}{\max }\left|y\left(x\right)-y_A\left(x\right)\right|.$

The adaptive enrichment process generates a sequence of active spaces

$V_{A_1}\subseteq V_{A_2}\subseteq \cdots ,$

where new Fibonacci modes are retained only when they provide a sufficient reduction in the residual. Thus, the computational convergence criterion is expressed as

$E_R^{\left(s+1\right)}<E_R^{\left(s\right)},$

where $s$ denotes the adaptive iteration.

For sufficiently smooth solutions that admit a convergent Fibonacci expansion, increasing the effective approximation space is expected to reduce the truncation error. The pruning mechanism preserves this behavior by removing a mode only when its elimination produces no significant deterioration in the residual.

To quantify the efficiency of the sparse representation, the compression ratio is defined as

$C_R=1-\frac{\left|A^{*}\right|}{N_{\max }+1},$

where $\left|A^{*}\right|$ is the number of active Fibonacci modes in the final approximation. A larger $C_R$ indicates a more compact representation.

Accordingly, the proposed method is assessed through three complementary quantities:

$\boxed{E_{\infty },E_R,C_R}$

representing solution accuracy, equation consistency, and basis sparsity, respectively. These measures are used in the numerical experiments to evaluate the convergence and computational effectiveness of the adaptive Fibonacci approximation.

6. Numerical Experiments

The proposed Adaptive Sparse Fibonacci Collocation (ASFC) method is evaluated on three nonlinear differential equations with known exact solutions. The benchmark problems include a nonlinear oscillator, a quadratic Riccati equation, and a third-order nonlinear equation. These examples are selected to examine the accuracy of the proposed method for oscillatory, exponential, and strongly nonlinear solution behaviors. In addition to the solution error, the differential-equation residual and the sparsity of the final Fibonacci representation are considered.

6.1. Problem 1—Nonlinear Oscillator

Consider the nonlinear differential equation

$u^{″}\left(t\right)-u\left(t\right)+u^2\left(t\right)+{\left[u^{\prime }\left(t\right)\right]}^2-1=0,0\le t\le 1,$

subject to the initial conditions

$u\left(0\right)=2,u^{\prime }\left(0\right)=0.$

The exact solution is

$u_{\text{ex}}\left(t\right)=1+\cos t.$

For an active Fibonacci index set $\mathcal{A}$ , the ASFC approximation is written as

$u_N\left(t\right)={\displaystyle \sum _{j\in \mathcal{A}}{c}_j}{F}_j\left(t\right).$

The initial conditions are imposed directly on the approximation,

${\displaystyle \sum _{j\in \mathcal{A}}{c}_j}{F}_j\left(0\right)=2,{\displaystyle \sum _{j\in \mathcal{A}}{c}_j}{F^{\prime }}_j\left(0\right)=0.$

The remaining coefficients are determined from the nonlinear collocation residual

$R_{1,N}\left(t\right)={u^{″}}_N\left(t\right)-u_N\left(t\right)+u_N^2\left(t\right)+{\left[{u^{\prime }}_N\left(t\right)\right]}^2-1,$

evaluated at the prescribed collocation points. The adaptive procedure introduces additional Fibonacci modes when the residual reduction obtained from an inactive candidate mode satisfies the enrichment criterion. Conversely, insignificant active modes are removed when their elimination does not produce a significant deterioration of the residual.

This problem is particularly useful for examining sparse basis selection because the exact solution contains a constant component and an oscillatory component while the governing equation includes quadratic nonlinearities.

6.2. Problem 2—Quadratic Riccati Equation

The second benchmark is a quadratic Riccati equation defined on

$0\le t\le 1$

by

$u^{\prime }\left(t\right)={\text{e}}^t-{\text{e}}^{3t}+2{\text{e}}^{2t}u\left(t\right)-{\text{e}}^t{u}^2\left(t\right),$ (18)

with the initial condition

$u\left(0\right)=1,$ (19)

and exact solution

$u_{\text{ex}}\left(t\right)={\text{e}}^t.$

Equivalently, the equation can be written in the standard quadratic Riccati form

$u^{\prime }\left(t\right)=P\left(t\right)u^2\left(t\right)+Q\left(t\right)u\left(t\right)+g\left(t\right),$

where

$P\left(t\right)=-{\text{e}}^t,Q\left(t\right)=2{\text{e}}^{2t},g\left(t\right)={\text{e}}^t-{\text{e}}^{3t}.$

Using the ASFC representation, the exact nonlinear differential equation is transformed into a finite-dimensional nonlinear collocation system. The coefficient associated with the initial condition is constrained through Equation (19), while the remaining coefficients are determined from the collocation equations generated by the residual of Equation (18).

Because the exponential exact solution contains both even and odd powers in its Taylor expansion, this benchmark provides a useful test of adaptive Fibonacci basis enrichment and the ability of the proposed method to construct a compact representation without imposing a fixed complete truncation.

6.3. Problem 3—Third-Order Nonlinear Equation

The third benchmark considers the nonlinear third-order differential equation

$u^{‴}\left(t\right)=-2{\text{e}}^{-3u\left(t\right)}+\frac{4}{{\left(1+t\right)}^3},0\le t\le 1,$

subject to

$u\left(0\right)=0,u^{\prime }\left(0\right)=1,u^{″}\left(0\right)=-1.$

The exact solution is

$u_{\text{ex}}\left(t\right)=\ln \left(1+t\right).$

Starting from the active Fibonacci representation,

$u_N\left(t\right)={\displaystyle \sum _{j\in \mathcal{A}}{c}_j}{F}_j\left(t\right),$

the required derivatives are evaluated using the Fibonacci differentiation operator described in Section 4. The three initial conditions are imposed through

$u_N\left(0\right)=0,{u^{\prime }}_N\left(0\right)=1,{u^{″}}_N\left(0\right)=-1.$

The remaining coefficients are obtained from the nonlinear collocation residual

$R_{3,N}\left(t\right)={u^{‴}}_N\left(t\right)+2{\text{e}}^{-3u_N\left(t\right)}-\frac{4}{{\left(1+t\right)}^3}.$

Additional Fibonacci modes are introduced when they provide a meaningful reduction of the residual, whereas insignificant modes are removed according to the pruning criterion. The process continues until the prescribed residual criterion is satisfied or no admissible enrichment provides a sufficient improvement.

This problem also provides a test of the proposed method for a higher-order nonlinear equation and explicitly verifies the treatment of the second derivative initial condition $u^{″}\left(0\right)=-1$ .

6.4. Exact-Solution Verification

Before computing the numerical errors, the exact solutions used in the three benchmark problems are independently verified against their corresponding governing differential equations and initial conditions. The verification is performed symbolically and is independent of the ASFC approximation. For each problem, the exact-solution verification residual is defined by substituting the exact solution into the corresponding differential equation.

6.4.1. Verification of Problem 1

For Problem 1,

$u_{\text{ex}}\left(t\right)=1+\cos t,$

and therefore

${u^{\prime }}_{\text{ex}}\left(t\right)=-\sin t,{u^{″}}_{\text{ex}}\left(t\right)=-\cos t.$

Substitution into the governing equation gives

$\begin{array}{c}{R}_{1,\text{ex}}\left(t\right)={u^{″}}_{\text{ex}}\left(t\right)-u_{\text{ex}}\left(t\right)+u_{\text{ex}}^2\left(t\right)+{\left[{u^{\prime }}_{\text{ex}}\left(t\right)\right]}^2-1\\ =-\cos t-\left(1+\cos t\right)+{\left(1+\cos t\right)}^2+{\sin }^{2}t-1\\ ={\cos }^{2}t+{\sin }^{2}t-1=0.\end{array}$

The initial conditions are also satisfied exactly:

$u_{\text{ex}}\left(0\right)=2,{u^{\prime }}_{\text{ex}}\left(0\right)=0.$

Consequently,

${\eta }_{ver,1}=\underset{0\le t\le 1}{\max }\left|R_{1,\text{ex}}\left(t\right)\right|=0.$

6.4.2. Verification of Problem 2

For Problem 2,

$u_{\text{ex}}\left(t\right)={\text{e}}^t,{u^{\prime }}_{\text{ex}}\left(t\right)={\text{e}}^t.$

Substitution into the governing Riccati equation yields

$\begin{array}{c}{R}_{2,\text{ex}}\left(t\right)={u^{\prime }}_{\text{ex}}\left(t\right)-{\text{e}}^t+{\text{e}}^{3t}-2{\text{e}}^{2t}{u}_{\text{ex}}\left(t\right)+{\text{e}}^t{u}_{\text{ex}}^2\left(t\right)\\ ={\text{e}}^t-{\text{e}}^t+{\text{e}}^{3t}-2{\text{e}}^{3t}+{\text{e}}^{3t}=0.\end{array}$

The prescribed initial condition is verified by

$u_{\text{ex}}\left(0\right)={\text{e}}^0=1.$

Hence,

${\eta }_{ver,2}=\underset{0\le t\le 1}{\max }\left|R_{2,\text{ex}}\left(t\right)\right|=0.$

6.4.3. Verification of Problem 3

For Problem 3,

$u_{\text{ex}}\left(t\right)=\ln \left(1+t\right),$

with derivatives

${u^{\prime }}_{\text{ex}}\left(t\right)=\frac{1}{1+t},{u^{″}}_{\text{ex}}\left(t\right)=-\frac{1}{{\left(1+t\right)}^2},{u^{‴}}_{\text{ex}}\left(t\right)=\frac{2}{{\left(1+t\right)}^3}.$

Furthermore,

${\text{e}}^{-3u_{\text{ex}}\left(t\right)}={\text{e}}^{-3\ln \left(1+t\right)}\frac{1}{{\left(1+t\right)}^3}.$

Therefore,

$R_{3,\text{ex}}\left(t\right)={u^{‴}}_{\text{ex}}\left(t\right)+2{\text{e}}^{-3u_{\text{ex}}\left(t\right)}-\frac{4}{{\left(1+t\right)}^3}=\frac{2}{{\left(1+t\right)}^3}+\frac{2}{{\left(1+t\right)}^3}-\frac{4}{{\left(1+t\right)}^3}=0.$

The initial conditions are satisfied exactly:

$u_{\text{ex}}\left(0\right)=0,{u^{\prime }}_{\text{ex}}\left(0\right)=1,{u^{″}}_{\text{ex}}\left(0\right)=-1.$

Thus,

${\eta }_{ver,3}=\underset{0\le t\le 1}{\max }\left|R_{3,\text{ex}}\left(t\right)\right|=0.$

The symbolic verification results are summarized in Table 1.

Table 1. Symbolic verification of the exact solutions for the three benchmark problems.

Problem	Exact solution	Verification residual	Initial-condition mismatch
Problem 1	$1+\cos t$	0	0
Problem 2	${\text{e}}^t$	0	0
Problem 3	$\ln \left(1+t\right)$	0	0

These results confirm that the three exact solutions satisfy their respective governing equations and prescribed initial conditions identically. The subsequent numerical error calculations are therefore performed relative to independently verified reference solutions.

6.5. Error Measures

Let $u_N\left(t\right)$ denote the ASFC approximation and $u_{\text{ex}}\left(t\right)$ the corresponding exact solution. The pointwise absolute error is defined by

$E\left(t\right)=\left|u_{\text{ex}}\left(t\right)-u_N\left(t\right)\right|.$

The maximum absolute error over the computational interval is

$E_{\infty }=\underset{t\in \left[a,b\right]}{\max }\left|u_{\text{ex}}\left(t\right)-u_N\left(t\right)\right|.$

For numerical evaluation, the maximum is calculated on the prescribed error-evaluation grid ${\left\{t_q\right\}}_{q=1}^{N_e}$ ,

$E_{\infty }=\underset{1\le q\le N_e}{\max }\left|u_{\text{ex}}\left(t_q\right)-u_N\left(t_q\right)\right|.$

In addition to the solution error, the maximum absolute differential-equation residual is used to measure the consistency of the numerical approximation with the governing equation. For a general differential operator $ℱ$ ,

$R_N\left(t\right)=ℱu_N,$

and

$R_{\infty }=\underset{t\in \left[a,b\right]}{\max }\left|R_N\left(t\right)\right|.$

On the numerical evaluation grid, this becomes

$R_{\infty }=\underset{1\le i\le N_c}{\max }\left|R_N\left(t_i\right)\right|.$

For the three benchmark problems, the corresponding residuals are

$R_{1,N}\left(t\right)={u^{″}}_N\left(t\right)-u_N\left(t\right)+u_N^2\left(t\right)+{\left[{u^{\prime }}_N\left(t\right)\right]}^2-1,$

$R_{2,N}\left(t\right)={u^{\prime }}_N\left(t\right)-{\text{e}}^t+{\text{e}}^{3t}-2{\text{e}}^{2t}{u}_N\left(t\right)+{\text{e}}^t{u}_N^2\left(t\right),$

and

$R_{3,N}\left(t\right)={u^{‴}}_N\left(t\right)+2{\text{e}}^{-3u_N\left(t\right)}-\frac{4}{{\left(1+t\right)}^3}.$

Accordingly, the problem-dependent maximum residual is defined as

$R_{\infty ,k}\left(t\right)=\underset{1\le i\le N_c}{\max }\left|R_{k,N}\left(t_i\right)\right|,k=1,2,3.$

The sparsity of the final Fibonacci representation is quantified by the number of retained active modes,

$N_{\text{act}}=\left|\mathcal{A}\right|,$

where $\mathcal{A}$ is the final active Fibonacci index set. The corresponding compression or sparsity ratio is defined by

$S=1-\frac{N_{\text{act}}}{N_{\text{cand}}},$

where $N_{\text{cand}}$ denotes the number of candidate Fibonacci modes considered by the adaptive procedure. A larger value of $S$ corresponds to a more compact final representation.

Thus, the numerical assessment is based on three complementary measures: $E_{\infty }$ , $R_{\infty }$ and $N_{\text{act}}$ (or $S$ ). The first measures the deviation from the independently verified exact solution, the second measures satisfaction of the governing differential equation, and the third quantifies the compactness of the adaptive Fibonacci representation.

Figures 1(a)-(c) provide graphical comparisons of the ASFC results for the three benchmark problems with respect to the independent variable$t\in \left[0,1\right]$ . For each individual problem, the numerical approximation $u\left(t\right)$ , pointwise absolute error $\left|e\left(t\right)\right|$ , and differential-equation residual $\left|R\left(t\right)\right|$ are systematically displayed alongside the Bernoulli polynomial results [1]. These subfigures serve as an explicit visual complement to the quantitative error and residual measures defined above, while the detailed discussion and interpretation of the numerical performance are presented in Section 7.

(a)

(b)

(c)

Figure 1. (a) Numerical solution, pointwise absolute error, and differential-equation residual for Problem 1 (Nonlinear Oscillator) obtained by the proposed Adaptive Sparse Fibonacci Collocation (ASFC) method and the Bernoulli polynomial method [1]; (b) Numerical solution, pointwise absolute error, and differential-equation residual for Problem 2 (Riccati Equation) obtained by the proposed Adaptive Sparse Fibonacci Collocation (ASFC) method and the Bernoulli polynomial method [1]; (c) Numerical solution, pointwise absolute error, and differential-equation residual for Problem 3 (Third-Order Nonlinear Equation) obtained by the proposed Adaptive Sparse Fibonacci Collocation (ASFC) method and the Bernoullipoly nomial method [1].

7. Results and Discussion

The numerical performance of the proposed Adaptive Sparse Fibonacci Collocation (ASFC) method is assessed using the three benchmark problems introduced in Section 6. The comparison considers three complementary indicators: the maximum absolute error $E_{\infty }$ , the maximum differential-equation residual $R_{\infty }$ , the number of active Fibonacci modes $N_{\text{act}}$ , and the corresponding compression ratio. The exact solutions used for the error calculations were independently verified against the governing equations and prescribed initial conditions, as reported in Section 6.4.

7.1. Summary of Numerical Results

The main numerical characteristics of the three benchmark problems are summarized in Table 2.

Table 2. Numerical performance of the proposed ASFC method.

Problem	Solution type	Candidate modes $N_{\text{cand}}$	Active modes $N_{\text{act}}$	Compression ratio	$E_{\infty }$	$R_{\infty }$
Problem 1	Oscillatory	25	7	72.0%	$3.86\times {10}^{-9}$	$4.80\times {10}^{-8}$
Problem 2	Exponential/Riccati	25	8	68.0%	$1.01\times {10}^{-11}$	$1.02\times {10}^{-10}$
Problem 3	Logarithmic, third-order	25	10	60.0%	$9.65\times {10}^{-10}$	$8.70\times {10}^{-9}$

The results show that the adaptive procedure does not retain the complete candidate Fibonacci basis. For Problem 1, only seven modes are retained, whereas eight modes are sufficient for Problem 2 and ten modes are retained for Problem 3. This behavior illustrates the principal sparse-approximation feature of the proposed formulation.

For Problem 2, the retention of eight modes from a 25-mode candidate basis corresponds to a compression ratio of 68.0%. Thus, approximately one-third of the candidate modes are retained in the final representation. The result is particularly relevant because the exact solution is exponential rather than polynomial, providing a useful test of the ability of the Fibonacci basis to represent smooth non-polynomial behavior.

7.2. Solution Accuracy

Figures 1(a)-(c) compares the numerical and exact solutions for all three benchmark problems together with the corresponding pointwise absolute errors and differential-equation residuals.

For Problem 1, the computed solution reproduces the oscillatory behavior of the exact solution over the complete interval. The use of seven active Fibonacci modes demonstrates that the oscillatory solution can be represented without retaining the entire candidate basis.

Problem 2 provides a different test because its exact solution is $u\left(t\right)={\text{e}}^t$ . Despite the exponential character of the solution and the quadratic Riccati nonlinearity, the adaptive procedure retains only eight modes. The resulting approximation follows the exact solution throughout the computational interval, while the pointwise error remains controlled.

Problem 3 represents the most demanding benchmark in terms of differential order, since the governing equation contains a third derivative and a nonlinear exponential term. The exact solution $u\left(t\right)=\ln \left(1+t\right)$ is also non-polynomial. The final approximation uses ten active Fibonacci modes and reproduces the qualitative behavior of the reference solution over $0\le t\le 1$ .

The maximum absolute errors reported in Table 2 provide the quantitative counterpart of these observations. In contrast to a visual comparison alone, $E_{\infty }$ measures the largest deviation between the numerical and independently verified exact solutions over the prescribed computational grid.

7.3. Residual Accuracy

Solution error alone is not sufficient to characterize the performance of a nonlinear collocation method. Therefore, the differential-equation residual is evaluated independently using

$R_N\left(t\right)=ℱu_N.$

The corresponding maximum residual $R_{\infty }$ measures the extent to which the computed approximation satisfies the governing differential equation.

The residual curves depicted in Figure 1(a3), Figure 1(b3), and Figure 1(c3) remain strictly bounded across the computational domain $t\in \left[0,1\right]$ for all three benchmark problems. In each case, the residual magnitude $\left|R\left(t\right)\right|$ associated with the ASFC method remains consistently suppressed below 10⁻⁷, demonstrating substantial consistency with the governing differential equations compared to the Bernoulli polynomial approach [1]. This behavior is directly aligned with the residual-controlled adaptive strategy introduced in Section 4, wherein inactive Fibonacci modes are selectively incorporated into the active basis only when they provide an admissible reduction in the differential residual.

The residual results are particularly important for Problem 3. Because the equation is third-order and contains the stiff nonlinear exponential term$-2{\text{e}}^{-3u\left(t\right)}$ , subsequent differentiation and nonlinear evaluations may significantly amplify numerical approximation errors. As demonstrated in Figure 1(c3), the residual $\left|R\left(t\right)\right|$ remains bounded on the order of 10⁻⁸, thereby providing an independent confirmation that the obtained sparse Fibonacci approximation rigorously satisfies the original governing equation rather than merely reproducing the qualitative shape of the exact solution.

7.4. Adaptive Basis Selection and Sparsity

The numerical results also illustrate the role of adaptive basis selection. The final approximation is written as

$u_N\left(t\right)={\displaystyle \sum _{j\in \mathcal{A}}{c}_j}{F}_j\left(t\right),$

where $\mathcal{A}$ is the adaptively determined active index set.

The number of retained modes is seven, eight, and ten for Problems 1-3, respectively. These values indicate that the required approximation dimension depends on the structure and differential order of the problem rather than being fixed uniformly in advance.

The compression ratio provides an additional measure of this behavior. In particular, the 68.0% compression obtained for Problem 2 demonstrates that the candidate space can be substantially reduced while maintaining the required numerical accuracy. The same interpretation applies to Problems 1 and 3, for which the same 25-mode candidate space is reduced to 7 and 10 active modes, corresponding to compression ratios of 72.0% and 60.0% respectively (Table 2).

An important consequence is that increasing the candidate basis does not automatically imply that all additional modes must be retained. Instead, the adaptive procedure attempts to identify modes that contribute meaningfully to residual reduction. This distinguishes the proposed ASFC formulation from conventional fixed-order Fibonacci collocation approaches in which the complete truncated basis is retained.

7.5. Comparison with the Bernoulli Polynomial Method

Figures 1(a)-(c) also provides a comparison with the Bernoulli polynomial results reported by Bicer and Sezer [1]. The comparison includes the numerical solution, pointwise absolute error, and differential-equation residual for the benchmark problems.

The comparison should be interpreted using the reported numerical quantities rather than the solution curves alone. In particular, $E_{\infty }$ characterizes approximation accuracy relative to the exact solution, while $R_{\infty }$ measures consistency with the governing differential equation. The number of active Fibonacci modes additionally characterizes the compactness of the proposed representation.

The ASFC formulation differs from the Bernoulli approach primarily in the construction of the approximation space. The Bernoulli formulation uses a prescribed polynomial representation, whereas the present method introduces a residual-controlled mechanism for selecting an active subset of Fibonacci modes. Consequently, the comparison provides information not only about approximation accuracy but also about the effect of adaptive basis selection.

7.6. Overall Numerical Assessment

Taken together, the three experiments cover oscillatory, exponential, and logarithmic solution behaviors and include first- and third-order nonlinear differential operators. The numerical results indicate that the proposed ASFC framework can represent these solutions using relatively small active Fibonacci spaces.

Three observations are particularly relevant. First, the numerical solutions remain consistent with the independently verified exact solutions. Second, the differential-equation residual provides an additional measure of numerical consistency and forms the basis of the adaptive enrichment mechanism. Third, the final active basis contains fewer modes than the corresponding candidate space, demonstrating the sparse character of the proposed approximation.

The results therefore support the central computational concept of the method: rather than prescribing a fixed approximation order and retaining every available basis function, the ASFC procedure seeks an active Fibonacci representation adapted to the nonlinear residual structure of the problem.

A complete assessment of computational efficiency would additionally require reporting implementation-dependent quantities such as CPU time, nonlinear-solver iterations, conditioning measures, and the number of enrichment/pruning iterations. These quantities are not introduced here unless they are directly available from the numerical implementation.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Bicer, K. and Sezer, M. (2019) A Computational Method for Solving Differential Equations with Quadratic Nonlinearity by Using Bernoulli Polynomials. Thermal Science, 23, 275-283.[CrossRef]
[2]	Chen, B., García-Bolós, R., Jódar, L. and Roselló, M.D. (2005) Chebyshev Polynomial Approximations for Nonlinear Differential Initial Value Problems. Nonlinear Analysis: Theory, Methods & Applications, 63, e629-e637.[CrossRef]
[3]	Abbasbandy, S. (2007) A New Application of He’s Variational Iteration Method for Quadratic Riccati Differential Equation by Using Adomian’s Polynomials. Journal of Computational and Applied Mathematics, 207, 59-63.[CrossRef]
[4]	Lakestani, M. and Dehghan, M. (2010) Numerical Solution of Riccati Equation Using the Cubic B-Spline Scaling Functions and Chebyshev Cardinal Functions. Computer Physics Communications, 181, 957-966.[CrossRef]
[5]	Abbasbandy, S. (2006) Homotopy Perturbation Method for Quadratic Riccati Differential Equation and Comparison with Adomian’s Decomposition Method. Applied Mathematics and Computation, 172, 485-490.[CrossRef]
[6]	El-Tawil, M.A., Bahnasawi, A.A. and Abdel-Naby, A. (2004) Solving Riccati Differential Equation Using Adomian’s Decomposition Method. Applied Mathematics and Computation, 157, 503-514.[CrossRef]
[7]	Akyüz-Daşcıoğlu, A. and Çerdi˙k-Yaslan, H. (2011) The Solution of High-Order Nonlinear Ordinary Differential Equations by Chebyshev Series. Applied Mathematics and Computation, 217, 5658-5666.[CrossRef]
[8]	Yüzbaşı, Ş. (2012) A Numerical Approximation Based on the Bessel Functions of First Kind for Solutions of Riccati Type Differential-Difference Equations. Computers & Mathematics with Applications, 64, 1691-1705.[CrossRef]
[9]	Yüzbaşı, Ş. (2017) A Numerical Scheme for Solutions of a Class of Nonlinear Differential Equations. Journal of Taibah University for Science, 11, 1165-1181.[CrossRef]
[10]	Ghomanjani, F. and Khorram, E. (2017) Approximate Solution for Quadratic Riccati Differential Equation. Journal of Taibah University for Science, 11, 246-250.[CrossRef]
[11]	Behiry, S.H., Hashish, H., El-Kalla, I.L. and Elsaid, A. (2007) A New Algorithm for the Decomposition Solution of Nonlinear Differential Equations. Computers & Mathematics with Applications, 54, 459-466.[CrossRef]
[12]	Barari, A., Omidvar, M., Ghotbi, A.R. and Ganji, D.D. (2008) Application of Homotopy Perturbation Method and Variational Iteration Method to Nonlinear Oscillator Differential Equations. Acta Applicandae Mathematicae, 104, 161-171.[CrossRef]
[13]	Bhrawy, A.H. and Abd-Elhameed, W.M. (2011) New Algorithm for the Numerical Solutions of Nonlinear Third‐Order Differential Equations Using Jacobi‐Gauss Collocation Method. Mathematical Problems in Engineering, 2011, Article ID: 837218.[CrossRef]
[14]	Doha, E.H., Abd-Elhameed, W.M. and Youssri, Y.H. (2013) Second Kind Chebyshev Operational Matrix Algorithm for Solving Differential Equations of Lane-Emden Type. New Astronomy, 23, 113-117.[CrossRef]
[15]	Fornberg, B. (1998) A Practical Guide to Pseudospectral Methods. Cambridge University Press.
[16]	Trefethen, L.N. (2000) Spectral Methods in MATLAB. Society for Industrial and Applied Mathematics.[CrossRef]
[17]	Falcón, S. and Plaza, Á. (2009) On K-Fibonacci Sequences and Polynomials and Their Derivatives. Chaos, Solitons & Fractals, 39, 1005-1019.[CrossRef]
[18]	Mirzaee, F. and Hoseini, S.F. (2016) Application of Fibonacci Collocation Method for Solving Volterra-Fredholm Integral Equations. Applied Mathematics and Computation, 273, 637-644.[CrossRef]
[19]	Abd-Elhameed, W. and Youssri, Y. (2016) A Novel Operational Matrix of Caputo Fractional Derivatives of Fibonacci Polynomials: Spectral Solutions of Fractional Differential Equations. Entropy, 18, Article 345.[CrossRef]
[20]	Çakmak, M. (2019) Fibonacci Operational Matrix Algorithm for Solving Differential Equations of Lane-Emden Type. Sakarya University Journal of Science, 23, 478-485.[CrossRef]
[21]	Ali, I., Haq, S., Nisar, K.S. and Arifeen, S.U. (2021) Numerical Study of 1D and 2D Advection-Diffusion-Reaction Equations Using Lucas and Fibonacci Polynomials. Arabian Journal of Mathematics, 10, 513-526.[CrossRef]
[22]	Postavaru, O. (2023) An Efficient Numerical Method Based on Fibonacci Polynomials to Solve Fractional Differential Equations. Mathematics and Computers in Simulation, 212, 406-422.[CrossRef]
[23]	S M, S., Kumar, P., Govindaraj, V., Qahiti, R.A., Hamali, W. and Mutum, Z.M. (2024) An Operational Matrix Approach with Vieta-Fibonacci Polynomial for Solving Generalized Caputo Fractal-Fractional Differential Equations. Ain Shams Engineering Journal, 15, Article 102678.[CrossRef]
[24]	Çakmak, M. and Alkan, S. (2022) Fibonacci Collocation Method for Solving a Class of Nonlinear Differential Equations. Turkish Journal of Mathematics and Computer Science, 15, 42-53.[CrossRef]
[25]	Canuto, C., Hussaini, M.Y., Quarteroni, A. and Zang, T.A. (2006) Spectral Methods: Fundamentals in Single Domains. Springer.
[26]	Tibshirani, R. (1996) Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58, 267-288.[CrossRef]
[27]	Donoho, D.L. (2006) Compressed Sensing. IEEE Transactions on Information Theory, 52, 1289-1306.[CrossRef]
[28]	Tropp, J.A. (2004) Greed Is Good: Algorithmic Results for Sparse Approximation. IEEE Transactions on Information Theory, 50, 2231-2242.[CrossRef]
[29]	Bungartz, H. and Griebel, M. (2004) Sparse Grids. Acta Numerica, 13, 147-269.[CrossRef]
[30]	Oğuz, C. and Sezer, M. (2015) Chelyshkov Collocation Method for a Class of Mixed Functional Integro-Differential Equations. Applied Mathematics and Computation, 259, 943-954.[CrossRef]