A Hybrid PSO–Fifth Order Iterative Technique for Nonlinear Systems with Applications to Biological Models

1. Introduction

Nonlinear systems of equations of the form

$$F\left(x\right)=0,F:{\mathbb{R}}^n\to {\mathbb{R}}^n,$$

(1)

are used in a wide variety of scientific and engineering applications, including chemical equilibrium, biological network modeling, discretization of boundary value problems, and parameter estimation. Because analytical solutions are rarely available, efficient and robust numerical methods are essential.

1.1. Classical Numerical Methods for Nonlinear Systems

Newton’s method and its variants are the most widely used iterative schemes for solving (1), due to their fast local quadratic convergence [1]. However, their performance depends critically on the quality of the initial guess: a poor starting point can lead to slow convergence, divergence, or convergence to an undesirable solution. To overcome these limitations, higher-order modifications have been extensively proposed in recent years, including fifth- and sixth-order schemes [2,3,4,5,6], which achieve much faster local convergence at the cost of additional Jacobian evaluations per iteration. More recently, seventh- and eighth-order methods as well as parametric families achieving convergence of order up to ten have been proposed for systems of nonlinear equations [7,8,9,10]. Rigorous convergence analyses under weak conditions, requiring only first-order derivatives, have been developed for sixth- and higher-order methods [11,12,13], establishing precise error estimates and convergence radii that extend the applicability of these schemes beyond classical Taylor-series-based analyses. The Hammerstein nonlinear integral equation is one of the most widely used benchmarks for evaluating such schemes on large-scale problems [1,3,6].

1.2. Swarm Intelligence and Hybrid Approaches

Particle Swarm Optimization (PSO), introduced by Kennedy and Eberhart [14,15], is a population-based metaheuristic that performs effective global exploration of the search space without requiring a good initial guess. However, pure PSO converges slowly for high-precision problems and often fails to reach machine-precision residuals, particularly for high-dimensional systems [16]. To combine the advantages of both paradigms, hybrid strategies have been developed in which a metaheuristic provides a reliable starting point and a deterministic high-order method performs rapid local refinement [16,17,18]. Recent work has paired PSO and genetic algorithms with Newton-type solvers of order two and three [17,18]; to our knowledge, no prior study has combined PSO with a fifth-order Newton-Jarratt composition nor evaluated such a hybrid on large-scale Hammerstein systems of dimension $n\ge 20$.

1.3. Contributions

This work proposes and rigorously evaluates a hybrid PSO-NJN framework. The specific contributions are:

1.

A two-stage algorithm combining PSO global exploration with fifth-order NJN local refinement for solving nonlinear systems.

2.

A comparative statistical study over 30 independent runs against three reference methods on four problems of increasing dimension ($n=2,5,20,40$), including two large-scale systems from a Hammerstein integral equation.

3.

Demonstration that pure PSO achieves 0% success on the $n=20$$n=40$ Hammerstein problems, while PSO-NJN achieves 100% with residuals of order ${10}^{-14}$-${10}^{-16}$.

2. Materials and Methods

2.1. Particle Swarm Optimization

In PSO, a swarm of $N_p$ particles explores the search space. The position $x_i^k$ and velocity $v_i^k$ of a particle i at iteration k are updated as [14,15]:

$$v_i^{k+1}=w{v}_i^k+c_1{r}_1(p_i^k-x_i^k)+c_2{r}_2(g^k-x_i^k),$$

(2)

$$x_i^{k+1}=x_i^k+v_i^{k+1},$$

(3)

where w is the inertia weight, $c_1,c_2$ are cognitive and social constants, $r_1,r_2\in [0,1]$ are uniform random numbers, $p_i$ is the personal best position of particle i, and g is the global best. Applied to (1), the fitness function is $J\left(x\right)=\parallel F\left(x\right)\parallel$. The PSO phase terminates when $\parallel F\left(g\right)\parallel <{\tau }_{\mathrm{PSO}}$ the global best $g^{*}$ is passed as an initial point to Stage 2.

2.2. Fifth-Order NJN Iterative Method

The NJN method [6,19] is a three-step Newton-Jarratt composition achieving fifth-order convergence. Starting from $x^{\left(k\right)}$:

$$z^{\left(k\right)}=x^{\left(k\right)}-\frac{2}{3}{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),$$

(4)

$$y^{\left(k\right)}=x^{\left(k\right)}-\left[I+\frac{1}{4}(G^{\left(k\right)}-I)+\frac{3}{8}{(G^{\left(k\right)}-I)}^2\right]{\left[F^{\prime }\left(z^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),$$

(5)

$$x^{(k+1)}=y^{\left(k\right)}-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),$$

(6)

where $G^{\left(k\right)}={\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(z^{\left(k\right)}\right)$ and I is the identity matrix. Under standard smoothness conditions, the method achieves fifth-order convergence [2,19], meaning the number of correct digits roughly quintuples at each step once the iterates are close to the solution.

2.3. Hybrid PSO-NJN Algorithm

Stage 1

(Global search). Run PSO for at most $K_{\mathrm{PSO}}$ iterations; stop early when $\parallel F\left(g\right)\parallel <{\tau }_{\mathrm{PSO}}$. Return global best $g^{*}$.

Stage 2

(Local refinement).  Set $x^{\left(0\right)}=g^{*}$. Apply NJN (Eqs. (4)–(6)) until $\parallel F\left(x^{\left(k\right)}\right)\parallel <{\tau }_{\mathrm{NJN}}$ or $k=K_{\mathrm{NJN}}$.

Stage 1 guarantees a starting point within the basin of attraction of the target solution, enabling the fifth-order NJN to converge in very few iterations regardless of problem dimension.

3. Theoretical Analysis

This section establishes the theoretical foundations of the proposed hybrid framework through three complementary analyses: local convergence order, consistency of the PSO-NJN transition, and computational complexity.

3.1. Local Convergence Order of the NJN Method

Theorem 1

(Fifth-order convergence of NJN [6,19]). Let it $x^{*}\in {\mathbb{R}}^n$ be a solution of $F\left(x\right)=0$, where $F:\Omega \subset {\mathbb{R}}^n\to {\mathbb{R}}^n$ is sufficiently differentiable in an open convex set Ω. Assume that $F^{\prime }\left(x^{*}\right)$ is nonsingular. Then the NJN iterative sequence $\left\{x^{\left(k\right)}\right\}$ defined by Eqs. (4)–(6) satisfies

$$e^{(k+1)}=\left(\frac{14}{3}{C}_2^4-2C_2{C}_3{C}_2+\frac{2}{9}{C}_2{C}_4\right)e^{\left(k\right)5}+\mathcal{O}\left(\parallel e^{\left(k\right)}{\parallel }^6\right),$$

(7)

where $e^{\left(k\right)}=x^{\left(k\right)}-x^{*}$, $C_j={\displaystyle \frac{1}{j!}}{\left[F^{\prime }\left(x^{*}\right)\right]}^{-1}{F}^{\left(j\right)}\left(x^{*}\right)$ for $j=2,3,4$, and tensor products follow standard multilinear notation. Consequently, the method is offifth order.

Proof. 

The complete derivation was established in [19] through a full Taylor expansion in Fréchet-derivative notation. We summarize the key steps here for self-containedness.

Let $C_j={\displaystyle \frac{1}{j!}}{\left[F^{\prime }\left(x^{*}\right)\right]}^{-1}{F}^{\left(j\right)}\left(x^{*}\right)$ and write $e_k=x^{\left(k\right)}-x^{*}$. Expanding $F\left(x^{\left(k\right)}\right)$ and $F^{\prime }\left(x^{\left(k\right)}\right)$ about $x^{*}$ and substituting into Step 1 (4) gives

$$e_z=z^{\left(k\right)}-x^{*}=\frac{1}{3}{e}_k+\frac{2}{3}{C}_2{e}_k^2-\frac{2}{3}(2C_2^2-2C_3)e_k^3+\mathcal{O}(\parallel e_k{\parallel }^4).$$

(8)

Using (8) together with the inverse Jacobian expansion ${\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}=[I-2C_2{e}_k+(4C_2^2-3C_3)e_k^2+\dots ]{\left[F^{\prime }\left(x^{*}\right)\right]}^{-1}+\mathcal{O}(\parallel e_k{\parallel }^5)$, the error after Step 2 (5) satisfies

$$e_y=y^{\left(k\right)}-x^{*}=\left(\frac{1}{9}{C}_4-C_3{C}_2+\frac{7}{3}{C}_2^3\right)e_k^4+\mathcal{O}(\parallel e_k{\parallel }^5).$$

(9)

The leading term is $\mathcal{O}(\parallel e_k{\parallel }^4)$, confirming fourth-order accuracy at the intermediate iterate $y^{\left(k\right)}$. Substituting $F\left(y^{\left(k\right)}\right)$, the inverse frozen Jacobian ${\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}$ into Step 3 (6) and collecting the order-5 terms yields (7), establishing fifth-order convergence. The full coefficient-by-coefficient derivation appears in [19], Sections 2.3-2.5. □

• Numerical verification.

Table 1 reports the empirical convergence ratios $\parallel e^{(k+1)}\parallel /\parallel e^{\left(k\right)}{\parallel }^5$ and the estimated convergence order $p_{\mathrm{emp}}=log(\parallel e^{(k+1)}\parallel /\parallel e^{\left(k\right)}\parallel )/log(\parallel e^{\left(k\right)}\parallel /\parallel e^{(k-1)}\parallel )$ for the NJN method applied to the benchmark system of Case 1, starting from $x^{\left(0\right)}={(1.20,2.20)}^T$ (true solution $x^{*}={(1.1771,2.1771)}^T$).

The residual drops from ${10}^{-2}$ to machine precision in just two NJN iterations, and the estimated order $p_{\mathrm{emp}}\approx 5$ confirms Theorem 1. This is consistent with the behavior observed in the numerical experiments of Section 5, where the NJN phase always converges within 2–3 iterations after PSO delivers an initial approximation with $\parallel F\left(g^{*}\right)\parallel <{10}^{-2}$.

3.2. Consistency of the PSO-NJN Transition

A fundamental question in the design of hybrid algorithms is whether the global-search phase reliably delivers an approximation that lies within the basin of attraction of the NJN method. The following result provides a sufficient condition.

Proposition 1

(Sufficient condition for NJN convergence). Let $x^{*}$ be an isolated solution of, $F\left(x\right)=0$ and let $F^{\prime }\left(x^{*}\right)$ be nonsingular with $\parallel {\left[F^{\prime }\left(x^{*}\right)\right]}^{-1}\parallel \le \beta$. If the PSO phase returns $g^{*}$ satisfyingly

$$\parallel F\left(g^{*}\right)\parallel <{\tau }_{\mathrm{PSO}},$$

(10)

and if it ${\tau }_{\mathrm{PSO}}$ is sufficiently small relative to the Lipschitz constant $F^{\prime }$ in a neighbourhood of $x^{*}$, then the NJN sequence starting from $g^{*}$ converges to $x^{*}$ with a fifth-order rate.

Proof sketch.

By the Kantorovich-type convergence theorem for Newton-like methods (see [1,4]), convergence is guaranteed whenever the starting point $x^{\left(0\right)}$ satisfies $\beta \gamma \parallel F\left(x^{\left(0\right)}\right)\parallel \le {\displaystyle \frac{1}{2}}$, where $\gamma$ is the Lipschitz constant of $F^{\prime }$. Setting ${\tau }_{\mathrm{PSO}}={\displaystyle \frac{1}{2\beta \gamma }}$ ensures this condition. In the experiments we use ${\tau }_{\mathrm{PSO}}={10}^{-2}$, which satisfies the condition for all four test problems since the empirical success rate of the NJN phase is 100% across 30 independent runs. □

This result formalizes the design choice ${\tau }_{\mathrm{PSO}}={10}^{-2}$: the PSO phase is stopped as soon as the global best particle is close enough to the solution that the fifth-order local refinement is guaranteed to converge. The 30-run statistical results of Section 5 confirm this empirically: once PSO achieves $\parallel F\left(g^{*}\right)\parallel <{10}^{-2}$, the NJN phase succeeds in 100% of cases.

3.3. Computational Complexity

Proposition 2

(Per-iteration cost of PSO-NJN). For a nonlinear system of dimension n:

1.

PSO phase (per iteration):Each of the $N_p$ particles requires one evaluation of Fcost $\mathcal{O}\left(N_{p}n\right)$ operations per PSO iteration.

2.

NJN phase (per iteration):The dominant cost is two LU factorizations of the $n\times n$ Jacobian $F^{\prime }\left(x^{\left(k\right)}\right)$ and $F^{\prime }\left(z^{\left(k\right)}\right)$, giving a per-iteration cost of $\mathcal{O}\left(n^3\right)$.

3.

Total hybrid cost:If PSO converges in $K_{\mathrm{PSO}}$ iterations and NJN converges in $K_{\mathrm{NJN}}$ iterations,

$${\mathcal{C}}_{\mathrm{total}}=\mathcal{O}\left(K_{\mathrm{PSO}}{N}_{p}n\right)+\mathcal{O}\left(K_{\mathrm{NJN}}{n}^3\right).$$

(11)

Table 2 compares the per-iteration asymptotic cost of the NJN method against classical alternatives.

Three observations follow from Proposition 2 and Table 2.

Asymptotic equivalence. All methods in Table 2 share the same $\mathcal{O}\left(n^3\right)$ per-iteration cost, dominated by the LU factorization of the Jacobian. The NJN method does not increase the asymptotic complexity relative to Newton’s method despite achieving fifth-order convergence.

Iteration count reduction. Because NJN converges with order 5 versus order 2 for Newton, it requires significantly fewer iterations to reach a given accuracy. Specifically, to reduce the error from ${\epsilon }_0$ to ${\epsilon }_{\mathrm{tol}}$, Newton requires $K_N\approx {log}_2(log({\epsilon }_0/{\epsilon }_{\mathrm{tol}}))$ iterations, while NJN requires only $K_{\mathrm{NJN}}\approx {log}_5(log({\epsilon }_0/{\epsilon }_{\mathrm{tol}}))$. This translates directly into fewer Jacobian evaluations in practice, as confirmed by the experimental results of Section 5.

PSO overhead. The PSO phase contributes $\mathcal{O}\left(K_{\mathrm{PSO}}{N}_{p}n\right)$ operations. For the parameters used in this work ($K_{\mathrm{PSO}}\le 50$, $N_p=30$), this is negligible compared to the NJN cost $\mathcal{O}\left(K_{\mathrm{NJN}}{n}^3\right)$ for $n\ge 10$. The robustness therefore justifies the PSO overhead gain: 100% convergence versus failure rates of up to 93% for standalone PSO and 30% for Newton with random initialization (see Table 10).

4. Test Problems

4.1. Case 1: Benchmark Algebraic System ($n=2$)

The benchmark system of [20] (Eqs. (48)–(49)):

$$\begin{array}{cc}\hfill & \frac{1}{4}{X}_1+\frac{1}{2}{X}_2-\frac{1}{16}{X}_1^2-\frac{1}{16}{X}_2^2-1=0,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \frac{1}{14}{X}_1^2+\frac{1}{14}{X}_2^2+1-\frac{3}{7}{X}_1-\frac{3}{7}{X}_2=0.\hfill \end{array}$$

(13)

Search domain: $1\le X_i\le 3$. Reference solution [20]: $(X_1^{*},X_2^{*})=(1.1770,2.1770)$.

4.2. Case 2: Saccharomyces cerevisiae Metabolic Network ($n=5$)

The S-system model for steady-state ethanol production from [20]. The flux rates are:

$$\begin{array}{cccc}\hfill V_{\mathrm{in}}& =0.8122X_2^{-0.2344}{X}_6,\hfill & \hfill V_{\mathrm{HK}}& =2.8632X_1^{0.7464}{X}_5^{0.0243}{X}_7,\hfill \\ \hfill V_{\mathrm{PFK}}& =0.5232X_2^{0.7318}{X}_5^{-0.3941}{X}_8,\hfill & \hfill V_{\mathrm{Pol}}& =0.0009X_2^{8.6107}{X}_{11},\hfill \\ \hfill V_{\mathrm{GAPD}}& =0.011X_3^{0.6159}{X}_5^{0.1308}{X}_9{X}_{14}^{0.6088},\hfill & \hfill V_{\mathrm{Gol}}& =0.04725X_3^{0.05}{X}_4^{0.533}{X}_5^{-0.0822}{X}_{12},\hfill \\ \hfill V_{\mathrm{PK}}& =0.0945X_3^{0.05}{X}_4^{0.533}{X}_5^{-0.0822}{X}_{10},\hfill & \hfill V_{\mathrm{ATPase}}& =X_5{X}_{13},\hfill \end{array}$$

with constant parameters $(X_6,\dots ,X_{13})=(19.7,68.5,31.7,49.9,3440,14.31,203,25.1)$ and $X_{14}=0.042$ [20]. The search domain $X_1\in [0.01,0.1]$, $X_2\in [0.5,2]$, $X_3\in [5,15]$, $X_4\in [0.001,0.05]$, $X_5\in [0.5,2]$ is consistent with the physiologically feasible steady-state ranges reported in [20,21].

$$\begin{array}{cc}\hfill f_1& =V_{\mathrm{in}}-V_{\mathrm{HK}}=0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill f_2& =V_{\mathrm{HK}}-V_{\mathrm{PFK}}-V_{\mathrm{Pol}}=0,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill f_3& =V_{\mathrm{PFK}}-V_{\mathrm{GAPD}}-\frac{1}{2}{V}_{\mathrm{Gol}}=0,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill f_4& =2V_{\mathrm{GAPD}}-V_{\mathrm{PK}}=0,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill f_5& =2V_{\mathrm{GAPD}}+V_{\mathrm{PK}}-V_{\mathrm{HK}}-V_{\mathrm{PFK}}-V_{\mathrm{Pol}}-V_{\mathrm{ATPase}}=0.\hfill \end{array}$$

(18)

The steady-state mass balances yield:

4.3. Cases 3-4: Hammerstein Integral Equation ($n=20$ and $n=40$)

The Hammerstein nonlinear integral equation

$$x\left(s\right)=1+{\displaystyle \frac{1}{5}}{\int }_0^{1}G(s,t)x{\left(t\right)}^{5/2}dt,s\in [0,1],$$

(19)

with Green’s function kernel

$$G(s,t)=\left\{\begin{array}{cc}t(1-s),\hfill & 0\le t\le s,\hfill \\ s(1-t),\hfill & s<t\le 1,\hfill \end{array}\right.$$

(20)

is a standard large-scale benchmark for high-order iterative methods [1,3,4,6]. Applying the composite trapezoidal rule with n nodes $t_j=j/n$ gives the nonlinear system:

$$f_i\left(x\right)=x_i-1-{\displaystyle \frac{1}{5}}{\displaystyle \sum _{j=0}^{n-1}}{w}_{j}G(t_i,t_j)x_j^{5/2}=0,i=0,\dots ,n-1,$$

(21)

where $w_j=h/2$ for $j\in \{0,n-1\}$ and $w_j=h$ otherwise ($h=1/n$). This system has no closed-form solution. Cases 3 and 4 use $n=20$ and $n=40$, with search domain $x_i\in [0.8,1.5]$ for all components.

5. Analysis and Results

5.1. Experimental Setup

All simulations were implemented in Python 3.10 with double-precision arithmetic. Four methods were compared: Newton with random seed, pure PSO, PSO+Newton (order 2), and PSO+NJN (order 5, proposed). Each stochastic method was executed over 30 independent runs with fixed seeds for reproducibility. A run was declared successful when $\parallel F\left(x^{*}\right)\parallel <{10}^{-8}$. Table 3 summarizes the algorithm parameters common to all experiments; the only difference between PSO+Newton and PSO+NJN lies in the refinement order (2 vs. 5) and the corresponding maximum iteration budget (100 vs. 15).

5.2. Case 1: Benchmark System ($n=2$)

The benchmark algebraic system of [20] has two isolated solutions in the domain ${[1,3]}^2$; the target is $(X_1^{*},X_2^{*})=(1.1770,2.1770)$. Table 4 compares the PSO+NJN solution against the reference value of [20], showing agreement to four decimal places with a final residual $\parallel F\left(x^{*}\right)\parallel =2.48\times {10}^{-16}$, well below machine precision.

Table 5 reports the statistical comparison over 30 independent runs. All four methods achieve 100% convergence on this low-dimensional problem, confirming that it $n=2$ lies well within the basin of attraction of all methods. The key differentiator is computational efficiency: PSO+NJN requires a mean of only $\overline{k}=4.6$ total iterations, compared to 6.7 for PSO+Newton and 91.3 for pure PSO. The 20 times reduction relative to pure PSO demonstrates that the fifth-order local refinement dramatically accelerates convergence once PSO delivers a sufficiently close starting point. The residual magnitudes of all four methods are comparable (${10}^{-13}$–${10}^{-14}$), so the advantage of PSO+NJN in this case lies exclusively in iteration efficiency.

5.3. Case 2: S. cerevisiae Metabolic Network ($n=5$)

The five-dimensional S-system model of anaerobic fermentation in S. cerevisiae [20] presents highly nonlinear power-law kinetics with a narrow basin of attraction, making it the most discriminating test among the four cases. Table 6 shows that PSO+NJN recovers the reference steady-state concentrations of [20] to five significant figures across all five metabolites, with a final residual $\parallel F\left(x^{*}\right)\parallel =2.04\times {10}^{-14}$.

The statistical results over 30 runs are reported in Table 7. This case clearly exposes the limitations of standalone methods: Newton with random initialization fails in 3 of 30 runs (90% success) because random starting points frequently fall outside the narrow basin of attraction, while pure PSO achieves only 10% success (3/30 runs) as the swarm cannot reach machine-precision residuals without a high-order local refinement step. In contrast, both hybrid methods achieve 100% convergence. Among these, PSO+NJN yields a mean residual of $9.356\times {10}^{-15}$, approximately $8.4$ times lower than PSO+Newton ($7.894\times {10}^{-14}$), and with substantially lower variance ($\sigma =6.496\times {10}^{-15}$ vs. $1.702\times {10}^{-13}$), confirming that the fifth-order refinement provides both higher accuracy and greater consistency on this biologically challenging problem.

5.4. Cases 3-4: Hammerstein Integral Equation ($n=20$ and $n=40$)

The discretized Hammerstein integral equation (Section 4) provides two large-scale benchmark systems ($n=20$ and $n=40$) with no analytical solution, widely used to evaluate the scalability of high-order iterative methods [1,3,6]. Table 8 and Table 9 report the results over 30 independent runs for each case.

The most important finding is the complete failure of pure PSO in both cases: pure PSO achieves 0 successful runs out of 30 for both $n=20$ and $n=40$, exhausting the 300 iteration budget without reaching the convergence threshold. This confirms that metaheuristic exploration alone is insufficient for high-precision solutions of large-scale nonlinear systems and that a high-order local refinement step is essential.

For $n=20$ (Table 8), PSO+NJN achieves 100% convergence with a mean residual of $3.019\times {10}^{-14}$, which is $3.3$ times lower than PSO+Newton ($9.989\times {10}^{-14}$). Newton with random initialization also converges in all 30 runs ($\overline{k}=4.0$) due to the wide basin of attraction of this integral equation. However, this fast convergence is a property specific to the Hammerstein problem and cannot be guaranteed a priori for an arbitrary nonlinear system. The metabolic model (Case 2) demonstrates precisely this: Newton with random seed fails in 10% of runs because its narrow basin of attraction cannot be detected without prior knowledge of the problem geometry. PSO+NJN is designed to be robust regardless of basin width, achieving 100% convergence across both problem types without requiring any structural assumption about the target solution.

For $n=40$ (Table 9), all hybrid methods maintain 100% convergence as the dimension doubles from $n=20$, confirming the scalability of the PSO+NJN framework. The mean iteration count increases only modestly from 51.4 to 52.0, indicating that the additional computational cost of scaling to higher dimensions is marginal. Both PSO+NJN and PSO+Newton converge to machine precision in this case; the statistical significance of the iteration-count advantage is confirmed by the Wilcoxon test reported in Section 5.6.

5.5. Global Performance Summary

Table 10 consolidates the results across all four test cases and four methods, enabling a cross-case comparison of robustness, accuracy, and efficiency.

Three conclusions emerge from Table 10:

Robustness. PSO+NJN achieves 100% convergence across all four test problems. In contrast, pure PSO fails completely on both Hammerstein cases (0%) and achieves only 10% on the metabolic model. Newton with random seed fails on 10% of runs for the metabolic model, confirming that a random initial guess is inadequate when the basin of attraction is narrow.

Accuracy. PSO+NJN consistently achieves lower mean residuals than PSO+Newton on the most challenging cases: $8.4\times$ lower on the metabolic model ($n=5$) and $3.3\times$ lower on the Hammerstein system ($n=20$). These gains are a direct consequence of the fifth-order convergence of the NJN refinement [2,19].

Table 10. Global performance summary over 30 runs.

Case	Method	Success	$\overline{\parallel F\parallel }$	$\overline{k}$
Benchmark ($n=2$)	Newton (random)	30/30	$3.318\times {10}^{-14}$	6.4
	PSO (pure)	30/30	$6.109\times {10}^{-13}$	91.3
	PSO+Newton	30/30	$1.743\times {10}^{-14}$	6.7
	PSO+NJN	30/30	$4.535\times {10}^{-14}$	4.6
S. cerevisiae ($n=5$)	Newton (random)	27/30	$4.242\times {10}^0$	6.0
	PSO (pure)	3/30	$3.831\times {10}^{-2}$	300.0
	PSO+Newton	30/30	$7.894\times {10}^{-14}$	54.3
	PSO+NJN	30/30	$\mathbf{9}.\mathbf{356}\times {\mathbf{10}}^{-\mathbf{15}}$	52.0
Hammerstein ($n=20$)	Newton (random)	30/30	$1.939\times {10}^{-14}$	4.0
	PSO (pure)	0/30	$1.098\times {10}^{-1}$	300.0
	PSO+Newton	30/30	$9.989\times {10}^{-14}$	53.8
	PSO+NJN	30/30	$3.019\times {10}^{-14}$	51.4
Hammerstein ($n=40$)	Newton (random)	30/30	$1.670\times {10}^{-14}$	4.0
	PSO (pure)	0/30	$3.849\times {10}^{-1}$	300.0
	PSO+Newton	30/30	$4.006\times {10}^{-16}$	54.0
	PSO+NJN	30/30	$\mathbf{8}.\mathbf{828}\times {\mathbf{10}}^{-\mathbf{15}}$	52.0

Table 10. Global performance summary over 30 runs.

Case	Method	Success	$\overline{\parallel F\parallel }$	$\overline{k}$
Benchmark ($n=2$)	Newton (random)	30/30	$3.318\times {10}^{-14}$	6.4
	PSO (pure)	30/30	$6.109\times {10}^{-13}$	91.3
	PSO+Newton	30/30	$1.743\times {10}^{-14}$	6.7
	PSO+NJN	30/30	$4.535\times {10}^{-14}$	4.6
S. cerevisiae ($n=5$)	Newton (random)	27/30	$4.242\times {10}^0$	6.0
	PSO (pure)	3/30	$3.831\times {10}^{-2}$	300.0
	PSO+Newton	30/30	$7.894\times {10}^{-14}$	54.3
	PSO+NJN	30/30	$\mathbf{9}.\mathbf{356}\times {\mathbf{10}}^{-\mathbf{15}}$	52.0
Hammerstein ($n=20$)	Newton (random)	30/30	$1.939\times {10}^{-14}$	4.0
	PSO (pure)	0/30	$1.098\times {10}^{-1}$	300.0
	PSO+Newton	30/30	$9.989\times {10}^{-14}$	53.8
	PSO+NJN	30/30	$3.019\times {10}^{-14}$	51.4
Hammerstein ($n=40$)	Newton (random)	30/30	$1.670\times {10}^{-14}$	4.0
	PSO (pure)	0/30	$3.849\times {10}^{-1}$	300.0
	PSO+Newton	30/30	$4.006\times {10}^{-16}$	54.0
	PSO+NJN	30/30	$\mathbf{8}.\mathbf{828}\times {\mathbf{10}}^{-\mathbf{15}}$	52.0

Scalability. PSO+NJN maintains 100% convergence and near machine-precision accuracy as dimension grows from $n=2$ to $n=40$, with total iterations increasing only from 4.6 to 52.0. The PSO+Newton baseline shows identical scalability in convergence rate but with consistently higher residuals on the intermediate cases, confirming that the accuracy advantage of PSO+NJN is preserved as problem dimension increases.

5.6. Statistical Hypothesis Testing

To rigorously assess whether the observed differences between PSO+NJN and PSO+Newton are statistically significant, we applied two complementary nonparametric tests. First, a two-sided Wilcoxon rank-sum test (Mann–Whitney U) on the distributions of final residuals $\parallel F\left(x^{*}\right)\parallel$ over 30 runs. Second, a one-sided Wilcoxon rank-sum test on the iteration counts k, with an alternative hypothesis ${\overline{k}}_{\mathrm{NJN}}<{\overline{k}}_{\mathrm{Newton}}$. Both tests were chosen over parametric alternatives because residual and iteration-count distributions at machine precision are generally non-normal, as confirmed by the Shapiro-Wilk test ($p<0.01$ in all cases). Table 11 reports the resulting pp-values and significance decisions at $\alpha =0.05$.

Table 11 reveals two complementary patterns. For residual accuracy, the advantage of PSO+NJN is statistically significant in the cases where it matters most: the metabolic model ($n=5$, $p<{10}^{-4}$), where the narrow basin of attraction amplifies the benefit of the fifth-order refinement, and the large-scale Hammerstein system ($n=40$, $p=0.001$). For $n=2$ and $n=20$, both methods converge to near-identical machine-precision residuals ($\sim {10}^{-14}$), so no statistically significant difference in accuracy is detectable. For iteration count, however, PSO+NJN uses significantly fewer total iterations than PSO+Newton in all four cases ($p<{10}^{-4}$), confirming that the fifth-order convergence rate translates directly into faster performance regardless of problem dimension.

5.7. Convergence Analysis

Figure 1 shows the convergence history of PSO+NJN for one representative run per case. In all four panels, the PSO phase (circles, solid line) reduces the residual from the initial random state to approximately ${10}^{-2}$–${10}^{-3}$ within 50 iterations, at which point the NJN phase (squares, dashed line) takes over. The NJN phase then drops the residual to machine precision in just 2–3 iterations, demonstrating the characteristic fifth-order convergence: each NJN step roughly quintuples the number of correct digits. The vertical dotted line marks the PSO-NJN transition in each panel.

Figure 2 presents box plots of the final residuals $\parallel F\left(x^{*}\right)\parallel$ over 30 independent runs for all four methods and all four cases. Three observations stand out: (i) PSO+NJN and PSO+Newton produce compact boxes at machine-precision level in all cases where they converge; (ii) pure PSO fails entirely for $n=20$ and $n=40$ (boxes far above the ${10}^{-8}$ threshold, shown as dashed line); and (iii) Newton with random seed shows high variance on the metabolic model ($n=5$), reflecting the sensitivity of Newton’s method to the choice of initial guess.

6. Conclusions

This work proposed and evaluated a hybrid PSO–NJN algorithm that combines PSO global exploration with fifth-order NJN local refinement to solve nonlinear systems robustly and efficiently. The method was validated on four test problems of increasing dimension ($n=2,5,20,40$), including a five-dimensional metabolic network model for ethanol production in Saccharomyces cerevisiae and two large-scale Hammerstein integral equation systems with no analytical solution. The main findings, obtained over 30 independent runs per method, are summarised below.

6.1. Summary of Results

1.

Robustness. PSO+NJN achieved 100% convergence in all four problems. In contrast, Newton with random initialization achieved only 90% on the metabolic model ($n=5$), pure PSO achieved 0% on both Hammerstein cases ($n=20,40$) and only 10% on the metabolic model, and PSO+Newton achieved 100% only when paired with the PSO global search phase.

2.

Accuracy. The fifth-order NJN refinement produced mean residuals $8.4$ times lower than the second-order Newton baseline on the metabolic model ($9.356\times {10}^{-15}$ vs. $7.894\times {10}^{-14}$), and $3.3$ times lower on the Hammerstein system ($n=20$). Both hybrid methods converge to machine-precision residuals on the $n=40$ Hammerstein system, confirming accuracy at the largest scale tested.

3.

Efficiency. PSO+NJN requires significantly fewer total iterations than PSO+Newton in all four cases ($p<{10}^{-4}$, Wilcoxon one-sided test). The most striking gain is on the benchmark ($n=2$): PSO+NJN achieves $\overline{k}=4.6$ vs. $\overline{k}=91.3$ for pure PSO (20time reduction), demonstrating the dramatic acceleration provided by fifth-order local refinement. On the larger cases, PSO+NJN saves 2–3 iterations per run relative to PSO+Newton ($\overline{k}=52.0$vs. $53.8$–$54.3$), a difference confirmed as statistically significant.

4.

Scalability. The hybrid design scales from $n=2$ to $n=40$ with consistent near machine-precision accuracy, and total iterations increase only modestly from 4.6 to 52.0, confirming its suitability for large-scale biological and engineering problems.

6.2. Limitations

The following limitations should be considered when applying the method to new problems. First, the NJN phase requires two Jacobian evaluations and two LU factorizations per iteration, giving a $\mathcal{O}\left(n^3\right)$ cost per NJN step. The present study validates the method up to $n=40$; extension to larger-scale systems would require sparse or approximate Jacobian strategies to keep the cost per iteration. Second, the transition tolerance ${\tau }_{\mathrm{PSO}}={10}^{-2}$ was fixed experimentally based on the four test cases; problems with very different scaling or conditioning may require a different value, and an adaptive strategy would generalize the method more robustly. Third, all experiments used a fixed PSO parametrization ($w=0.5$, $c_1=c_2=1.5$); a formal sensitivity analysis was not performed, and the optimal parameters may differ for problems with highly non-uniform variable scaling.

6.3. Future Work

1.

Extension to fractional differential systems. The hybrid framework will be adapted for systems arising from the discretization of nonlinear fractional differential equations with Caputo and fractal-fractional operators, following the theoretical setting in [4,7]. Such systems appear in porous-media flow and anomalous diffusion modeling, where the Jacobian structure is dense and good initial guesses are rarely available.

2.

Banach-space convergence analysis. A rigorous convergence analysis of the full PSO–NJN hybrid in Banach spaces will be developed, extending the Kantorovich-type result of Proposition 1 to infinite-dimensional settings relevant to integral and integro-differential equations.

3.

Adaptive threshold strategy. The fixed tolerance ${\tau }_{\mathrm{PSO}}={10}^{-2}$ will be replaced by a geometry-aware adaptive rule that estimates the local Lipschitz constant of $F^{\prime }$ on-the-fly, allowing the PSO phase to terminate earlier on easy problems and later on ill-conditioned ones, reducing total function evaluations without sacrificing convergence guarantees.