Paper deep dive
Physics-Informed Stochastic Configuration Machine: A Backpropagation-Free Neural Network with Fast Training for Nonlinear Differential Equations
Yuehao Song, Zhong Chen, Lihui Cen, Liang Wu, Kai Zhang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/28/2026, 3:51:29 AM
Summary
The paper introduces the Physics-Informed Stochastic Configuration Machine (PI-SCM), a backpropagation-free neural network framework designed to solve forward and inverse problems in nonlinear differential equations. Unlike traditional Physics-Informed Neural Networks (PINNs) that rely on gradient descent and automatic differentiation, PI-SCM uses analytical evaluation of local Jacobians to linearize the physical loss. This allows for the explicit determination of optimal network weights via generalized linear least squares solvers, significantly accelerating training and avoiding non-convex optimization traps. The framework includes three progressive algorithms (PI-SC-I, PI-SC-II, PI-SC-III) and demonstrates high-fidelity accuracy and robust parameter identification.
Entities (10)
Relation Signals (10)
PI-SCM → isa → Backpropagation-free Neural Network
confidence 98% · we propose the Physics-Informed Stochastic Configuration Machine (PI-SCM), a novel backpropagation-free framework
PI-SC-III → ispartof → PI-SCM
confidence 95% · global updating (PI-SC-III)
PI-SC-I → ispartof → PI-SCM
confidence 95% · progressive algorithmic suite comprising localized construction (PI-SC-I)
PI-SC-II → ispartof → PI-SCM
confidence 95% · sliding-window updating (PI-SC-II)
PINNs → relieson → Backpropagation
confidence 95% · reliance on backpropagation-based gradient descent
PI-SCM → solves → Nonlinear Differential Equations
confidence 95% · framework for both forward and inverse problems in differential equations
PI-SCM → uses → Analytical Jacobians
confidence 92% · core mathematical contribution lies in the analytical evaluation of local Jacobians
PI-SCM → uses → Least Squares Solvers
confidence 92% · explicit determination of optimal network weights via a sequence of generalized linear least squares solvers
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:While Physics-Informed Neural Networks (PINNs) have emerged as a transformative paradigm for solving complex differential equations, their reliance on backpropagation-based gradient descent and automatic differentiation (AD) imposes significant computational bottlenecks and severe non-convex optimization challenges. To overcome these fundamental limitations, we propose the Physics-Informed Stochastic Configuration Machine (PI-SCM), a novel backpropagation-free framework for both forward and inverse problems in differential equations. The core mathematical contribution lies in the analytical evaluation of local Jacobians for nonlinear differential operators, which facilitates a linearized representation of the physical loss and projects it into a unified, linearized algebraic subspace. This reformulation allows for the explicit determination of optimal network weights via a sequence of generalized linear least squares solvers, effectively bypassing the iterative traps of traditional nonlinear optimizers. We develop a progressive algorithmic suite comprising localized construction (PI-SC-I), sliding-window updating (PI-SC-II), and global updating (PI-SC-III), and rigorously establish their universal approximation properties. Extensive experiments demonstrate that PI-SCM achieves high-fidelity predictive accuracy and robust parameter identification while accelerating the training process by orders of magnitude compared to standard PINNs. Our work provides a highly efficient and scalable foundation for next-generation, real-time Scientific Machine Learning applications.
Tags
Links
- Source: https://arxiv.org/abs/2608.26549v1
- Canonical: https://arxiv.org/abs/2608.26549v1
Trouble viewing inline? Open PDF directly →
Full Text
114,299 characters extracted from source content.
Expand or collapse full text
Physics-Informed Stochastic Configuration Machine: A Backpropagation-Free Neural Network with Fast Training for Nonlinear Differential Equations Yuehao Song Zhong Chen Lihui Cen Liang Wu and Kai Zhang †thanks: Yuehao Song, Zhong Chen, and Lihui Cen are with the School of Automation, Central South University, 932 South Lushan Road, Changsha 410083, China (e-mails: yh.song@csu.edu.cn; zhongchen@csu.edu.cn; lhcen@csu.edu.cn).†thanks: Liang Wu is with Johns Hopkins University, Baltimore, MD 21218, USA (e-mail: wliang14@jh.edu).†thanks: Kai Zhang is with the State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China (e-mail: zhangkai@iwhr.com).†thanks: Corresponding author: Lihui Cen (e-mail: lhcen@csu.edu.cn). Abstract While Physics-Informed Neural Networks (PINNs) have emerged as a transformative paradigm for solving complex differential equations, their reliance on backpropagation-based gradient descent and automatic differentiation (AD) imposes significant computational bottlenecks and severe non-convex optimization challenges. To overcome these fundamental limitations, we propose the Physics-Informed Stochastic Configuration Machine (PI-SCM), a novel backpropagation-free framework for both forward and inverse problems in differential equations. The core mathematical contribution lies in the analytical evaluation of local Jacobians for nonlinear differential operators, which facilitates a linearized representation of the physical loss and projects it into a unified, linearized algebraic subspace. This reformulation allows for the explicit determination of optimal network weights via a sequence of generalized linear least squares solvers, effectively bypassing the iterative traps of traditional nonlinear optimizers. We develop a progressive algorithmic suite—comprising localized construction (PI-SC-I), sliding-window updating (PI-SC-I), and global updating (PI-SC-I)—and rigorously establish their universal approximation properties. Extensive experiments demonstrate that PI-SCM achieves high-fidelity predictive accuracy and robust parameter identification while accelerating the training process by orders of magnitude compared to standard PINNs. Our work provides a highly efficient and scalable foundation for next-generation, real-time Scientific Machine Learning applications. Index Terms: Nonlinear Differential Equations, Scientific Machine Learning, Physics-Informed Neural Networks, Stochastic Configuration Networks, Backpropagation-free. I Introduction Nonlinear differential equations, encompassing both ordinary differential equations (ODEs) and partial differential equations (PDEs), serve as the fundamental mathematical language for describing complex dynamical systems across diverse scientific and engineering disciplines. Traditionally, solving these governing equations relies on classical numerical techniques, such as the finite element method (FEM) and high-order Runge-Kutta solvers. While mathematically rigorous, these traditional solvers often encounter significant computational bottlenecks when confronted with the “curse of dimensionality” or complex geometric domains. To circumvent these limitations, Scientific Machine Learning (SciML) has emerged as a transformative computational paradigm. As a burgeoning field, SciML seeks to integrate machine learning-derived methodologies into traditional engineering frameworks, thereby leveraging domain knowledge and known physical constraints to inform and regularize the learning process [1]. At the forefront of this revolution are Physics-Informed Neural Networks (PINNs) [2]. By embedding the physical governing equations directly into the loss function of deep neural networks, PINNs utilize established physical laws to drive the optimization process, enabling robust, mesh-free approximations. Consequently, PINNs have seen widespread and highly successful applications across varied engineering domains. For instance, Schiassi et al. [3] employed PINNs to learn the state–costate dynamics satisfying optimality conditions, enabling efficient solutions to complex optimal control problems such as orbital transfers. Shukla et al. [4] utilized PINNs to identify and characterize surface-breaking cracks in metal plates. Berkhahn et al. [5] utilized PINNs to model COVID-19 infection and hospitalization scenarios. Diao et al. [6] applied PINNs to tackle intricate boundary conditions in multi-material problems within solid mechanics. However, despite these remarkable successes, the standard PINN architecture is fundamentally constrained by severe computational bottlenecks. Relying heavily on backpropagation-based gradient descent algorithms—such as Adam and L-BFGS [7]—to minimize a composite loss function, this optimization paradigm encounters two critical mathematical challenges when applied to highly nonlinear dynamical systems. First, from an optimization perspective, the incorporation of nonlinear differential operators induces a highly nonconvex optimization landscape [8]. This intricate topology frequently triggers severe gradient pathologies, where competing gradient magnitudes and directions from different loss components hinder effective weight updates. Consequently, iterative gradient-based optimizers are highly susceptible to becoming trapped in sub-optimal local minima, manifesting as residual stagnation and degraded predictive accuracy. Second, from a computational perspective, the inherent mechanics of backpropagation impose a severe efficiency bottleneck. Evaluating the physical loss necessitates computing derivatives with respect to the spatio-temporal coordinates via automatic differentiation (AD). Continuously traversing the deep computational graph backward to update weights entails a massive memory footprint and incurs an exorbitant computational cost per iteration [8]. This fundamental limitation results in slow convergence speeds and prohibitive overall training times, rendering standard PINNs impractical for real-time applications or large-scale parametric studies. Therefore, developing a robust, backpropagation-free learning paradigm that bypasses non-convex optimization traps and the heavy computational overhead imposed by iterative gradient updates and AD to guarantee numerical stability and efficiency remains an urgent challenge in the SciML community. In this context, randomized learning methods have emerged as a promising alternative, offering a fundamental departure from the gradient descent paradigms used in conventional neural networks. These models typically rely on a two-step training strategy: first, hidden-layer weights and biases are randomly assigned; second, the output weights are determined through analytical computation. Early manifestations of this concept include the Random Vector Functional Link Network (RVFLN) [9] and feed forward neural networks with random weights proposed by Schmidt et al. [10]. Subsequently, Extreme Learning Machines (ELM) [11] gained widespread traction in machine vision and industrial processes due to their remarkable learning speed and ease of implementation. To address the issues of node quality and network scale configuration inherent in static randomized networks, incremental construction methods such as Incremental Random Weight Neural Networks (IRWNN) [12] were developed to build architectures node by node. Building upon these foundations, Stochastic Configuration Networks (SCNs) [13] introduced a rigorous supervisory mechanism for parameter allocation. Unlike earlier randomized models, SCNs randomly configure hidden parameters within data-dependent adaptive ranges under an inequality constraint, which guarantees the universal approximation property. By combining this constructive strategy with an analytical linear least squares solver for output weights, SCNs enable exceptionally fast modeling while maintaining strong generalization capabilities. The efficacy of this framework has inspired numerous advancements, including DeepSCN [14], algorithmic optimizations via genetic algorithms [15] or Monte Carlo tree search [16, 17], and successful applications across diverse domains [18, 19, 20, 21, 22]. However, standard SCNs are inherently purely data-driven architectures, limiting their direct applicability in SciML scenarios where high-quality data is scarce and systems are governed by strict physical laws. To bridge this gap, Xu et al. [23] introduced a Physics-Informed Stochastic Configuration Network (PISCN). Yet, this framework suffers from two fundamental limitations. First, during the critical hidden node configuration phase, the node selection process does not incorporate physical information constraints. Second, the subsequent training process reverts to traditional gradient descent, abandoning the analytical linear least squares solver of the original SCN framework. This significantly degrades its modeling efficiency and negates the primary computational advantage of randomized networks. In addition, operator-theoretic learning approaches, such as Koopman operator theory, provide computationally efficient methods for modeling large-scale nonlinear systems. Within this framework, Wu et al. proposed a parallelizable multi-step EDMD method that decomposes the least-squares identification problem across prediction horizons and state coordinates, enabling efficient offline modeling [24]. However, the method remains data-driven and does not incorporate physical constraints. To retain the computational efficiency of SCNs in unsupervised and semi-supervised physics-informed learning, this paper proposes a novel, backpropagation-free framework: the Physics-Informed Stochastic Configuration Machine (PI-SCM). Our approach fundamentally restructures the optimization paradigm for solving complex differential equations. The main contributions of this work are summarized as follows: 1) At the core of PI-SCM is a local linearization of the nonlinear differential operator. The operator Jacobians with respect to the state variables and their derivatives are evaluated analytically, allowing the physical residual to be expressed in terms of increments in the network approximation. The resulting linearized physical equations are then combined with the observational constraints to form a unified linearized algebraic system. 2) At each construction or update step, the output coefficients are computed using a pseudoinverse-based linear least-squares solver. Because the required derivatives are evaluated analytically, the procedure does not rely on automatic differentiation. PI-SCM therefore provides a backpropagation-free training mechanism while retaining the governing equations as explicit constraints on the learned solution. 3) The framework comprises three progressive algorithms. PI-SC-I performs a localized update when a new hidden node is added, PI-SC-I updates the nodes within a sliding window, and PI-SC-I globally updates all constructed nodes. These algorithms offer different trade-offs between computational cost and approximation accuracy. Under the stated assumptions and acceptance conditions, their update rules ensure monotonic residual reduction, while the accompanying universal approximation property provides a theoretical foundation for the framework. 4) The same linearized formulation is also extended to inverse problems. By updating physical parameters together with the network approximation, PI-SCM supports state reconstruction and parameter identification without introducing a separate backpropagation-based optimization procedure. The remainder of this paper is organized as follows. Section I introduces the preliminary concepts, defining the network architecture and formulating the forward and inverse physical problems. Section I establishes the theoretical and algorithmic foundations of the PI-SCM framework, rigorously proving its universal approximation property and outlining the hyperparameter configuration strategy. Section IV presents comprehensive numerical experiments to validate the predictive accuracy and computational efficiency of the proposed framework against backpropagation-based PINNs on benchmark physical systems. Finally, Section V summarizes the core findings of this study and discusses potential avenues for future research. I Preliminaries of the PI-SCM Framework and Problem Statement I-A Network Architecture of PI-SCM The fundamental distinction between the proposed Physics-Informed Stochastic Configuration Machine (PI-SCM) and conventional backpropagation-based networks lies in its constructive, backpropagation-free architecture. Unlike standard deep neural networks with fixed topologies, PI-SCM constructs a single-hidden-layer neural network incrementally node-by-node. Suppose the network has generated L hidden nodes. For a given input coordinate vector =(x1,x2,…,xn)⊤x=(x_1,x_2,…,x_n) , the output of the k-th hidden node (k=1,2,…,Lk=1,2,…,L) is defined as hk()=g(k⊤+bk),h_k(x)=g(w_k x+b_k), (1) where g(⋅)g(·) denotes an activation function, k∈ℝnw_k ^n are the hidden weights, and bk∈ℝb_k is the corresponding hidden bias. Let k∈ℝm β_k ^m denote the output weights. The global state approximation of the PI-SCM with L hidden nodes, denoted as L()u_L(x), is computed as: L()=∑k=1Lkhk().u_L(x)= _k=1^L β_kh_k(x). (2) In standard PINNs, the parameters k,bk,k\w_k,b_k, β_k\ are updated simultaneously via gradient descent. In contrast, the PI-SCM framework randomly assigns the hidden parameters kw_k and bkb_k within a problem-dependent adaptive range [−τ,τ][-τ,τ]. These parameters are fixed once assigned. The optimal output weights k β_k are determined analytically via pseudo-inverse solvers. This architecture serves as the foundation for solving the problems defined below. I-B Forward Problem We formulate the proposed framework for a general nonlinear differential system of order q≥1q≥ 1. Consider an open bounded domain Ω⊂ℝn ^n governing the underlying physical mechanisms, and let ∂Ω∂ denote its boundary. The physical state of the system is characterized by an unknown vector-valued function :Ω→ℝmu: ^m. To represent its high-order differential state, we introduce a multi-index =(γ1,γ2,…,γn)∈ℕ0n γ=( _1, _2,…, _n) _0^n, with ||=∑j=1nγj| γ|= _j=1^n _j, and define: D()=∂||()∂x1γ1∂x2γ2⋯∂xnγn,D()=().D γu(x)= ∂^| γ|u(x)∂ x_1 _1∂ x_2 _2·s∂ x_n _n, D^0u(x)=u(x). (3) For any nonnegative integer s, let ℐs=∈ℕ0n:||≤sI_s=\ γ _0^n:| γ|≤ s\ denote the set of all derivative indices up to order s. The corresponding generalized differential state is defined as (s)():=D()∈ℐs.D^(s)u(x):= \D γu(x) \_ γ _s. (4) Consequently, the generalized mapping of the nonlinear high-order physical system can be formulated as (,(q)())=,∀∈Ω,F\! (x,D^(q)u(x) )=0, ∈ , (5) subject to the initial/boundary conditions ℬ(,(qb)())=,∀∈∂Ω,B\! (x,D^(q_b)u(x) )=0, ∈∂ , (6) where 0≤qb≤q0≤ q_b≤ q, F is a nonlinear differential operator, and ℬB represents the specific boundary or initial constraint operators. Evaluating these continuous residuals computationally requires discretization. Let Ωp⊂Ω _p⊂ denote the finite set of interior physical collocation points, Ωb⊂∂Ω _b⊂∂ denote the boundary/initial training points, and Ωd⊂Ω _d⊂ denote the set of observational data points (in semi-supervised scenarios). Consequently, the objective of the forward problem is to determine the optimal set of network parameters k,bk,kk=1L\w_k,b_k, β_k\_k=1^L such that the state approximation L()u_L(x) simultaneously satisfies the governing physical equations (Eq. (5)) over Ωp _p, the initial/boundary conditions (Eq. (6)) over Ωb _b, and the empirical measurements over Ωd _d. I-C Inverse Problem In many practical scientific and engineering scenarios, the governing differential equations are characterized by unknown physical parameters that need to be inferred from empirical measurements. This constitutes the inverse problem. To formulate the inverse physical system, we modify the nonlinear high-order differential operator to explicitly incorporate an unknown parameter vector ∈ℝnp λ ^n_p. The governing equation is extended as: (,(q)(),)=,∀∈Ω.F\! (x,D^(q)u(x); λ )=0, ∈ . (7) Accordingly, the inverse problem tasks the network with simultaneously determining the optimal network parameters k,bk,kk=1L\w_k,b_k, β_k\_k=1^L and identifying the underlying physical parameter λ. The optimal joint configuration must ensure that the state approximation satisfies the parameterized physical laws governing the system while accurately fitting the empirical measurements. I Physics-Informed Stochastic Configuration Machine I-A Forward Problem The fundamental challenge in optimizing standard PINNs lies in the highly non-convex loss landscape induced by nonlinear differential operators. PI-SCM circumvents this bottleneck by decoupling the nonlinearity through local Taylor expansions, transforming the optimization into a sequence of linear least squares problems. We begin with the base localized optimization algorithm, denoted as PI-SC-I. Suppose the network has generated L−1L-1 hidden nodes, and the current state approximation is given by L−1()=∑k=1L−1khk()u_L-1(x)= _k=1^L-1 β_kh_k(x), where k∈ℝm β_k ^m denotes the output weight vector, hk()=g(k⊤+bk)h_k(x)=g(w_k x+b_k) represents the output of the k-th hidden node, g(⋅)g(·) is the activation function, and k,bkw_k,b_k are the randomly configured input weights and bias, respectively. When evaluating the addition of the L-th node hL()h_L(x), the updated state becomes L()=L−1()+LhL()u_L(x)=u_L-1(x)+ β_Lh_L(x). Evaluating at the unified collocation points denoted by x, let us define the physical residual at the L-th step as: Lp:=(,(q)L).e_L^p:=F\! (x,D^(q)u_L ). (8) Similarly, the corresponding initial/boundary condition residual evaluated at x is defined as: Lb:=ℬ(,(qb)L).e_L^b:=B\! (x,D^(q_b)u_L ). (9) By examining the residual difference between step L and step L−1L-1, we obtain: Lp−L−1p=(,(q)(L−1+LhL))−(,(q)L−1),Lb−L−1b=ℬ(,(qb)(L−1+LhL))−ℬ(,(qb)L−1). splite_L^p-e_L-1^p&=F\! (x,D^(q)(u_L-1+ β_Lh_L) )\\ & -F\! (x,D^(q)u_L-1 ),\\ e_L^b-e_L-1^b&=B\! (x,D^(q_b)(u_L-1+ β_Lh_L) )\\ & -B\! (x,D^(q_b)u_L-1 ). split (10) It becomes evident that if the underlying PDE system F is nonlinear with respect to any component of the generalized differential state (q)D^(q)u, the residual difference in the above equation is inherently nonlinear with respect to the newly added output weight L β_L. This nonlinearity directly precludes the use of standard linear least squares algorithms for explicitly determining L β_L. To decouple this nonlinearity and project it into a linear subspace, we apply the functional first-order Taylor expansion around the previously established generalized state (q)L−1D^(q)u_L-1. By analytically evaluating the local Jacobians of the differential operator with respect to every derivative component up to order q, the physical residual is linearized as Lp=(,(q)L)≈(,(q)L−1)+∑∈ℐqD(L−1)(DL−DL−1), splite_L^p&=F\! (x,D^(q)u_L ) \! (x,D^(q)u_L-1 )\\ &+ _ γ _q _D γu(u_L-1) (D γu_L-D γu_L-1 ), split (11) where the local Jacobian matrices evaluated at the complete previous differential state are mathematically denoted as: D(L−1)=∂(D)|(q)=(q)L−1,∈ℐq,J_D γu(u_L-1)= ∂(D γu) |_D^(q)u=D^(q)u_L-1, γ _q, (12) each yielding a matrix in ℝm×mR^m× m. In particular, the zero-order case = γ=0 gives D=J_D^0u=J_u. Because L β_L is independent of x, the differential increment associated with every multi-index satisfies: DL−DL−1=LDhL(),∈ℐq.D γu_L-D γu_L-1= β_LD γh_L(x), γ _q. (13) Substituting (13) into the linearized expansion yields Lp≈L−1p+L()L,e_L^p _L-1^p+A_L(x) β_L, (14) where L()=∑∈ℐqD(L−1)DhL().A_L(x)= _ γ _qJ_D γu(u_L-1)D γh_L(x). (15) The highly nonlinear high-order local residual evaluation is thus reduced to a purely linear mapping with respect to L β_L. Regarding the initial/boundary conditions, we assume ℬB can be formulated as a generalized high-order boundary operator: ℬ(,(qb))=∑∈ℐqb()D−bc()=,B\! (x,D^(q_b)u )= _ γ _q_bP_ γ(x)D γu-g_bc(x)=0, (16) where ()P_ γ(x) denotes the coefficient matrix associated with the derivative DD γu, and bc()g_bc(x) represents the prescribed boundary data. For qb=1q_b=1, selecting the derivative coefficients along the outward normal direction yields the standard Robin condition. Taking the difference between the boundary residuals of step L and step L−1L-1, we have: Lb−L−1b=(∑∈ℐqb()DhL())L.e_L^b-e_L-1^b= ( _ γ _q_bP_ γ(x)D γh_L(x) ) β_L. (17) It is evident that this formulation is strictly linear with respect to L β_L, which yields Lb=L−1b+L()L,e_L^b=e_L-1^b+C_L(x) β_L, (18) where L()=∑∈ℐqb()DhL()C_L(x)= _ γ _q_bP_ γ(x)D γh_L(x). Furthermore, to empower the framework with semi-supervised learning capabilities by integrating sparse observational data, the data residual evaluated at the points denoted by x is defined as Ld:=data()−L(),e_L^d:=u_data(x)-u_L(x), (19) where data()u_data(x) denotes the observations. The difference between the data residuals at step L and step L−1L-1 is computed as: Ld−L−1d=data()−L−1()−LhL()−(data()−L−1())=−hL()mL. splite_L^d-e_L-1^d&=u_data(x)-u_L-1(x)- β_Lh_L(x)\\ & - (u_data(x)-u_L-1(x) )\\ &=-h_L(x)I_m β_L. split (20) Thus, we have: Ld=L−1d+L()L,e_L^d=e_L-1^d+D_L(x) β_L, (21) where L()=−hL()mD_L(x)=-h_L(x)I_m. Assuming that at the L-th incremental step, we aim to simultaneously minimize the physical residual Lpe_L^p, the boundary residual Lbe_L^b, and the data residual Lde_L^d, the comprehensive optimization objective for the PI-SC-I algorithm is formulated as the following linear least squares loss function: ℒ(L) ( β_L) =1Np‖L()L+L−1p‖Ωp2 = 1N_p \|A_L(x) β_L+e_L-1^p \|_ _p^2 (22) +ωbcNb‖L()L+L−1b‖Ωb2 + _bcN_b \|C_L(x) β_L+e_L-1^b \|_ _b^2 +ωdataNd‖L()L+L−1d‖Ωd2, + _dataN_d \|D_L(x) β_L+e_L-1^d \|_ _d^2, where ωbc _bc and ωdata _data are the penalty parameters balancing the initial/boundary conditions and the observational data against the physical domain. The notations ∥⋅∥Ωp2\|·\|_ _p^2, ∥⋅∥Ωb2\|·\|_ _b^2, and ∥⋅∥Ωd2\|·\|_ _d^2 represent the standard discrete squared L2L_2 norms evaluated over Ωp=ipi=1Np _p=\x_i^p\_i=1^N_p, Ωb=ibi=1Nb _b=\x_i^b\_i=1^N_b, and Ωd=idi=1Nd _d=\x_i^d\_i=1^N_d, respectively, where NpN_p, NbN_b, and NdN_d are the corresponding numbers of physical, boundary, and data points. At the L-th step, the weighted generalized residual vector is defined as ~L=[1Npcoli=1,…,Np(Lp(ip))ωbcNbcoli=1,…,Nb(Lb(ib))ωdataNdcoli=1,…,Nd(Ld(id))]. e_L= bmatrix 1 N_p i=1,…,N_pcol (e_L^p(x_i^p) )\\ _bcN_b i=1,…,N_bcol (e_L^b(x_i^b) )\\ _dataN_d i=1,…,N_dcol (e_L^d(x_i^d) ) bmatrix. (23) where col(⋅)col(·) denotes the vertical concatenation operator. The three residual blocks in (22) can be combined into a unified linear least squares system. Let us define the weighted design matrix LM_L and the corresponding target vector Lv_L as L _L =[1Npcoli=1,…,Np(L(ip))ωbcNbcoli=1,…,Nb(L(ib))ωdataNdcoli=1,…,Nd(L(id))], = bmatrix 1 N_p i=1,…,N_pcol (A_L(x_i^p) )\\ _bcN_b i=1,…,N_bcol (C_L(x_i^b) )\\ _dataN_d i=1,…,N_dcol (D_L(x_i^d) ) bmatrix, L _L =−~L−1. =- e_L-1. (24) Accordingly, the objective in (22) can be written as ℒ(L)=‖LL−L‖22.L( β_L)= \|M_L β_L-v_L \|_2^2. (25) Consequently, the minimum-norm solution of the linearized least squares problem is explicitly computed as ^L=L†L β_L=M_L v_L. To control the Taylor truncation error of the true nonlinear high-order residual, the adopted output-weight increment is defined as L∗=αL†L, β_L^*= _L v_L, (26) where † denotes the Moore-Penrose pseudo-inverse operator, and α∈(0,1]α∈(0,1] acts as a damping factor. When α=1α=1, (26) coincides with the optimal linearized solution; when α<1α<1, it represents a damped step along the same analytical direction. As will be rigorously established in the subsequent theoretical analysis, introducing this parameter guarantees the existence of a feasible parameter configuration that enforces strict monotonic descent of the true residual. Furthermore, to ensure numerical stability against potential ill-conditioning of the localized Jacobian evaluations, we compute the pseudo-inverse via the Truncated Singular Value Decomposition (TSVD). By dynamically filtering out unstable singular components below a relative condition threshold η (e.g., η=10−10η=10^-10), this approach effectively circumvents the numerical instability induced by matrix ill-conditioning. To rigorously guarantee the effectiveness of the configured nonlinear parameters and the asymptotic convergence of the algorithm, we introduce a supervisory node selection mechanism. Using the weighted residual vector defined in (23), the comprehensive residual norm at the L-th step is given by ℰL2=‖~L‖2E_L^2=\| e_L\|^2. In the PI-SC-I algorithm, we establish a contraction evaluation metric ξL _L: ξL=(r+μL)ℰL−12−ℰL2, _L=(r+ _L)E_L-1^2-E_L^2, (27) where r∈(0,1)r∈(0,1), and μL=1−rL+1 _L= 1-rL+1. During the generation of the L-th hidden node, we perform TmaxT_max independent trials by randomly sampling the hidden parameters (L,bLw_L,b_L) and calculating their corresponding optimal output weights L∗ β_L^*. Full candidate configurations satisfying ξL>0 _L>0 are deposited into a feasible pool. The algorithm ultimately selects the optimal parameter set that maximizes ξL _L, thereby completing the algorithm of PI-SC-I and enforcing the steepest residual descent. Remark 1 It is imperative to emphasize that the physical residual Lpe_L^p utilized in computing ℰL2E_L^2 represents the true physical residual, rather than its linearized Taylor approximation. As will be rigorously proven in Theorem 1, adhering to this true residual is the fundamental prerequisite for guaranteeing the universal approximation property of the PI-SC-I framework. In the original SCN framework [13], a global evaluation of the output weights was designed after the local node selection to overcome the slow convergence caused by the strictly local constructive scheme. Here, we adopt an analogous global updating mechanism for the proposed framework, establishing the PI-SC-I algorithm. Let L()=[1()m,2()m,…,L()m]H_L(x)=[h_1(x)I_m,h_2(x)I_m,…,h_L(x)I_m] denote the global hidden layer output matrix, and the global output weight vector is denoted as glob=[1⊤,2⊤,…,L⊤]⊤ β_glob=[ β_1 , β_2 ,…, β_L ] . The globally updated state approximation is denoted as L()=L()globu_L(x)=H_L(x) β_glob. The physical residual at the L-th step is thus given by: Lp=(,(q)(L()glob)).e_L^p=F\! (x,D^(q)(H_L(x) β_glob) ). (28) Similarly, considering the nonlinearity of the differential operator F with respect to glob β_glob, we apply the first-order Taylor expansion around the locally constructed generalized baseline state (q)locD^(q)u_loc. Specifically, let loc∗=[1⊤,2⊤,…,(L∗)⊤]⊤ β_loc^*=[ β_1 , β_2 ,…,( β_L^*) ] denote the configuration achieved by the PI-SC-I algorithm at the current step, yielding the local state approximation loc()=L()loc∗u_loc(x)=H_L(x) β_loc^* and its corresponding physical residual locpe_loc^p. Recalling the local Jacobian matrices defined in Eq. (12), we evaluate them directly at the current localized configuration loc()u_loc(x). The physical residual is thereby linearized via these Jacobian operators as Lp≈locp+∑∈ℐqD(loc)(DL()glob−Dloc())=glob()glob−(), splite_L^p& _loc^p+ _ γ _qJ_D γu(u_loc) (D γH_L(x) β_glob-D γu_loc(x) )\\ &=A_glob(x) β_glob-b(x), split (29) where glob() _glob(x) =∑∈ℐqD(loc)DL(), = _ γ _qJ_D γu(u_loc)D γH_L(x), (30) () (x) =∑∈ℐqD(loc)Dloc()−locp. = _ γ _qJ_D γu(u_loc)D γu_loc(x)-e_loc^p. Following the same rationale, the initial/boundary condition residual is globally formulated as Lb=∑∈ℐqb()DL()glob−bc()=glob()glob−bc(), splite_L^b&= _ γ _q_bP_ γ(x)D γH_L(x) β_glob-g_bc(x)\\ &=C_glob(x) β_glob-g_bc(x), split (31) where glob()=∑∈ℐqb()DL()C_glob(x)= _ γ _q_bP_ γ(x)D γH_L(x). The data residual is globally formulated as: Ld=data()−L()glob=data()−glob()glob, splite_L^d&=u_data(x)-H_L(x) β_glob\\ &=u_data(x)-D_glob(x) β_glob, split (32) where glob()=L()D_glob(x)=H_L(x). Building upon these global algebraic formulations, the comprehensive objective function for the PI-SC-I algorithm is constructed as a linear least squares problem: ℒglob(glob)=1Np‖glob()glob−()‖Ωp2+ωbcNb‖glob()glob−bc()‖Ωb2+ωdataNd‖glob()glob−data()‖Ωd2. splitL_glob( β_glob)&= 1N_p \|A_glob(x) β_glob-b(x) \|_ _p^2\\ & + _bcN_b \|C_glob(x) β_glob-g_bc(x) \|_ _b^2\\ & + _dataN_d \|D_glob(x) β_glob-u_data(x) \|_ _d^2. split (33) The three terms in (33) can be combined into a unified linear least squares system. The corresponding weighted global design matrix globM_glob and target vector globv_glob are given below: glob=[1Npcoli=1,…,Np(glob(ip))ωbcNbcoli=1,…,Nb(glob(ib))ωdataNdcoli=1,…,Nd(glob(id))].M_glob= bmatrix 1 N_p i=1,…,N_pcol (A_glob(x_i^p) )\\ _bcN_b i=1,…,N_bcol (C_glob(x_i^b) )\\ _dataN_d i=1,…,N_dcol (D_glob(x_i^d) ) bmatrix. (34a) glob=[1Npcoli=1,…,Np((ip))ωbcNbcoli=1,…,Nb(bc(ib))ωdataNdcoli=1,…,Nd(data(id))].v_glob= bmatrix 1 N_p i=1,…,N_pcol (b(x_i^p) )\\ _bcN_b i=1,…,N_bcol (g_bc(x_i^b) )\\ _dataN_d i=1,…,N_dcol (u_data(x_i^d) ) bmatrix. (34b) Accordingly, the objective in (33) can be written as ℒglob(glob)=‖globglob−glob‖22.L_glob( β_glob)= \|M_glob β_glob-v_glob \|_2^2. (35) In practice, as the hidden layer incrementally expands, the global design matrix globM_glob becomes increasingly susceptible to severe ill-conditioning. Consequently, the output weights for the entire PI-SC-I network are determined by extending the aforementioned TSVD strategy to evaluate the global pseudo-inverse: ~glob∗=glob†glob. β_glob^*=M_glob v_glob. (36) This unified pseudo-inverse approach thoroughly recalculates the entire hidden layer’s output configuration. However, because ~glob∗ β_glob^* is analytically optimal exclusively within the linearized subspace, the nonlinear nature of the underlying PDE operator might occasionally induce a severe Taylor truncation penalty. Applying this full step blindly could lead to an unexpected surge in the true nonlinear residual. To guarantee the monotonic convergence of the network within PI-SC-I, we introduce a discrete line search mechanism. We define the global search direction as Δ=~glob∗−loc∗ β= β_glob^*- β_loc^*, yielding a parameterized weight trajectory for a scaling factor α∈[0,1]α∈[0,1]: (α)=loc∗+αΔ. β(α)= β_loc^*+α β. (37) Let ℰL2(α)E_L^2(α) denote the comprehensive true residual norm evaluated with the candidate weights (α) β(α). The algorithm computes in parallel a discrete set of uniformly distributed candidate step sizes, denoted as ⊂[0,1]A⊂[0,1], which includes the boundary limits 1.01.0 and 0.00.0. The best scaling factor α∗α^* within A is selected as α∗=argminα∈ℰL2(α).α^*= _α E_L^2(α). (38) The configuration corresponding to this α∗α^* is then adopted as the final updated weight vector ∗=(α∗) β^*= β(α^*). Although the global updating mechanism of PI-SC-I achieves higher accuracy, the computational complexity grows rapidly with the number of hidden nodes. Consequently, for large-scale data analysis or complex physical simulations requiring dense node configurations, this exhaustive recalculation incurs substantial computational overhead. To strike a balance between approximation accuracy and computational efficiency, a sliding-window strategy, as originally proposed in the SCN literature [13], can be integrated into the framework to formulate the PI-SC-I algorithm. Specifically, given a window size W<LW<L, the update in PI-SC-I is restricted to the most recent W nodes while freezing the preceding L−WL-W nodes. Following the conventional SC-I partitioning scheme, we decompose the hidden layer matrix and weights as L()=[pre(),win()]H_L(x)=[H_pre(x),H_win(x)] and glob=[pre⊤,win⊤]⊤ β_glob=[ β_pre , β_win ] , respectively. The frozen state is thus fixed as L−W()=pre()preu_L-W(x)=H_pre(x) β_pre. To extend this block-wise projection principle to our physics-informed setting, we first systematically analyze the remaining discrepancies across the computational domains. Instead of re-evaluating the entire network, the active window aims to fit the gap between the global targets (b, bcg_bc, datau_data) and the contributions already provided by the frozen state L−Wu_L-W. These residual target vectors are explicitly evaluated over their respective domains as: L−Wp() _L-W^p(x) =()−∑∈ℐqD(loc)DL−W(), =b(x)- _ γ _qJ_D γu(u_loc)D γu_L-W(x), (39) L−Wb() _L-W^b(x) =bc()−∑∈ℐqb()DL−W(), =g_bc(x)- _ γ _q_bP_ γ(x)D γu_L-W(x), L−Wd() _L-W^d(x) =data()−L−W(). =u_data(x)-u_L-W(x). Concurrently, applying the linearized operators exclusively to the windowed hidden nodes winH_win yields: W()=∑∈ℐqD(loc)Dwin(),W()=∑∈ℐqb()Dwin(),W()=win(). splitA_W(x)&= _ γ _qJ_D γu(u_loc)D γH_win(x),\\ C_W(x)&= _ γ _q_bP_ γ(x)D γH_win(x),\\ D_W(x)&=H_win(x). split (40) The comprehensive objective function is constructed as a unified least squares problem: ℒwin(win)=1Np‖W()win−L−Wp()‖Ωp2+ωbcNb‖W()win−L−Wb()‖Ωb2+ωdataNd‖W()win−L−Wd()‖Ωd2. splitL_win( β_win)&= 1N_p \|A_W(x) β_win-r_L-W^p(x) \|_ _p^2\\ & + _bcN_b \|C_W(x) β_win-r_L-W^b(x) \|_ _b^2\\ & + _dataN_d \|D_W(x) β_win-r_L-W^d(x) \|_ _d^2. split (41) The three terms in (41) can be combined into a unified linear least squares system. The corresponding weighted windowed design matrix WM_W and target vector Wv_W are given below: W=[1Npcoli=1,…,Np(W(ip))ωbcNbcoli=1,…,Nb(W(ib))ωdataNdcoli=1,…,Nd(W(id))].M_W= bmatrix 1 N_p i=1,…,N_pcol (A_W(x_i^p) )\\ _bcN_b i=1,…,N_bcol (C_W(x_i^b) )\\ _dataN_d i=1,…,N_dcol (D_W(x_i^d) ) bmatrix. (42a) W=[1Npcoli=1,…,Np(L−Wp(ip))ωbcNbcoli=1,…,Nb(L−Wb(ib))ωdataNdcoli=1,…,Nd(L−Wd(id))].v_W= bmatrix 1 N_p i=1,…,N_pcol (r_L-W^p(x_i^p) )\\ _bcN_b i=1,…,N_bcol (r_L-W^b(x_i^b) )\\ _dataN_d i=1,…,N_dcol (r_L-W^d(x_i^d) ) bmatrix. (42b) Accordingly, the objective in (41) can be written as ℒwin(win)=‖Wwin−W‖22.L_win( β_win)= \|M_W β_win-v_W \|_2^2. (43) Consistent with our strategy to suppress potential ill-conditioning, the output weights for the current window are evaluated via the TSVD pseudo-inverse: win∗=W†W. β_win^*=M_W v_W. (44) Finally, the complete window-updated weight vector is explicitly reassembled as: ~win∗=[prewin∗]. β_win^*= bmatrix β_pre\\ β_win^* bmatrix. (45) This sliding-window scheme bounds the number of columns of the design matrix WM_W to mWmW, regardless of the total network width L. By dynamically replacing the full-batch recalculation with this windowed block-wise pseudo-inverse, the PI-SC-I framework reduces the computational burden while still reliably preserving the accelerated convergence benefits. However, similar to the global approach in PI-SC-I, because ~win∗ β_win^* is analytically optimal exclusively within the linearized subspace, applying this windowed step blindly could induce a severe Taylor truncation penalty. To guarantee the monotonic convergence of the network within PI-SC-I, we introduce the same discrete line search mechanism. The search direction and its associated weight trajectory are defined as Δ=~win∗−loc∗,(α)=loc∗+αΔ. β= β_win^*- β_loc^*, β(α)= β_loc^*+α β. (46) Once this windowed direction is constructed, the step size α∗α^* and the final updated weight vector ∗=(α∗) β^*= β(α^*) are determined using the same method described for PI-SC-I. I-B Inverse Problem To extend the PI-SCM framework to inverse problems, we define the governing equation as (,(q),)=F(x,D^(q)u; λ)=0. Within the localized construction phase (PI-SC-I), the incremental updating objective is formulated to jointly estimate the local output weight L β_L and the physical parameter λ. To achieve this joint evaluation, we build upon the high-order linearization strategy established in Section I-A. By performing a simultaneous first-order Taylor expansion with respect to both (q)D^(q)u and λ around the previous configuration ((q)L−1,L−1)(D^(q)u_L-1, λ_L-1), the linearized physical equation governing the current step is explicitly derived as L()L+λ(,L−1,L−1)(L−L−1)≈−L−1p,A_L(x) β_L+J_λ(x;u_L-1, λ_L-1)( λ_L- λ_L-1)≈-e_L-1^p, (47) where L−1pe_L-1^p and LA_L are defined in (14) and (15), and the parameter Jacobian evaluated at the previous generalized state is defined as λ(,L−1,L−1)=∂|(q)=(q)L−1,=L−1J_λ(x;u_L-1, λ_L-1)= ∂ λ |_D^(q)u=D^(q)u_L-1, λ= λ_L-1. To evaluate the parameter increment directly, we define ΔL=L−L−1 λ_L= λ_L- λ_L-1. This yields the linearized unified system: L()L+λ(,L−1,L−1)ΔL≈−L−1p.A_L(x) β_L+J_λ(x;u_L-1, λ_L-1) λ_L≈-e_L-1^p. (48) We define the localized joint update vector as L=[L⊤,ΔL⊤]⊤ θ_L=[ β_L , λ_L ] . Casting the equation into a unified matrix-vector product yields the augmented design matrix L,aug()=[L(),λ(,L−1,L−1)]A_L,aug(x)=[A_L(x),J_λ(x;u_L-1, λ_L-1)]. Since the boundary constraints and observational data are inherently independent of the physical parameters, their forward design matrices (LC_L and LD_L, inherited from (18) and (21)) are systematically zero-padded, yielding L,aug()=[L(),]C_L,aug(x)=[C_L(x),0] and L,aug()=[L(),]D_L,aug(x)=[D_L(x),0]. The joint parameter identification is rigorously formulated as a linear least-squares minimization problem: ℒ(L) ( θ_L) =1Np‖L,aug()L+L−1p‖Ωp2 = 1N_p\|A_L,aug(x) θ_L+e_L-1^p\|_ _p^2 (49) +ωbcNb‖L,aug()L+L−1b‖Ωb2 + _bcN_b\|C_L,aug(x) θ_L+e_L-1^b\|_ _b^2 +ωdataNd‖L,aug()L+L−1d‖Ωd2. + _dataN_d\|D_L,aug(x) θ_L+e_L-1^d\|_ _d^2. Solving this linear least squares problem via TSVD yields the candidate joint solution ^L θ_L. Consistent with the forward PI-SC-I algorithm, the adopted configuration is damped as L∗=α^L θ_L^*=α θ_L, from which L∗ β_L^* and ΔL∗ λ_L^* are extracted. Subsequently, these candidate configurations are rigorously filtered through the established supervisory node selection mechanism to guarantee monotonic residual descent. Upon successful addition of the L-th node, the absolute local baseline is established as (loc,loc)(u_loc, λ_loc), where loc=L−1+ΔL∗ λ_loc= λ_L-1+ λ_L^*. Updating λ strictly based on a localized node addition restricts identification accuracy. To overcome this limitation, we extend the established PI-SC-I algorithm for the precise joint identification of global output weights and parameters. Specifically, we perform the Taylor expansion around the newly established intermediate baseline ((q)loc,loc)(D^(q)u_loc, λ_loc). Let Δglob=−loc λ_glob= λ- λ_loc denote the parameter increment. The global augmented joint vector is defined as glob=[glob⊤,Δglob⊤]⊤ θ_glob=[ β_glob , λ_glob ] . Building upon the algebraic formulations derived in the forward PI-SC-I algorithm, the expanded global linearized physical equation becomes: glob()glob+λ(,loc,loc)Δglob≈(),A_glob(x) β_glob+J_λ(x;u_loc, λ_loc) λ_glob (x), (50) where globA_glob and b are defined in (29) and (30), and λ(,loc,loc)=∂|(q)=(q)loc,=locJ_λ(x;u_loc, λ_loc)= ∂ λ |_D^(q)u=D^(q)u_loc, λ= λ_loc. Let λ,loc()=λ(,loc,loc)J_λ,loc(x)=J_λ(x;u_loc, λ_loc). We define the global augmented matrices as glob,aug()=[glob(),λ,loc()]A_glob,aug(x)=[A_glob(x),J_λ,loc(x)], glob,aug()=[glob(),]C_glob,aug(x)=[C_glob(x),0], and glob,aug()=[glob(),]D_glob,aug(x)=[D_glob(x),0]. The joint update is formulated as the following linear least-squares problem: ℒglob(glob) _glob( θ_glob) =1Np‖glob,aug()glob−()‖Ωp2 = 1N_p\|A_glob,aug(x) θ_glob-b(x)\|_ _p^2 (51) +ωbcNb‖glob,aug()glob−bc()‖Ωb2 + _bcN_b\|C_glob,aug(x) θ_glob-g_bc(x)\|_ _b^2 +ωdataNd‖glob,aug()glob−data()‖Ωd2. + _dataN_d\|D_glob,aug(x) θ_glob-u_data(x)\|_ _d^2. Solving this linear least squares problem via TSVD yields the candidate joint solution ~glob∗=[(~glob∗)⊤,(Δ~glob∗)⊤]⊤ θ_glob^*=[( β_glob^*) ,( λ_glob^*) ] . To alleviate the computational bottleneck of this full-batch joint inversion, the PI-SC-I algorithm is similarly adapted. Given a predefined window size W, we define the windowed parameter increment as Δwin=−loc λ_win= λ- λ_loc and construct win=[win⊤,Δwin⊤]⊤ θ_win=[ β_win , λ_win ] . By extracting WA_W, WC_W, and WD_W defined in (40), we construct the windowed augmented matrices as W,aug()=[W(),λ,loc()]A_W,aug(x)=[A_W(x),J_λ,loc(x)], W,aug()=[W(),]C_W,aug(x)=[C_W(x),0], and W,aug()=[W(),]D_W,aug(x)=[D_W(x),0]. The joint update is formulated as the following linear least-squares problem: ℒwin(win) _win( θ_win) =1Np‖W,aug()win−L−Wp()‖Ωp2 = 1N_p\|A_W,aug(x) θ_win-r_L-W^p(x)\|_ _p^2 (52) +ωbcNb‖W,aug()win−L−Wb()‖Ωb2 + _bcN_b\|C_W,aug(x) θ_win-r_L-W^b(x)\|_ _b^2 +ωdataNd‖W,aug()win−L−Wd()‖Ωd2. + _dataN_d\|D_W,aug(x) θ_win-r_L-W^d(x)\|_ _d^2. Solving this linear least squares problem via TSVD yields the candidate joint solution win∗=[(win∗)⊤,(Δ~win∗)⊤]⊤ θ_win^*=[( β_win^*) ,( λ_win^*) ] . To guarantee the monotonic decrease of the residual sequence during this joint parameter-weight update, we extend the discrete line search mechanism to the inverse problem. We define the absolute joint configuration vector as =[⊤,⊤]⊤ =[ β , λ ] . Let loc∗=[loc∗⊤,loc⊤]⊤ _loc^*=[ β_loc^* , λ_loc ] denote the absolute baseline joint configuration secured by PI-SC-I. Depending on the utilized algorithm, ~∗ ^* is explicitly reassembled as: ~∗=[~glob∗loc+Δ~glob∗],for PI-SC-I[prewin∗loc+Δ~win∗],for PI-SC-I. ^*= cases bmatrix β_glob^*\\ λ_loc+ λ_glob^* bmatrix,&for PI-SC-I\\[15.00002pt] bmatrix β_pre\\ β_win^*\\ λ_loc+ λ_win^* bmatrix,&for PI-SC-I cases. (53) We then establish a linear search path governed by a scaling factor α∈[0,1]α∈[0,1] as: (α)=loc∗+α(~∗−loc∗). (α)= _loc^*+α( ^*- _loc^*). (54) By computing in parallel the comprehensive true residual norm ℰL2(α)E_L^2(α) across the discrete candidate set A, the optimal scaling factor α∗α^* is determined via explicit minimization: α∗=argminα∈ℰL2(α).α^*= _α E_L^2(α). (55) The optimally scaled absolute configuration is then adopted as the final updated vector ∗=(α∗) ^*= (α^*). To synthesize the mathematical formulations derived for both forward and inverse learning tasks, the comprehensive training strategies of the PI-SCM framework are summarized in the following algorithms. Algorithm 1 details the fundamental localized construction process (PI-SC-I), while Algorithm 2 and Algorithm 3 describe the advanced sliding-window (PI-SC-I) and global updating (PI-SC-I) mechanisms, respectively. Algorithm 1 Training strategy of PI-SC-I 0: Given coordinates =(x1,…,xn)⊤x=(x_1,…,x_n) partitioned into collocation sets Ωp,Ωb,Ωd _p, _b, _d; Sparse observational target datau_data; The differential operator F and initial/boundary operator ℬB; Initialized Lmax,TmaxL_max,T_max; Error tolerance ϵε; Three sets of scalars Υ=τ1,…,τend,ℛ=r1,…,rend,1=α1,1:Δα1:α1,end =\ _1,…, _end\,R=\r_1,…,r_end\,A_1=\ _1,1: _1: _1,end\; Penalty weights ωbc,ωdata _bc, _data. 0: Hidden parameters L∗,bL∗L=1Lend\w_L^*,b_L^*\_L=1^L_end, output weights L∗L=1Lend\ β_L^*\_L=1^L_end (and identified physical parameter ∗ λ^* for inverse tasks). 1: Initialize L=1,0=L=1,u_0=0. If inverse problem, initialize 0 λ_0. 2: Evaluate initial residuals 0p,0b,0de_0^p,e_0^b,e_0^d and comprehensive norm ℰ02E_0^2. 3: while L≤LmaxL≤ L_max and ℰL−12>ϵE^2_L-1>ε do 4: Initialize two empty sets P and Ω1 _1. 5: for α1∈1 _1 _1 do 6: for τ∈Υτ∈ do 7: for k=1,2,…,Tmaxk=1,2,…,T_max do 8: Randomly sample hidden parameters L∼[−τ,τ]nw_L [-τ,τ]^n and bL∼[−τ,τ]b_L [-τ,τ]. 9: For forward tasks, assemble L,LM_L,v_L by (24) and compute L∗ β_L^* by (26). For inverse tasks, solve (49) via TSVD to obtain ^L θ_L, and set L∗=α1^L θ_L^*= _1 θ_L. 10: Extract L∗ β_L^* (and Δ∗ λ^* for inverse tasks). 11: Update L=L−1+L∗hLu_L=u_L-1+ β_L^*h_L. 12: Evaluate ℰL2E_L^2. 13: Calculate ξL=(r+μL)ℰL−12−ℰL2 _L=(r+ _L)E_L-1^2-E_L^2, where μL=(1−r)/(L+1) _L=(1-r)/(L+1). 14: if ξL>0 _L>0 then 15: Save L,bL,L∗w_L,b_L, β_L^* (and Δ∗ λ^* for inverse tasks) in P, ξL _L in Ω1 _1. 16: end if 17: end for 18: if P is not empty then 19: Find L,bL,L∗w_L,b_L, β_L^* (and Δ∗ λ^* for inverse tasks) that maximize ξL _L in Ω1 _1. 20: Update Lu_L, Lp,Lb,Lde_L^p,e_L^b,e_L^d, ℰL2E_L^2 (and set L∗=L−1+ΔL∗ λ_L^*= λ_L-1+ λ_L^*). 21: Break (go to Step 27). 22: else 23: Update r to the subsequent value in ℛR. 24: Continue 25: end if 26: end for 27: end for 28: L←L+1L← L+1. 29: end while 30: return L∗,bL∗L=1Lend\w_L^*,b_L^*\_L=1^L_end, L∗L=1Lend\ β_L^*\_L=1^L_end (and ∗ λ^* for inverse tasks). Algorithm 2 Training strategy of PI-SC-I 0: The inputs of Algorithm 1; Window size W<LmaxW<L_max; 2=α2,1:Δα2:α2,endA_2=\ _2,1: _2: _2,end\. 0: Hidden parameters L∗,bL∗L=1Lend\w_L^*,b_L^*\_L=1^L_end, output weights L∗L=1Lend\ β_L^*\_L=1^L_end (and identified parameter ∗ λ^* for inverse tasks). 1: Initialize L=1,0=L=1,u_0=0. If inverse problem, initialize 0 λ_0. 2: Evaluate initial residuals 0p,0b,0de_0^p,e_0^b,e_0^d and comprehensive norm ℰ02E_0^2. 3: while L≤LmaxL≤ L_max and ℰL−12>ϵE^2_L-1>ε do 4: Execute Step 4-27 of Algorithm 1 to secure L∗,bL∗w_L^*,b_L^* and loc∗ β_loc^* (or loc∗ _loc^* for inverse tasks). 5: if L≤WL≤ W then 6: For forward tasks, assemble glob,globM_glob,v_glob by (34) and compute ~∗ β^* by (36). For inverse tasks, solve (51) via TSVD to obtain ~glob∗ θ_glob^*. 7: else 8: Retrieve pre=[1∗⊤,…,L−W∗⊤]⊤ β_pre=[ β_1^* ,…, β_L-W^* ] . 9: Calculate L−Wp,L−Wb,L−Wdr_L-W^p,r_L-W^b,r_L-W^d by (39). 10: For forward tasks, assemble W,WM_W,v_W by (42) and compute win∗ β_win^* by (44). 11: For inverse tasks, solve (52) via TSVD to obtain win∗ θ_win^*. 12: For forward tasks, let ~∗=[pre⊤,(win∗)⊤]⊤ β^*=[ β_pre ,( β_win^*) ] . 13: end if 14: Initialize an empty set Ω2 _2. 15: Obtain Δ=~∗−loc∗ β= β^*- β_loc^* (or obtain ~∗ ^* by (53), Δ=~∗−loc∗ = ^*- _loc^* for inverse tasks). 16: for α2∈2 _2 _2 do 17: Calculate (α2)=loc∗+α2Δ β( _2)= β_loc^*+ _2 β (or (α2) ( _2) for inverse tasks) and corresponding ℰL2(α2)E_L^2( _2). Save ℰL2(α2)E_L^2( _2) to Ω2 _2. 18: end for 19: Find α2∗α^*_2 that minimizes ℰL2(α2)E_L^2( _2) in Ω2 _2 and adopt ∗=(α2∗) β^*= β(α^*_2) (or ∗=(α2∗) ^*= (α^*_2) for inverse tasks). 20: Extract ∗ β^* (and ∗ λ^* for inverse tasks). 21: Update Lu_L, Lp,Lb,Lde_L^p,e_L^b,e_L^d, ℰL2E_L^2 (and set L∗=∗ λ_L^*= λ^*). 22: L←L+1L← L+1. 23: end while 24: return L∗,bL∗L=1Lend\w_L^*,b_L^*\_L=1^L_end, L∗L=1Lend\ β_L^*\_L=1^L_end (and ∗ λ^* for inverse tasks). Algorithm 3 Training strategy of PI-SC-I 0: The inputs of Algorithm 1; 2=α2,min:Δα2:α2,maxA_2=\ _2,min: _2: _2,max\. 0: Hidden parameters L∗,bL∗L=1Lend\w_L^*,b_L^*\_L=1^L_end, output weights L∗L=1Lend\ β_L^*\_L=1^L_end (and identified parameter ∗ λ^* for inverse tasks). 1: Initialize L=1,0=L=1,u_0=0. If inverse problem, initialize 0 λ_0. 2: Evaluate initial residuals 0p,0b,0de_0^p,e_0^b,e_0^d and comprehensive norm ℰ02E_0^2. 3: while L≤LmaxL≤ L_max and ℰL−12>ϵE^2_L-1>ε do 4: Execute Step 4-27 of Algorithm 1 to secure L∗,bL∗w_L^*,b_L^* and loc∗ β_loc^* (or loc∗ _loc^* for inverse tasks). 5: For forward tasks, assemble glob,globM_glob,v_glob by (34) and compute ~glob∗ β_glob^* by (36). For inverse tasks, solve (51) via TSVD to obtain ~glob∗ θ_glob^*. 6: Initialize an empty set Ω2 _2. 7: Obtain Δ=~glob∗−loc∗ β= β_glob^*- β_loc^* (or obtain ~∗ ^* by (53), Δ=~∗−loc∗ = ^*- _loc^* for inverse tasks). 8: for α2∈2 _2 _2 do 9: Calculate (α2)=loc∗+α2Δ β( _2)= β_loc^*+ _2 β (or (α2) ( _2) for inverse tasks) and corresponding ℰL2(α2)E_L^2( _2). Save ℰL2(α2)E_L^2( _2) to Ω2 _2. 10: end for 11: Find α2∗α^*_2 that minimizes ℰL2(α2)E_L^2( _2) in Ω2 _2 and adopt ∗=(α2∗) β^*= β(α^*_2) (or ∗=(α2∗) ^*= (α^*_2) for inverse tasks). 12: Extract ∗ β^* (and ∗ λ^* for inverse tasks). 13: Update Lu_L, Lp,Lb,Lde_L^p,e_L^b,e_L^d, ℰL2E_L^2 (and set L∗=∗ λ_L^*= λ^*). 14: L←L+1L← L+1. 15: end while 16: return L∗,bL∗L=1Lend\w_L^*,b_L^*\_L=1^L_end, L∗L=1Lend\ β_L^*\_L=1^L_end (and ∗ λ^* for inverse tasks). I-C Universal Approximation Property Building upon the algorithmic formulations established in Section I-A and Section I-B, this section constructs a unified theoretical foundation for the PI-SCM framework. We rigorously establish the universal approximation property of the proposed algorithms, proving that PI-SC-I, PI-SC-I, and PI-SC-I inherently guarantee the asymptotic convergence of the nonlinear residuals to zero. While the ensuing explicit proofs are formulated within the context of the forward problem, the residual-convergence guarantees extend to the inverse problem under the same augmented-feature conditions. As established in Section I-B, defining L θ_L maps the joint parameter identification into an identical linearized subspace. This algebraic structure preserves the underlying contraction argument, provided that the augmented Jacobian satisfies the required non-orthogonality condition. Unique recovery of λ additionally requires the conventional local identifiability condition on the parameter Jacobian. Before proceeding, we formally articulate a standard regularity assumption regarding the differential operator. Assumption 1 Let =col∈ℐq(D)=col((q)).z= γ _qcol (D γu )=col\! (D^(q)u ). (56) The activation function satisfies g∈Cq(ℝ)g∈ C^q(R). The differential operator F is assumed to be twice continuously Fréchet differentiable (C2C^2) with respect to z. Consequently, its second Fréchet derivative is locally bounded. As the convergence of entire PI-SCM architecture relies on the stability of incremental node additions, we initially prove the universal approximation property for the PI-SC-I algorithm, which serves as the convergent baseline for all subsequent schemes. To facilitate a concise derivation, we recall the weighted generalized residual vector ~L e_L, its comprehensive norm ℰL2E_L^2, the weighted design matrix LM_L, and the target vector L=−~L−1v_L=- e_L-1 established in Section I-A. By applying the first-order Taylor expansion to ~L e_L, we have ~L=~L−1+LL+vec(L)=LL−L+vec(L) e_L= e_L-1+M_L β_L+r_vec( β_L)=M_L β_L-v_L+r_vec( β_L). The Taylor truncation error originates exclusively from the nonlinear physical residual. The true nonlinear residual norm at the L-th step is thus given by: ℰL2=‖LL−L‖2+ℛTaylor(L),E_L^2=\|M_L β_L-v_L\|^2+R_Taylor( β_L), (57) where ℛTaylor(L)R_Taylor( β_L) is defined as: ℛTaylor(L)=2(LL−L)⊤vec(L)+‖vec(L)‖2.R_Taylor( β_L)=2(M_L β_L-v_L) r_vec( β_L)+\|r_vec( β_L)\|^2. (58) As shown in (58), ℛTaylor(L)R_Taylor( β_L) represents the omitted second-order Taylor remainder. Beyond operator regularity, convergence requires the admissible hidden-node family Γ to provide a nonorthogonal linearized direction for every nonzero discrete residual, as formalized below. Assumption 2 At each incremental step L, the admissible hidden-node family Γ is residual-complete in the sense that, for every nonzero vector ~ e in the discrete generalized residual space, there exists an h∈Γh∈ whose weighted design matrix L(h)M_L(h) satisfies L(h)⊤~≠.M_L(h) e 0. (59) Under these assumptions, we now establish the universal approximation property of the PI-SC-I algorithm. Theorem 1 Suppose that Assumptions 1 and 2 hold, and the generalized residuals are bounded. Given r∈(0,1)r∈(0,1) and a nonnegative real number sequence μL\ _L\ with limL→+∞μL=0 _L→+∞ _L=0 and μL≤(1−r) _L≤(1-r), let the threshold sequence be defined as δL=(1−r−μL)ℰL−12 _L=(1-r- _L)E_L-1^2. Then, there exists a candidate hidden node hL∈Γh_L∈ and a scaling factor α¯∈(0,1] α∈(0,1], such that for any α∈(0,α¯)α∈(0, α), the output weight vector, defined as L(α)=αL†L, β_L(α)= _L v_L, (60) guarantees that the following inequality holds: ℰL−12−ℰL2≥δL.E_L-1^2-E_L^2≥ _L. (61) By accepting the node and the output weights that satisfy (61), the sequence of ℰL2\E_L^2\ monotonically decreases, satisfying limL→+∞ℰL2=0 _L→+∞E_L^2=0. Proof: If ~L−1= e_L-1=0, the desired residual convergence has already been achieved. Otherwise, ~L−1≠ e_L-1 0. According to Assumption 2, there exists an admissible hidden node hL∈Γh_L∈ such that L⊤~L−1≠M_L e_L-1 0. Since L=−~L−1v_L=- e_L-1, it follows that L⊤L≠M_L v_L 0. Then, we consider the output weight vector L(α) β_L(α) defined in (60) and define L=LL†P_L=M_LM_L . It satisfies L⊤=LP_L =P_L and L2=LP_L^2=P_L. Hence, LL(α)−L=(αL−)LM_L β_L(α)-v_L=( _L-I)v_L. Substituting this relation into the true residual expansion, we define Δℰ2(α) ^2(α) as the reduction in the true residual norm: Δℰ2(α) ^2(α) =ℰL−12−ℰL2 =E_L-1^2-E_L^2 (62) =α(2−α)‖LL‖2−ℛTaylor(αL†L). =α(2-α)\|P_Lv_L\|^2-R_Taylor( _L v_L). Because L⊤L≠M_L v_L 0, LLP_Lv_L is nonzero, guaranteeing that ‖LL‖2>0\|P_Lv_L\|^2>0. According to Assumption 1, the remainder is bounded. Since the generalized residuals are assumed to be globally bounded, the absolute magnitude of the expanded scalar penalty, denoted as |ℛTaylor||R_Taylor|, inherently inherits this quadratic bound. Thus, for every accepted hidden node there exists a finite constant C~L>0 C_L>0 yielding |ℛTaylor|≤α2C~L‖L†L‖2|R_Taylor|≤α^2 C_L\|M_L v_L\|^2. We apply −ℛTaylor≥−|ℛTaylor|-R_Taylor≥-|R_Taylor| to (62), yielding: Δℰ2(α)≥α(2−α)‖LL‖2−α2C~L‖L†L‖2. ^2(α)≥α(2-α)\|P_Lv_L\|^2-α^2 C_L\|M_L v_L\|^2. (63) Since α>0α>0, to guarantee Δℰ2(α)>0 ^2(α)>0, the inequality requires: (2−α)‖LL‖2−αC~L‖L†L‖2>0.(2-α)\|P_Lv_L\|^2-α C_L\|M_L v_L\|^2>0. (64) Because we have ‖LL‖2>0\|P_Lv_L\|^2>0, ‖LL‖2+C~L‖L†L‖2\|P_Lv_L\|^2+ C_L\|M_L v_L\|^2 is guaranteed to be positive. Thus, we can obtain the upper bound for α, denoted as α¯ α: α¯=min1,2‖LL‖2‖LL‖2+C~L‖L†L‖2,0<α<α¯. α= \1, 2\|P_Lv_L\|^2\|P_Lv_L\|^2+ C_L\|M_L v_L\|^2 \, 0<α< α. (65) This construction inherently ensures α¯>0 α>0. Consequently, for any assigned α within the constructed interval α∈(0,α¯)α∈(0, α), we have Δℰ2(α)=c>0 ^2(α)=c>0. As the relaxation parameter r approaches 11, the threshold δL=(1−r−μL)ℰL−12 _L=(1-r- _L)E_L-1^2 approaches 00. Thus, since c>0c>0, there always exists an r∈(0,1)r∈(0,1) such that δL≤c _L≤ c. This ensures that the descent inequality Δℰ2(α)≥δL ^2(α)≥ _L can always be rigorously satisfied. Consequently, we obtain: ℰL2≤ℰL−12−δL=(r+μL)ℰL−12.E_L^2 _L-1^2- _L=(r+ _L)E_L-1^2. (66) Note that limL→+∞μL=0 _L→+∞ _L=0. By utilizing (66), it strictly follows that limL→+∞ℰL2=0 _L→+∞E_L^2=0. ∎ Remark 2 In PI-SC-I, explicitly computing α¯ α is unnecessary. Instead, α can be dynamically adapted utilizing a decay strategy analogous to the relaxation parameter r. Notably, in most empirical scenarios—particularly during the initial phases of training—the first-order linear descent heavily dominates the second-order remainder bound, yielding ‖LL‖2≫C~L‖L†L‖2\|P_Lv_L\|^2 C_L\|M_L v_L\|^2. Under this condition, the optimal α converges to 11. Consequently, the framework simply initializes α=1α=1 and decays it exclusively when the candidate nodes fail to achieve residual descent, thereby guaranteeing feasible node configuration with minimal computational overhead. Building upon the universal approximation property of the PI-SC-I algorithm established in Theorem 1, we now extend our theoretical analysis to the globally updated PI-SC-I algorithm. Recalling the optimal localized configuration loc∗ β_loc^* and the global configuration ~glob∗ β_glob^* established in Section I-A, let ℰL2()=‖~L()‖2E_L^2( β)=\| e_L( β)\|^2 denote the comprehensive true residual norm evaluated at any candidate configuration β. Consequently, the local residual norm is explicitly defined as ℰL,loc2=ℰL2(loc∗)E_L,loc^2=E_L^2( β_loc^*). Theorem 2 Suppose the conditions of Theorem 1 hold. We define the global search direction as Δ=~glob∗−loc∗ β= β_glob^*- β_loc^*. Then, there exists a scaling factor α¯∈[0,1] α∈[0,1] such that for any α∈[0,α¯]α∈[0, α], the adopted configuration, defined as =loc∗+αΔ, β= β_loc^*+α β, (67) guarantees that the following inequality holds: ℰL2()≤ℰL,loc2.E_L^2( β) _L,loc^2. (68) By accepting the configuration that satisfies this inequality, the sequence of ℰL2\E_L^2\ monotonically decreases, satisfying limL→+∞ℰL2=0 _L→+∞E_L^2=0. Proof: If ~glob∗=loc∗ β_glob^*= β_loc^*, we have Δ= β=0. Consequently, the step size inherently yields ℰL2()=ℰL,loc2E_L^2( β)=E_L,loc^2 for any assigned α. In this scenario, taking α∗=0α^*=0 constitutes a valid trivial solution that satisfies the acceptance criterion. If ~glob∗≠loc∗ β_glob^*≠ β_loc^*, Δ β constitutes a strict descent direction for the linearized least-squares objective. Analogous to the proof of Theorem 1, the first-order linear descent dominates the bounded second-order remainder for sufficiently small α>0α>0. This guarantees the existence of a positive upper bound α¯∈(0,1] α∈(0,1] such that for all α∈(0,α¯]α∈(0, α], we have ℰL2()<ℰL,loc2E_L^2( β)<E_L,loc^2. If the two configurations attain the same minimum value of the linearized objective, taking α¯=0 α=0 satisfies the acceptance criterion. Since the boundary case α=0α=0 inherently satisfies the equality ℰL2()=ℰL,loc2E_L^2( β)=E_L,loc^2, the inequality described in (68) holds for all α∈[0,α¯]α∈[0, α]. Recalling from (66) that ℰL,loc2≤(r+μL)ℰL−12E_L,loc^2≤(r+ _L)E_L-1^2. Thus, we have: ℰL2≤ℰL,loc2≤(r+μL)ℰL−12.E_L^2 _L,loc^2≤(r+ _L)E_L-1^2. (69) Using the same arguments established in the proof of Theorem 1, we conclude that limL→+∞ℰL2=0 _L→+∞E_L^2=0. ∎ Following the theoretical guarantees of PI-SC-I, we now extend our theoretical analysis to the PI-SC-I algorithm. Since this algorithm employs an identical discrete line search mechanism, for brevity, we omit the proof and directly present its universal approximation theorem. We recall the windowed configuration ~win∗ β_win^* defined in Section I-A. Theorem 3 Suppose the conditions of Theorem 1 hold. We define the windowed search direction as Δ=~win∗−loc∗ β= β_win^*- β_loc^*. Then, there exists a scaling factor α¯∈[0,1] α∈[0,1] such that for any α∈[0,α¯]α∈[0, α], the adopted configuration, defined as =loc∗+αΔ, β= β_loc^*+α β, (70) guarantees that the following inequality holds: ℰL2()≤ℰL,loc2.E_L^2( β) _L,loc^2. (71) By accepting the configuration that satisfies this inequality, the sequence of ℰL2\E_L^2\ monotonically decreases, satisfying limL→+∞ℰL2=0 _L→+∞E_L^2=0. I-D Hyperparameter Configuration Strategy Implementing the PI-SCM framework relies on the appropriate configuration of the hyperparameter sets, particularly the candidate pool for the stochastic parameter sampling range Υ=τ1,τ2,…,τend =\ _1, _2,…, _end\. Unlike conventional backpropagation-based deep learning where weight initialization merely dictates the starting point of optimization, the sampling bounds L∼[−τ,τ]nw_L [-τ,τ]^n and bL∼[−τ,τ]b_L [-τ,τ] govern the basis function space generated at each incremental step. Consistent with the original SCN theory [13], adopting a progressively larger τ empowers the randomized parameters Lw_L and bLb_L to explore a broader configuration space, thereby generating basis functions with stronger nonlinear expressivity. However, extending this mechanism to PI-SCM introduces a structural bottleneck absent in purely data-driven configurations. Taking the localized PI-SC-I algorithm as an illustrative example, the construction of L()A_L(x) in (14) necessitates the analytical evaluation of the hidden-node derivatives up to order q. Given the node output hL()=g(L⊤+bL)h_L(x)=g(w_L x+b_L) and an arbitrary multi-index ∈ℐq γ _q, repeated application of the chain rule yields: DhL()=g(||)(L⊤+bL)∏j=1nwL,jγj,D γh_L(x)=g^(| γ|)(w_L x+b_L) _j=1^nw_L,j _j, (72) where wL,jw_L,j is the j-th component of the hidden weights Lw_L, and g(||)g^(| γ|) denotes the corresponding activation derivative. Recalling (15), we have: L()=∑∈ℐqD(L−1)g(||)(L⊤+bL)∏j=1nwL,jγj.A_L(x)= _ γ _qJ_D γu(u_L-1)g^(| γ|)(w_L x+b_L) _j=1^nw_L,j _j. (73) For the smooth bounded activation functions adopted in this work, the magnitude of each derivative contribution within L()A_L(x) is governed primarily by the product ∏j=1nwL,jγj _j=1^nw_L,j _j. Consequently, a derivative block of order ||| γ| scales as (τ||)O(τ^| γ|). When a disproportionately large τ is assigned, the highest-order derivative terms can dominate the lower-order state-dependent terms. This column-wise magnitude disparity inflates the condition number of LM_L defined in (24), and the effect becomes progressively stronger as q increases. Consequently, a large τ exacerbates matrix ill-conditioning. While the TSVD approach introduced in Section I-A mitigates the adverse effects of matrix ill-conditioning, feeding it highly ill-conditioned matrices forces the truncation of excessive singular components, destroying the physical information embedded within the Jacobian operators. Conversely, restricting τ to small values guarantees well-conditioned matrices but restricts the expressive capacity of the network. For complex problems, relying exclusively on such constrained spaces requires a large number of hidden nodes to achieve acceptable accuracy, leading to slow convergence and structural redundancy. To resolve the tension among convergence rate, structural exploratory capacity, and numerical matrix stability, the candidate pool is formulated as a parametrically scaled monotonically increasing sequence, defined as Υ=a×1,2,4,8,16,32,64 =a×\1,2,4,8,16,32,64\. Here, a>0a>0 serves as a problem-dependent characteristic scaling factor. The magnitude of a is specifically tailored to the governing differential equations to equilibrate the numerical magnitudes among the state-dependent Jacobian J_u and the derivative-dependent Jacobians DJ_D γu for 1≤||≤q1≤| γ|≤ q. PI-SCM processes this scaled sequence in ascending order. The framework samples from the narrowest bounds first, escalating to the subsequent broader candidate range only when the current configuration space fails to secure a valid residual descent. This strategy effectively balances numerical stability with the structural expressivity required for complex physical dynamics. IV Experimental Results This section evaluates PI-SCM along three complementary axes. First, forward problems assess state approximation for the Van der Pol oscillator, the two-dimensional Helmholtz equation, and the Allen–Cahn equation. Second, inverse problems assess joint state reconstruction and parameter identification on the Van der Pol and Helmholtz systems. Third, an activation-function study isolates the effect of the hidden-node basis. IV-A Experimental Systems IV-A1 Van der Pol oscillator The Van der Pol oscillator [25] is a nonlinear, non-conservative dynamical system used in circuit modeling, seismology, biological neuron modeling, and optimal control [26]. We consider the controlled second-order equation utt−μ(1−u2)ut+u−U=0,(t,U)∈[0,5]×[0,1],u_t-μ(1-u^2)u_t+u-U=0, (t,U)∈[0,5]×[0,1], (74) where μ is the nonlinear damping coefficient and U is the exogenous input. The initial conditions are u(0,U)=−0.25u(0,U)=-0.25 and ut(0,U)=−2.0u_t(0,U)=-2.0. The value μ=1μ=1 is prescribed in the forward task, whereas μ is identified in the inverse task. IV-A2 Two-dimensional Helmholtz equation The Helmholtz equation is a canonical elliptic model of time-harmonic wave propagation. The benchmark is uxx+uyy+k2u−q(x,y)=0,(x,y)∈[−1,1]×[−1,1],u_x+u_y+k^2u-q(x,y)=0, (x,y)∈[-1,1]×[-1,1], (75) with homogeneous Dirichlet conditions u(−1,y)=u(1,y)=u(x,−1)=u(x,1)=0.u(-1,y)=u(1,y)=u(x,-1)=u(x,1)=0. (76) A manufactured source is used: q(x,y)=−(a1π)2sin(a1πx)sin(a2πy)−(a2π)2sin(a1πx)sin(a2πy)+sin(a1πx)sin(a2πy), splitq(x,y)=&-(a_1π)^2 (a_1π x) (a_2π y)\\ &-(a_2π)^2 (a_1π x) (a_2π y)\\ &+ (a_1π x) (a_2π y), split (77) where a1=1a_1=1 and a2=4a_2=4. The value k2=1k^2=1 is prescribed in the forward task, whereas k2k^2 is identified in the inverse task. The corresponding analytical solution is u(x,y)=sin(a1πx)sin(a2πy).u(x,y)= (a_1π x) (a_2π y). (78) IV-A3 Allen–Cahn equation The Allen–Cahn equation models phase separation and presents a challenging combination of a stiff reaction term and sharp moving interfaces. We consider ut−ϵuxx+5u3−5u=0,(x,t)∈[−1,1]×[0,0.5],u_t-ε u_x+5u^3-5u=0, (x,t)∈[-1,1]×[0,0.5], (79) where ϵ=10−4ε=10^-4, u(x,0)=x2cos(πx)u(x,0)=x^2 (π x), and the periodic boundary constraints are u(−1,t)=u(1,t)u(-1,t)=u(1,t) and ux(−1,t)=ux(1,t)u_x(-1,t)=u_x(1,t). Because the small diffusion coefficient produces a large scale disparity between the derivative and reaction terms, we define X=ωnxX= _nx with ωn=1/ϵ=100 _n=1/ ε=100. The transformed equation is ∂u∂t−ϵωn2∂2u∂X2+5u3−5u=0. ∂ u∂ t-ε _n^2 ∂^2u∂ X^2+5u^3-5u=0. (80) Since ϵωn2=1ε _n^2=1, the coordinate stretching balances the state- and derivative-dependent Jacobian contributions in (73) and improves the conditioning of LM_L without changing the underlying physical solution. For each system, state accuracy is measured by the root mean square error (RMSE) over NtestN_test test inputs: RMSE=1Ntest∑i=1Ntest‖pred(i)−true(i)‖22,RMSE= 1N_test _i=1^N_test \|u_pred( ζ_i)-u_true( ζ_i) \|_2^2, (81) where i ζ_i denotes the iith test input, and predu_pred and trueu_true denote the predicted and reference states, respectively. All methods are implemented in Python and evaluated on a workstation equipped with an NVIDIA RTX 4060 Ti GPU, an Intel Core i5-14600KF CPU, and 16 GB of RAM. IV-B Forward Problems: Experimental Setup We consider both unsupervised and semi-supervised settings. For the Van der Pol system, we sample Np=1000N_p=1000 interior points and Ni=100N_i=100 initial points, with an additional Nd=100N_d=100 Runge–Kutta observations in the semi-supervised setting. The penalty weights are set to ωbc=1 _bc=1 and ωdata=1 _data=1. The constructive parameters are Υ=15×1,2,4,8,16,32,64 =15×\1,2,4,8,16,32,64\, Tmax=50T_max=50, Lmax=300L_max=300, and W=200W=200, with 1=1.0,0.9,…,0.1A_1=\1.0,0.9,…,0.1\, 2=1.0,0.9,…,0A_2=\1.0,0.9,…,0\, and ℛ=1−10−k∣k=1,…,7R=\1-10^-k k=1,…,7\. For the two-dimensional Helmholtz equation, we use Np=2000N_p=2000 interior collocation points, Nbc=1000N_bc=1000 boundary points, and Lmax=400L_max=400. The remaining PI-SCM settings follow those of the Van der Pol experiment. The baselines include shallow and deep PINNs with architectures 1×4001× 400 and 4×644× 64, respectively, together with the physics-informed extreme learning machine (PIELM) [27], which uses the same 1×4001× 400 architecture as PI-SCM. For the Allen–Cahn equation, we set Lmax=500L_max=500, while the other PI-SCM settings follow those of the Van der Pol experiment. The corresponding variable-scaling physics-informed neural network (VS-PINN) baselines [28] use shallow and deep architectures of 1×5001× 500 and 4×644× 64, respectively. Across all forward experiments, each PINN-based baseline uses the same physics collocation, initial/boundary, and observational points as the corresponding PI-SCM experiment. All PINN-based baselines employ the hyperbolic tangent activation and are trained with Adam for 5,000 epochs, with an initial learning rate of 10−210^-2 for the Van der Pol and Helmholtz equations and 10−310^-3 for the Allen–Cahn equation. For the Van der Pol equation, we additionally compare with PISCN [23], which adopts the same constructive settings as PI-SCM and is subsequently optimized with Adam for 100 epochs using a learning rate of 5×10−35× 10^-3. Since PISCN requires labeled data during construction, it is evaluated only in the semi-supervised setting. PI-SCM and PISCN use the sine activation. Each configuration is repeated 100 times. IV-C Forward Problems: Results All tabulated values are reported as mean ± standard deviation, and predictive accuracy is evaluated on 10,000 randomly sampled test coordinates. For each system, the table reports training time and RMSE under the unsupervised and semi-supervised paradigms, while the accompanying figure presents the semi-supervised PI-SC-I prediction, absolute error, and representative cross-sectional profiles. IV-C1 Van der Pol oscillator The corresponding results are presented in Table I and Fig. 1, respectively. TABLE I: Forward results for the Van der Pol oscillator. Algorithm Paradigm Time (s) RMSE PINN (1×3001× 300) Unsupervised 21.3±0.5721.3± 0.57 3.69e−1±9.18e−23.69e-1± 9.18e-2 Semi-supervised 27.5±0.8927.5± 0.89 8.05e−2±4.02e−38.05e-2± 4.02e-3 PINN (4×644× 64) Unsupervised 40.0±1.0240.0± 1.02 5.49e−3±3.03e−35.49e-3± 3.03e-3 Semi-supervised 50.8±1.0850.8± 1.08 1.34e−3±5.22e−41.34e-3± 5.22e-4 PISCN Unsupervised – – Semi-supervised 13.8±0.3713.8± 0.37 7.58e−2±2.59e−27.58e-2± 2.59e-2 PI-SC-I Unsupervised 0.53±0.020.53± 0.02 3.88e−1±1.41e−13.88e-1± 1.41e-1 Semi-supervised 0.60±0.010.60± 0.01 1.23e−1±3.82e−21.23e-1± 3.82e-2 PI-SC-I Unsupervised 1.37±0.041.37± 0.04 1.30e−2±5.73e−31.30e-2± 5.73e-3 Semi-supervised 1.57±0.051.57± 0.05 5.27e−3±2.63e−35.27e-3± 2.63e-3 PI-SC-I Unsupervised 1.44±0.041.44± 0.04 1.93e−3±9.86e−41.93e-3± 9.86e-4 Semi-supervised 1.64±0.041.64± 0.04 5.92e−4±1.76e−45.92e-4± 1.76e-4 Fig. 1: Forward solution of the Van der Pol oscillator obtained by semi-supervised PI-SC-I. Top left: predicted state over (t,U)(t,U). Top right: absolute pointwise error. Bottom: cross-sectional trajectories at U=0.3U=0.3 and U=0.6U=0.6. IV-C2 Two-dimensional Helmholtz equation The corresponding results are presented in Table I and Fig. 2, respectively. TABLE I: Forward results for the two-dimensional Helmholtz equation. Algorithm Paradigm Time (s) RMSE PINN (1×4001× 400) Unsupervised 40.5±1.2340.5± 1.23 2.62e−2±1.97e−22.62e-2± 1.97e-2 Semi-supervised 53.2±1.5153.2± 1.51 1.13e−2±4.50e−31.13e-2± 4.50e-3 PINN (4×644× 64) Unsupervised 75.6±1.1175.6± 1.11 3.28e−2±1.43e−23.28e-2± 1.43e-2 Semi-supervised 105±1.38105± 1.38 1.44e−2±5.06e−31.44e-2± 5.06e-3 PIELM Unsupervised 0.02±0.000.02± 0.00 9.09e−3±3.83e−39.09e-3± 3.83e-3 Semi-supervised 0.02±0.000.02± 0.00 6.61e−3±3.84e−36.61e-3± 3.84e-3 PI-SC-I Unsupervised 0.57±0.010.57± 0.01 1.91e−1±9.82e−21.91e-1± 9.82e-2 Semi-supervised 0.66±0.010.66± 0.01 1.57e−1±7.25e−21.57e-1± 7.25e-2 PI-SC-I Unsupervised 1.99±0.051.99± 0.05 5.19e−3±2.05e−35.19e-3± 2.05e-3 Semi-supervised 2.20±0.062.20± 0.06 4.71e−3±1.78e−34.71e-3± 1.78e-3 PI-SC-I Unsupervised 2.31±0.042.31± 0.04 1.02e−6±4.79e−71.02e-6± 4.79e-7 Semi-supervised 2.54±0.052.54± 0.05 1.04e−6±8.75e−71.04e-6± 8.75e-7 Fig. 2: Forward solution of the two-dimensional Helmholtz equation obtained by PI-SC-I. Top left: predicted state over (x,y)(x,y). Top right: absolute pointwise error. Bottom: cross-sectional profiles at x=−0.5x=-0.5 and x=0.5x=0.5. IV-C3 Allen–Cahn equation The corresponding results are presented in Table I and Fig. 3, respectively. TABLE I: Forward results for the Allen–Cahn equation. Algorithm Paradigm Time (s) RMSE VS-PINN (1×5001× 500) Unsupervised 31.8±0.7931.8± 0.79 5.40e−1±2.63e−25.40e-1± 2.63e-2 Semi-supervised 36.2±0.9436.2± 0.94 4.65e−1±4.67e−24.65e-1± 4.67e-2 VS-PINN (4×644× 64) Unsupervised 65.0±1.0665.0± 1.06 5.23e−1±4.14e−25.23e-1± 4.14e-2 Semi-supervised 71.6±1.3571.6± 1.35 7.75e−2±1.06e−17.75e-2± 1.06e-1 PI-SC-I Unsupervised 1.59±0.041.59± 0.04 5.28e−1±3.50e−35.28e-1± 3.50e-3 Semi-supervised 1.73±0.031.73± 0.03 3.16e−1±2.47e−23.16e-1± 2.47e-2 PI-SC-I Unsupervised 3.13±0.073.13± 0.07 1.71e−1±4.44e−21.71e-1± 4.44e-2 Semi-supervised 3.41±0.063.41± 0.06 2.25e−2±1.04e−22.25e-2± 1.04e-2 PI-SC-I Unsupervised 3.49±0.143.49± 0.14 6.16e−3±2.07e−36.16e-3± 2.07e-3 Semi-supervised 3.70±0.073.70± 0.07 1.57e−3±2.88e−41.57e-3± 2.88e-4 Fig. 3: Forward solution of the Allen–Cahn equation obtained by semi-supervised PI-SC-I. Top left: predicted state over (t,x)(t,x). Top right: absolute pointwise error. Bottom: spatial profiles at t=0.1t=0.1 and t=0.4t=0.4. IV-D Forward Problems: Result Analysis The forward results in Tables I–I reveal a consistent accuracy–efficiency hierarchy across the three forward benchmarks. PI-SC-I provides the lowest construction cost but may underfit problems involving strong oscillations, high spatial frequencies, or sharp transitions. Extending the correction scope in PI-SC-I substantially improves predictive accuracy, while the global correction in PI-SC-I consistently achieves the lowest errors among the three variants. Importantly, this gain preserves the main computational advantage of the framework, with PI-SC-I remaining substantially faster than the deep gradient-based baselines while achieving comparable or better accuracy. This trend is consistent across both ODE and PDE problems. For the unsupervised Van der Pol problem, PI-SC-I reduces the RMSE from 5.49×10−35.49× 10^-3 for the deep PINN to 1.93×10−31.93× 10^-3, while reducing the training time from 40.040.0 s to 1.441.44 s. For the Helmholtz equation, it achieves an RMSE on the order of 10−610^-6 under both training paradigms. On the semi-supervised Allen–Cahn problem, PI-SC-I reaches an RMSE of 1.57×10−31.57× 10^-3, reducing the error by approximately a factor of 4949 while being 19×19× faster than the corresponding deep VS-PINN. The particularly high accuracy observed for the Helmholtz problem can be attributed to its linear operator structure. Because the governing equation is linear in both u and its derivatives uxxu_x and uyyu_y, the state-space linearization used in PI-SCM is exact for this operator and introduces no nonlinear truncation error. This enables the global least-squares correction in PI-SC-I to recover the solution with very high precision. For the same reason, PIELM is considered only for the Helmholtz benchmark, since it is restricted to differential equations that are linear in the state and its derivatives. Although PIELM is faster for this problem, its RMSE remains more than three orders of magnitude higher than that of PI-SC-I. Sparse observations provide additional benefits when the physics constraints alone pose a challenging approximation problem, as demonstrated by the Van der Pol and Allen–Cahn results. Their contribution is limited for the Helmholtz problem because the unsupervised model already achieves a very low error level. Overall, the forward experiments demonstrate that PI-SCM provides the most favorable balance among predictive accuracy, computational efficiency, and the ability to incorporate sparse observations. IV-E Inverse Problems: Experimental Setup The inverse experiments jointly recover the state and one scalar physical coefficient. Each system compares the two PINN architectures with PI-SC-I, PI-SC-I, and PI-SC-I. For both inverse benchmarks, the unknown physical parameter has a ground-truth value of 1.01.0 and is initialized at 0.10.1: μtrue=1.0 _true=1.0 and μ(0)=0.1μ^(0)=0.1 for Van der Pol, and ktrue2=1.0k^2_true=1.0 and (k2)(0)=0.1 (k^2 )^(0)=0.1 for Helmholtz. Each experiment uses Nd=100N_d=100 observations, while all other sampling, construction, and baseline settings match the corresponding forward experiment. Each method is evaluated over 100 independent trials IV-F Inverse Problems: Results All tabulated values are reported as mean ± standard deviation, and state predictive accuracy is evaluated on 10,000 randomly sampled test coordinates. For each system, the table reports parameter identification error, state RMSE, and training time. IV-F1 Van der Pol oscillator The corresponding results are presented in Table IV. TABLE IV: Inverse results for the Van der Pol oscillator. Algorithm |μpred−μtrue|| _pred- _true| Time (s) State RMSE PINN (1×3001× 300) 2.88e−1±1.09e−12.88e-1± 1.09e-1 22.0±0.5122.0± 0.51 5.85e−2±1.35e−25.85e-2± 1.35e-2 PINN (4×644× 64) 2.35e−4±2.01e−42.35e-4± 2.01e-4 51.8±0.9451.8± 0.94 9.51e−4±3.42e−49.51e-4± 3.42e-4 PI-SC-I 4.01e−1±4.69e−24.01e-1± 4.69e-2 0.71±0.020.71± 0.02 1.08e−1±2.05e−21.08e-1± 2.05e-2 PI-SC-I 1.91e−3±1.59e−31.91e-3± 1.59e-3 1.64±0.031.64± 0.03 4.80e−3±2.10e−34.80e-3± 2.10e-3 PI-SC-I 5.08e−4±2.75e−45.08e-4± 2.75e-4 1.95±0.041.95± 0.04 5.96e−4±2.00e−45.96e-4± 2.00e-4 IV-F2 Two-dimensional Helmholtz equation The corresponding results are presented in Table V. TABLE V: Inverse results for the two-dimensional Helmholtz equation. Algorithm |kpred2−ktrue2||k^2_pred-k^2_true| Time (s) State RMSE PINN (1×4001× 400) 1.09e+1±2.87e+01.09e+1± 2.87e+0 50.6±1.1950.6± 1.19 2.82e−1±1.27e−12.82e-1± 1.27e-1 PINN (4×644× 64) 6.72e−1±6.27e−16.72e-1± 6.27e-1 101±2.13101± 2.13 2.36e−2±1.22e−22.36e-2± 1.22e-2 PI-SC-I 3.72e+1±1.93e+13.72e+1± 1.93e+1 1.49±0.031.49± 0.03 2.17e−1±9.90e−22.17e-1± 9.90e-2 PI-SC-I 7.14e−1±7.55e−17.14e-1± 7.55e-1 3.14±0.073.14± 0.07 5.04e−3±2.21e−35.04e-3± 2.21e-3 PI-SC-I 8.42e−4±1.60e−38.42e-4± 1.60e-3 3.51±0.103.51± 0.10 6.70e−5±1.15e−46.70e-5± 1.15e-4 IV-G Inverse Problems: Result Analysis The inverse results in Tables IV and V show the same update-scope hierarchy as in the forward experiments. PI-SC-I is less effective for joint state and parameter recovery, while PI-SC-I improves the reconstruction accuracy through windowed correction. PI-SC-I provides the most accurate overall recovery. Across both inverse benchmarks, PI-SC-I provides parameter-identification accuracy comparable to or better than that of the deep PINN. It also improves state reconstruction and requires substantially less training time. Overall, these results demonstrate that PI-SCM retains high predictive accuracy and computational efficiency in inverse problems. IV-H Activation-Function Study: Experimental Setup We isolate the influence of the hidden-node basis on the semi-supervised Van der Pol forward problem. PI-SC-I is evaluated with sine, hyperbolic tangent, and sigmoid activations. LmaxL_max is set to 1,000, while the sampling points, penalty weights, candidate ranges, and remaining constructive settings are kept identical to the semi-supervised Van der Pol configuration described above. Each activation is evaluated over 100 independent trials. IV-I Activation-Function Study: Results Table VI summarizes the state RMSE , while Fig. 4 shows the evolution of the average true residual ℰL2E_L^2 with the number of accepted hidden nodes. TABLE VI: Semi-supervised PI-SC-I results with different activation functions. Activation RMSE Sine 5.66e−4±3.59e−45.66e-4± 3.59e-4 Hyperbolic tangent 2.71e−3±8.92e−42.71e-3± 8.92e-4 Sigmoid 1.68e−2±7.96e−31.68e-2± 7.96e-3 Fig. 4: Evolution of the average true residual ℰL2E_L^2 with the number of accepted hidden nodes for three activation functions. IV-J Activation-Function Study: Result Analysis Sine activation achieves the highest predictive accuracy and the fastest residual convergence among the tested activation functions. As shown in Fig. 4, it reduces the nonlinear residual more rapidly and converges to the lowest residual level, supporting its use as the default activation in the reported PI-SCM experiments. V Conclusion This paper introduced PI-SCM, a backpropagation-free framework for nonlinear differential equations. Its progressive algorithmic suite, ranging from localized construction (PI-SC-I) through windowed correction (PI-SC-I) to global correction (PI-SC-I), converts nonlinear physics-informed learning into a sequence of analytically solvable least-squares problems. The universal approximation analysis establishes the theoretical basis for residual convergence, while the inverse formulation enables state reconstruction and physical-parameter identification within the same framework. The empirical study considered forward approximation on the Van der Pol, two-dimensional Helmholtz, and stiff Allen–Cahn systems; inverse identification on the Van der Pol and Helmholtz systems; and an activation-function ablation. Across these tasks, PI-SC-I consistently provided the strongest accuracy–efficiency trade-off among the proposed variants. It achieved competitive or lower errors than deep gradient-based PINNs while reducing training time by one to two orders of magnitude. The semi-supervised results also showed that PI-SCM can assimilate sparse observations when the physical constraints alone are insufficient, and the activation study identified sine as the most effective basis for the oscillatory benchmark. Future work should examine scalability to higher-dimensional geometries, more severe multiscale stiffness, noisy parameter-identification settings, and adaptive selection of the candidate range and update scope. Extending the current shallow constructive architecture while retaining analytical training and numerical stability is another important direction toward real-time scientific machine-learning applications. References [1] Baker N, Alexander F, Bremer T, et al. Workshop report on basic research needs for scientific machine learning: Core technologies for artificial intelligence[R]. USDOE Office of Science (SC), Washington, DC (United States), 2019. [2] Raissi M, Perdikaris P, Karniadakis G E. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[J]. Journal of Computational Physics, 2019, 378: 686-707. [3] Schiassi E , D’Ambrosio A , Drozd K ,et al. Physics-Informed Neural Networks for Optimal Planar Orbit Transfers[J]. Journal of Spacecraft and Rockets, 2022, 59(3):16.DOI:10.2514/1.A35138. [4] Shukla K, Di Leoni P C, Blackshire J, et al. Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks[J]. Journal of Nondestructive Evaluation, 2020, 39(3): 61. [5] Berkhahn S, Ehrhardt M. A physics-informed neural network to model COVID-19 infection and hospitalization scenarios[J]. Advances in Continuous and Discrete Models, 2022, 2022(1): 61. [6] Diao Y, Yang J, Zhang Y, et al. Solving multi-material problems in solid mechanics using physics-informed neural networks based on domain decomposition technology[J]. Computer Methods in Applied Mechanics and Engineering, 2023, 413: 116120. [7] Wang S, Sankaran S, Wang H, et al. An expert’s guide to training physics-informed neural networks. arXiv preprint arXiv:2308.08468, 2023. [8] Cuomo S, Di Cola V S, Giampaolo F, et al. Scientific machine learning through physics–informed neural networks: Where we are and what’s next[J]. Journal of Scientific Computing, 2022, 92(3): 88. [9] Pao Y H, Takefuji Y. Functional-link net computing: theory, system architecture, and functionalities[J]. Computer, 1992, 25(5): 76-79. [10] Schmidt W F, Kraaijveld M A, Duin R P W. Feed forward neural networks with random weights[C]//International Conference on Pattern Recognition. IEEE Computer Society Press, 1992: 1-4. [11] Huang G B, Zhu Q Y, Siew C K. Extreme learning machine: theory and applications[J]. Neurocomputing, 2006, 70(1-3): 489-501. [12] Liu D, Chang T S, Zhang Y. A constructive algorithm for feedforward neural networks with incremental training[J]. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 2002, 49(12): 1876-1879. [13] Wang D, Li M. Stochastic configuration networks: Fundamentals and algorithms[J]. IEEE Transactions on Cybernetics, 2017, 47(10): 3466-3479. [14] Wang D, Li M. Deep stochastic configuration networks with universal approximation property[C]//2018 International Joint Conference on Neural Networks (IJCNN). IEEE, 2018: 1-8. [15] Li K, Wang W, Lin S. Soft measurement of ammonia nitrogen concentration based on GA-SCN[C]. IEEE Symposium on Product Compliance Engineering-Asia (ISPCE-CN), IEEE, 2018: 1-4. [16] Felicetti M J, Wang D. Deep stochastic configuration networks with optimised model and hyper-parameters[J]. Information Sciences, 2022, 600: 431-441. [17] Felicetti M J, Wang D. Deep stochastic configuration networks with different random sampling strategies[J]. Information Sciences, 2022, 607: 819-830. [18] Pan J, Luan F, Gao Y, et al. FPGA-based implementation of stochastic configuration network for robotic grasping recognition[J]. IEEE Access, 2020, 8: 139966-139973. [19] Li W, Tao H, Li H, Chen K, et al. Greengage grading using stochastic configuration networks and a semi-supervised feedback mechanism[J]. Information Sciences, 2019, 488: 1-12. [20] Li W, Zhang Q, Wang D, et al. Stochastic configuration networks for self-blast state recognition of glass insulators with adaptive depth and multi-scale representation[J]. Information Sciences, 2022, 604: 61-79. [21] Chen K, An J, Fang Y, et al. Research on solid waste plastic bottle cognitive based on YOLOv5s and deep stochastic configuration network[C]. International Conference on Automation, Control and Robotics Engineering, 2022: 275-280. [22] Liu W, Ren C, Xu Y. PV generation forecasting with missing input data: a super-resolution perception approach[J]. IEEE Transactions on Sustainable Energy, 2021, 12(2): 1493-1496. [23] Xu L, Yang C, Xu X, et al. Physics-informed stochastic configuration network promoted model predictive control with multi-objective optimization[J]. Artificial Intelligence Review, 2025, 58(9): 281. [24] Wu L, Tan W G Y, Zhou L, Braatz R D, Drgoňa J. Least-squares multi-step Koopman operator learning for model predictive control. arXiv preprint arXiv:2601.11901, 2026. [25] Hafeez H Y, Ndikilar C E, Isyaku S. Analytical study of the van der pol equation in the autonomous regime[J]. Progress in Physics, 2015, 11: 252. [26] Andersson J, Åkesson J, Diehl M. Dynamic optimization with CasADi[C]//2012 IEEE 51st IEEE Conference on Decision and Control (CDC). IEEE, 2012: 681-686. [27] Dwivedi V, Srinivasan B. Physics-informed extreme learning machine (PIELM): A rapid method for the numerical solution of partial differential equations[J]. Neurocomputing, 2020, 391: 96-118. [28] Ko S, Park S. VS-PINN: A fast and efficient training of physics-informed neural networks using variable-scaling methods for solving PDEs with stiff behavior[J]. Journal of Computational Physics, 2025, 529: 113860.