Application of a Neural-Network Algorithm Based on Orthogonal Laguerre Polynomials to the Solution of the Cauchy Problem for Differential Equations with a Variable-Order Fractional Derivative on the Example of the Fractional Lane–Emden Equation

1. Introduction

Modeling physical processes using differential equations often encounters limitations of the classical mathematical apparatus. In recent decades, there has been increased interest in using fractional calculus to describe systems with memory and self-similarity effects [1,2]. Unlike classical integer-order derivatives, fractional operators allow one to take into account hereditary properties and nonlocal interactions, which makes them a powerful tool for modeling a wide range of phenomena—from viscoelasticity to anomalous diffusion to spread of viruses [3,4,5,6,7]. A special place among them is occupied by variable-order operators, where the exponent of the fractional derivative is itself a function of the spatial or temporal coordinate. Such operators provide a more flexible description of systems whose rheological or diffusive properties change depending on the state or position in space [8,9,10]. One of the actively developing directions is the application of fractional calculus to problems of mathematical physics and astrophysics.

One of the important models of mathematical physics is the Lane–Emden equation [11,12,13,14], which describes the equilibrium of spherical gas clouds in stellar astrophysics. The Lane–Emden equation arises in the study of polytropic stars, the density distribution in isothermal gas spheres, and in the theory of gravitational equilibrium [15]. The introduction of a fractional derivative into this equation allows one to take into account nonlocal properties of gravitational interaction and complex thermodynamic processes, which makes the model more adequate for describing real astrophysical objects. In recent years, a number of numerical and approximate methods have been proposed for solving fractional Lane–Emden equations, including variational-iteration methods [16], spectral methods based on Legendre polynomials [17], collocation methods with wavelets [18], and other approaches [19,20]. However, most of these works consider either an integer order of the derivative or a constant fractional order, which limits their applicability to problems with variable memory.

The analytical solution of the Cauchy problem for the fractional singular nonlinear Lane–Emden differential equation exists only in particular cases (for example, for an integer exponent), which necessitates the development of efficient numerical methods [21]. Traditional approaches, such as finite-difference schemes, encounter difficulties due to the singularity of the solution near zero and the nonlocal nature of the fractional derivative. This requires the construction of special grids and leads to an increase in computational costs. In the case of a variable order of the fractional derivative, the situation becomes even more complicated, since the kernel of the fractional operator depends on the coordinate, and classical difference schemes require adaptation to each value of the exponent [22,23,24,25,26]. The problem of increasing computational costs in the numerical solution of fractional differential equations of both constant and variable order can be partially solved by parallelization, implementing both CPU versions of such algorithms and hybrid CPU-GPU analogues [27]. Accelerating computations is important when solving coefficient inverse problems for equations with fractional derivatives [28]. However, if such a numerical method has an unavoidable sequential part, the parallelization efficiency may be low.

In parallel with the development of fractional calculus, there is a rapid development of machine-learning methods. Neural-network approaches to solving differential equations (Physics-Informed Neural Networks, PINNs) have shown their effectiveness for a wide class of problems [29,30,31]. They allow one to obtain an approximate solution in analytical form without requiring grid construction and easily adapt to various types of equations and boundary conditions [32]. PINNs are successfully applied in deformable solid mechanics, heat-transfer problems, aerodynamics, and other fields [33,34]. Among the most recent developments in this direction, we note the VooPINN architecture [35]—a physics-informed neural network with a variable-order operator designed for forward and inverse problems of time-fractional partial differential equations. VooPINN demonstrates that variable-order fractional operators can be embedded directly into a neural-network framework, but it relies on a deep multilayer architecture and does not address singular ordinary differential equations with a variable-order fractional derivative, which is the subject of the present work. However, classical neural-network architectures with multilayer perceptrons are not always optimal for problems with singularities or a specific domain of definition. In addition, training deep networks requires significant computational resources and can be unstable in the presence of singularities in the solution.

In this paper, we propose a combined approach that uses the advantages of both orthogonal polynomials and neural-network optimization methods. We develop a neural-network method (LaNN—Laguerre Neural Network), where orthogonal Laguerre polynomials act as activation functions [36,37,38,39], and the basis functions are shifted to the point $x_0$, which ensures exact satisfaction of the initial conditions of the Cauchy problem at $x_0$. This choice is not accidental: Laguerre polynomials are defined on a semi-infinite interval and form a complete orthogonal basis, which makes them an ideal tool for expanding the solution of the Cauchy problem [38]. Note that neural-network methods based on Laguerre polynomials have already been applied to solve ordinary differential equations [36], as well as equations of the Lane–Emden–Fowler type [40]. However, these works do not consider variable-order fractional derivatives and do not use basis shifting for exact satisfaction of the initial conditions at $x_0$. Unlike classical physics-informed neural networks (PINNs), which require a deep architecture, the proposed method is a network with a single “weight” layer, which simplifies the training process and interpretation of the results. The single-layer structure and the absence of separate nonlinear activations are a consequence of the fact that the Laguerre basis already provides the necessary nonlinearity. Close in spirit approaches based on wavelet neural networks also demonstrate effectiveness for fractional singular models [41], but the use of Laguerre polynomials provides a more natural representation of the solution on a semi-infinite interval and automatic satisfaction of the initial conditions. It is also worth noting a recent work [42], in which a Laguerre neural network optimized by a genetic algorithm was applied to fractal-fractional Bagley–Torvik equations. That study confirms the efficiency of the Laguerre basis for nonlocal fractional models, but it does not consider variable-order fractional derivatives, singular Lane–Emden-type nonlinearities, or exact satisfaction of Cauchy initial conditions via a shifted basis—the features that distinguish the present LaNN method.

The aim of this work is to develop, programmatically implement, and test a neural-network method based on orthogonal Laguerre polynomials with basis shifting for the numerical solution of the Cauchy problem for a fractional Lane–Emden differential equation of variable order. According to the results obtained in [36], among single-layer polynomial networks, LaNN occupies an intermediate position: it is more accurate than a network based on Chebyshev polynomials and faster than a network based on Legendre polynomials. For the class of singular fractional problems under consideration, where a combination of acceptable accuracy and computational efficiency is important, such a compromise seems justified. In addition, unlike ChNN, the Laguerre basis is defined on a half-axis, which agrees better with the nonlocality of the fractional operator. The paper describes the architecture of the LaNN network, the training algorithm based on the Adam method, and the loss function taking into account the equation residual, and provides a software implementation in the Python programming language, including the possibility of parallelization. The influence of the variable order of the fractional derivative on the complexity of the neural network is analyzed.

The paper is organized as follows. In the first section, preliminary information about orthogonal Laguerre polynomials and variable-order fractional calculus is presented. In the second section, the statement of the Cauchy problem for the fractional Lane–Emden equation is formulated. In the third section, the architecture of the LaNN network with a shifted basis, the loss function, and the training algorithm are described. In the fourth section, the software implementation of the method is considered, including vectorized and loop variants of computations. In the fifth section, the results of computational experiments are presented, and the accuracy, stability, and computational efficiency of the proposed approach are analyzed. In the conclusion, the main results are summarized and directions for further work are discussed.

2. Preliminary Information

Let us introduce some preliminary information that we will rely on in this paper.

2.1. Orthogonal Laguerre Polynomials

Laguerre polynomials are classical orthogonal polynomials widely used in mathematics, physics, and engineering. They are solutions of the Laguerre differential equation and form a complete orthogonal system on a semi-infinite interval [37,38].

Definition 1. Laguerre polynomials, denoted as $L_n\left(x\right)$, are defined as solutions of the second-order linear ordinary differential equation (Laguerre equation) [38]:

$$x{y}^{\prime \prime }+(1-x)y^{\prime }+ny=0,$$

where n is a nonnegative integer. They can be defined through the Rodrigues formula [38], which is useful for their generation:

$$L_n\left(x\right)={\displaystyle \frac{e^x}{n!}}{\displaystyle \frac{d^n}{d{x}^n}}\left(e^{-x}{x}^n\right),$$

and also represented as a closed sum [38]:

$$L_n\left(x\right)={\displaystyle \sum _{j=0}^n}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right){\displaystyle \frac{{(-1)}^j}{j!}}{x}^j.$$

Main properties [38].

1. Orthogonality and completeness. A key property of Laguerre polynomials is their orthogonality on the interval $[0,\infty )$ with weight function $e^{-x}$ [38]:

$${\int }_0^{\infty }{e}^{-x}{L}_n\left(x\right)L_m\left(x\right)dx=0,n\ne m,$$

as well as completeness, which means that any function $f\left(x\right)$ with finite norm ${\int }_0^{\infty }{f}^2\left(x\right)e^{-x}dx<\infty$ can be represented as a convergent series in Laguerre polynomials. These properties make them very useful for expanding functions into series.

2. Recurrence relations. Laguerre polynomials satisfy several recurrence relations that simplify their computation [38]:

$$(n+1)L_{n+1}\left(x\right)=(2n+1-x)L_n\left(x\right)-n{L}_{n-1}\left(x\right),L_0\left(x\right)=1,L_1\left(x\right)=1-x.$$

3. Relation with the derivative [38]:

$$x{L}_n^{\prime }\left(x\right)=n{L}_n\left(x\right)-n{L}_{n-1}\left(x\right).$$

4. Generating function [38]

$${\displaystyle \sum _{n=0}^{\infty }}{L}_n\left(x\right)t^n={\displaystyle \frac{1}{1-t}}{e}^{-xt/(1-t)}$$

allows one to obtain all polynomials simultaneously.

2.2. Elements of Fractional Calculus

Fractional calculus can be studied in more detail in the monographs [1,2]. We will focus on some definitions and properties of Riemann–Liouville and Gerasimov–Caputo variable-order fractional derivatives [19,20,21,22,23,24,25].

Definition 2. Let the function $f\left(x\right)$ be integrable on the interval $[0,1]$, and let the function $\alpha \left(x\right)\in C[0,1]$ be sign-positive. Then the Riemann–Liouville fractional integral of order $\alpha \left(x\right)$ is defined as follows:

$${}^{RL}I_{0x}^{\alpha \left(x\right)}f\left(x\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right(x\left)\right)}}{\int }_0^x{(x-\tau )}^{\alpha \left(x\right)-1}f\left(\tau \right)d\tau ,x>0,$$

where $\Gamma \left(\alpha \right(x\left)\right)$ is the Euler gamma function.

Definition 3. Let m be an integer such that $m-1<\alpha \left(x\right)<m$. Then the Riemann–Liouville fractional derivative of order $\alpha \left(x\right)$ is defined as follows:

$${}^{RL}D_{0x}^{\alpha \left(x\right)}f\left(x\right)={\displaystyle \frac{1}{\Gamma (m-\alpha (x\left)\right)}}{\displaystyle \frac{d^m}{d{x}^m}}{\int }_0^x{(x-\tau )}^{m-\alpha \left(x\right)-1}f\left(\tau \right)d\tau ,$$

(1)

where ${\displaystyle \frac{d^m}{d{x}^m}}$ is the ordinary integer-order derivative of order m.

Definition 4. The Gerasimov–Caputo fractional derivative of order $m-1<\alpha \left(x\right)<m$ is defined as follows:

$${}^{GC}D_{0x}^{\alpha \left(x\right)}f\left(x\right)={\displaystyle \frac{1}{\Gamma (m-\alpha (x\left)\right)}}{\int }_0^x{(x-\tau )}^{m-\alpha \left(x\right)-1}{f}^{\left(m\right)}\left(\tau \right)d\tau ,$$

(2)

where $f^{\left(m\right)}\left(x\right)={\displaystyle \frac{d^{m}f\left(x\right)}{d{x}^m}}$ is the ordinary integer-order derivative of order m.

Note an important property of the fractional derivative (2):

$${}^{GC}D_{0x}^{\alpha \left(x\right)}{x}^{\mu }=\left\{\begin{array}{cc}{\displaystyle \frac{\Gamma (\mu +1)}{\Gamma (\mu -\alpha (x)+1)}}{x}^{\mu -\alpha \left(x\right)},\hfill & \mu >m-1,\hfill \\ 0,\hfill & \mu \in \{0,\dots ,m-1\}.\hfill \end{array}\right.$$

(3)

From property (3) it follows that ${}^{GC}D_{0x}^{\alpha \left(x\right)}{x}^{\mu }=0$ for $\mu \in \{0,\dots ,m-1\}$. Therefore, in various applications, the Gerasimov–Caputo fractional derivative is usually chosen, since classical initial conditions are specified for it, and its derivative of a constant gives zero.

Remark on the relation between the Riemann–Liouville and Gerasimov–Caputo derivatives. The well-known relation

$${}^{GC}D_{0x}^{\alpha }f\left(x\right)={}^{RL}D_{0x}^{\alpha }f\left(x\right)-{\displaystyle \sum _{k=0}^{m-1}}{\displaystyle \frac{f^{\left(k\right)}\left(0\right)x^{k-\alpha }}{\Gamma (k-\alpha +1)}},$$

(4)

is valid for a constant order $\alpha$. In the case of a variable order $\alpha \left(x\right)$, it does not hold in general, since $\alpha \left(x\right)$ enters both the kernel and the gamma function, and the outer derivative $d^m/d{x}^m$ also acts on $\alpha \left(x\right)$. Therefore, in this paper we use definition (2) directly and rely on property (3), which remains valid for a variable order as well, since $\alpha \left(x\right)$ is constant with respect to the integration variable $\tau$.

In what follows, we will call differential equations with fractional derivatives (1) or (2) fractional differential equations.

3. Problem Statement

We consider the Cauchy problem for the fractional Lane–Emden equation:

$$\begin{array}{c}{}^{GC}D_{0x}^{\alpha \left(x\right)}y\left(x\right)+{\displaystyle \frac{2}{x}}{y}^{\prime }\left(x\right)+4(2e^{y\left(x\right)}+e^{y\left(x\right)/2})=f\left(x\right),\hfill \\ \hfill y\left(x_0\right)=y_0,y^{\prime }\left(x_0\right)=y_0^{\prime },x\in [x_0,1],\end{array}$$

(5)

where $y\left(x\right)\in C^2[x_0,1]$, the Gerasimov–Caputo variable-order fractional derivative $1<\alpha \left(x\right)<2$, $\alpha \left(x\right)\in C[x_0,1]$ has the form according to (2):

$${}^{GC}D_{0x}^{\alpha \left(x\right)}y\left(x\right)={\displaystyle \frac{1}{\Gamma (2-\alpha (x\left)\right)}}{\int }_0^x{(x-\tau )}^{1-\alpha \left(x\right)}{y}^{\prime \prime }\left(\tau \right)d\tau .$$

(6)

Note that in the case when $f\left(x\right)=0$, $\alpha \left(x\right)=2$, we obtain the Cauchy problem for the classical Lane–Emden equation with exponential nonlinearity [3,4,8,9].

Remark on the lower limit of the fractional derivative. When stating the Cauchy problem for an equation with a fractional derivative ${}^{GC}D_{0x}^{\alpha \left(x\right)}$, the lower limit of the fractional operator plays an important role. If the initial conditions are specified at a point $x_0>0$, then the operator ${}^{GC}D_{0x}^{\alpha \left(x\right)}$ for $x>x_0$ uses the values of $y^{\prime \prime }\left(\tau \right)$ on the entire interval $[0,x]$, including the segment $[0,x_0]$, which is not described by the initial conditions. In this case, the Cauchy problem at $x_0$ for an operator with lower limit 0 is not well-posed without specifying the prehistory of the solution on $[0,x_0]$.

In this paper, two scenarios are considered.

First scenario (original problem). In Section 6, for $f\left(x\right)=0$, we set $x_0=0$. The initial conditions are specified at $x=0$, and the lower limit of the fractional derivative coincides with the point where the initial conditions are specified. There is no prehistory, and the statement is well-posed.

Second scenario (test problem with a known exact solution). When verifying the method, the manufactured-solution technique is used: the equation is considered on the entire interval $[0,1]$ with a known exact solution $y_{ex}\left(x\right)$ that satisfies the operator ${}^{GC}D_{0x}^{\alpha \left(x\right)}$. The initial conditions at $x_0>0$ are taken from this exact solution, so that the prehistory on $[0,x_0]$ is implicitly specified by $y_{ex}\left(x\right)$ itself. This is a standard technique for testing numerical methods, and no contradiction with the statement arises.

It should be emphasized that for the neural-network approximation $y_{LaNN}\left(x\right)$, coincidence with $y_{ex}\left(x\right)$ on $[0,x_0]$ is not guaranteed: the basis functions are defined for all $x\ge 0$ through the analytic continuation ${\psi }_n(x-x_0)$, but the equation on $[0,x_0]$ is not included in the loss function. Therefore, examples with $x_0>0$ in Section 6 should be interpreted precisely as manufactured-solution tests, where the prehistory is known from the exact solution, and not as a demonstration of a well-posed Cauchy problem with arbitrary $x_0>0$. Practical problems with $x_0>0$ and unknown prehistory require a transition to the operator (7); the modification of the method in this direction was not carried out and is a subject of further work.

For practical problems in which the initial conditions are specified at a point $x_0>0$ and the prehistory on $[0,x_0]$ is unknown, one should use the fractional derivative with lower limit $x_0$:

$${}^{GC}D_{x_0,x}^{\alpha \left(x\right)}y\left(x\right)={\displaystyle \frac{1}{\Gamma (2-\alpha (x\left)\right)}}{\int }_{x_0}^x{(x-\tau )}^{1-\alpha \left(x\right)}{y}^{\prime \prime }\left(\tau \right)d\tau .$$

(7)

Such an operator takes into account only the part of the solution that is determined by the initial conditions at $x_0$ and leads to a well-posed Cauchy problem. The transition from ${}^{GC}D_{0x}^{\alpha \left(x\right)}$ to ${}^{GC}D_{x_0,x}^{\alpha \left(x\right)}$ in the proposed LaNN method reduces to replacing $x^{k-\alpha \left(x\right)}$ with ${(x-x_0)}^{k-\alpha \left(x\right)}$ in the formulas for the fractional derivative of the basis functions. This can be a subject of a separate study; in this paper, the main emphasis is on the case $x_0=0$ (original problem) and on the verification of the method using the manufactured solution (test examples).

The classical Lane–Emden equation with exponential nonlinearity arises in the description of hydrostatic equilibrium of spherically symmetric polytropic stars, isothermal gas spheres, and gravitationally bound clouds. In this context, x is a dimensionless radial coordinate, and $y\left(x\right)$ is a dimensionless Emden function, which can be interpreted as the logarithm of density, temperature, or as a dimensionless gravitational potential. The term ${\displaystyle \frac{2}{x}}{y}^{\prime }\left(x\right)$ reflects the spherical geometry of the Laplace operator: it appears in the transition from Cartesian coordinates to radially symmetric ones. The nonlinear term $4(2e^{y\left(x\right)}+e^{y\left(x\right)/2})$ models an effective source/self-action that depends exponentially on the state of the medium. If y is interpreted as the logarithm of dimensionless density, then $e^y$ is proportional to density, and $e^{y/2}$ is a weaker nonlinear factor that may correspond, for example, to heat release, reaction, or nonlinear coupling with the gravitational field.

The Gerasimov–Caputo variable-order fractional derivative ${}^{GC}D_{0x}^{\alpha \left(x\right)}y\left(x\right)$ in the model equation (5) generalizes the second derivative and introduces the effect of spatial nonlocality into the model. Unlike the classical second derivative, which is determined only by the local behavior of the function, the fractional operator takes into account the prehistory or distribution of the solution on the entire interval $[0,x]$. The variable order $\alpha \left(x\right)$ means that the degree of nonlocality changes with radius: in internal regions and on the periphery, the memory contribution may be different. This allows one to describe media whose rheological, thermal, or gravitational properties depend on the local state. The right-hand side $f\left(x\right)$ plays the role of an effective external source or sink: it can be a mass distribution, energy release, or an external field. The initial conditions $y\left(x_0\right)=y_0$, $y^{\prime }\left(x_0\right)=y_0^{\prime }$ specify the value of the function and its gradient at the starting radius.

4. LaNN Architecture and Its Description

The neural-network method allows one to obtain the solution of problem (5) in the form:

$$y_{LaNN}\left(x\right)=y_0+y_0^{\prime }(x-x_0)+{\displaystyle \sum _{n=0}^N}{w}_n{\psi }_n(x-x_0),$$

(8)

where ${\psi }_n\left(\xi \right)={\xi }^2{L}_n\left(\xi \right)={\sum }_{j=0}^n{l}_{n,j}{\xi }^{j+2},\xi =x-x_0$ are modified Laguerre polynomials with a shift to the point $x_0$ that satisfy the initial conditions (5):

$$y_{LaNN}\left(x_0\right)=y_0,y_{LaNN}^{\prime }\left(x_0\right)=y_0^{\prime },$$

$l_{n,j}=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{j}}\right){\displaystyle \frac{{(-1)}^j}{j!}}$ are binomial coefficients, N is the order of the Laguerre polynomial, and $w_n$ are weight coefficients that are refined during training. The basis shift to the point $x_0$ ensures exact satisfaction of the initial conditions of the Cauchy problem for any $x_0>0$, which is fundamentally important for singular problems where $x_0$ may be different from zero.

It should be emphasized that the basis shift to the point $x_0$ ensures exact satisfaction of the initial conditions but does not change the lower limit of the fractional operator. In this paper, the fractional derivative is taken with the operator ${}^{GC}D_{0x}^{\alpha \left(x\right)}$ with lower limit 0. This corresponds either to the case $x_0=0$ (original problem) or to the case when the solution on $[0,x_0]$ is known and taken into account in the right-hand side (test examples with manufactured solution). For practical problems with $x_0>0$ and unknown prehistory, one should use the operator ${}^{GC}D_{x_0,x}^{\alpha \left(x\right)}$; the corresponding formulas are obtained by replacing $x^{k-\alpha \left(x\right)}$ with ${(x-x_0)}^{k-\alpha \left(x\right)}$.

Training proceeds as follows. A training sample of fixed size M is chosen from $x\in [x_0,1]$. The points of this sample are called collocation points.

Remark 1. Collocation points can be chosen in different ways, for example, uniformly on the entire interval $[x_0,1]$ or subject to some distribution, for example, logarithmic or Chebyshev.

Then, at the chosen collocation points, the residual between the left-hand side of equation (5) with account for (8) and its right-hand side is calculated:

$$L\left(x_i\right)={}^{GC}D_{0x}^{\alpha \left(x\right)}{y}_{LaNN}\left(x_i\right)+{\displaystyle \frac{2}{x_i}}{y}_{LaNN}^{\prime }\left(x_i\right)+4\left(2e^{y_{LaNN}\left(x_i\right)}+e^{y_{LaNN}\left(x_i\right)/2}\right)-f\left(x_i\right),$$

(9)

$i=1,\dots ,M.$ Here, the Gerasimov–Caputo fractional derivative of the shifted basis is calculated taking into account the binomial expansion:

$${}^{GC}D_{0x}^{\alpha \left(x\right)}{y}_{LaNN}\left(x\right)={\displaystyle \sum _{n=0}^N}{w}_n{\displaystyle \sum _{j=0}^n}{l}_{n,j}{\displaystyle \sum _{k=2}^{j+2}}\left({\displaystyle \genfrac{}{}{0pt}{}{j+2}{k}}\right){(-x_0)}^{j+2-k}{\displaystyle \frac{\Gamma (k+1)}{\Gamma (k+1-\alpha (x\left)\right)}}{x}^{k-\alpha \left(x\right)}.$$

Note that the linear terms $y_0+y_0^{\prime }(x-x_0)$ give a zero contribution to the Gerasimov–Caputo fractional derivative of order $1<\alpha \left(x\right)<2$, since ${}^{GC}D_{0x}^{\alpha \left(x\right)}$ for $\mu \in \{0,1\}$.

Then we form the loss function $Loss\left(w\right)$ as the normalized weighted mean-square value of the residual at the collocation points (9) with the addition of $L_2$-regularization of the weights $w_n$ of solution (8):

$$Loss\left(w\right)={\displaystyle \frac{1}{M}}{\displaystyle \sum _{i=1}^M}{\displaystyle \frac{\phi \left(x_i\right)}{\overline{\phi }}}·L{\left(x_i\right)}^2+\lambda {\displaystyle \sum _{n=0}^N}{w}_n^2,$$

(10)

where $\lambda$ is the regularization parameter, $\phi \left(x_i\right)$ are the weights of the collocation points, and $\overline{\phi }$ is the mean value.

Remark 2. It should be noted that the weight coefficients $\phi \left(x_i\right)$ for the collocation points can be chosen to improve accuracy in the region of small values, where the solution has a singularity.

Remark 3. $L_2$-regularization allows one to avoid overfitting and excessively large weight coefficients $w_n$ of solution (8). The neural network can adapt to noise or specific features of a particular set of collocation points, which leads to nonsmooth or nonphysical solutions. $L_2$-regularization limits the growth of weight coefficients, promoting a smoother approximation. In addition, with a large number of parameters (high polynomial degree), the weight coefficients can become very large, causing numerical instability (Runge’s phenomenon). Regularization mitigates this effect. The loss function shows how strongly the equation is violated at the collocation points (10). The smaller the loss, the better the network “knows” the solution. To achieve this, it is necessary to reduce this loss by updating (choosing) the weight coefficients at each iteration.

Remark 4. The process of constructing the loss function (10) is called the forward pass. It forms a computational graph that stores the sequence of operations and is necessary for the subsequent calculation of derivatives (gradients).

The weights are updated in accordance with the Adam (Adaptive Moment Estimation) algorithm from the family of gradient descent methods.

The formulas for the weights $w=(w_0,\dots ,w_N)$ at iteration k are as follows:

1. Gradient

$$g_k\left(w\right)={\nabla }_{w}L\left(w_{k-1}\right),$$

(11)

where ${\nabla }_w=\left|{\displaystyle \frac{\partial }{\partial w_0}},{\displaystyle \frac{\partial }{\partial w_1}},\dots ,{\displaystyle \frac{\partial }{\partial w_N}}\right|$, L is the loss function (10).

Remark 5. It should be noted that usually the gradients (11) are calculated using the backpropagation algorithm. The backward pass through the computational graph from the loss function node to the model weights is performed using the chain rule. Therefore, the process of constructing gradients is called the backward pass.

2. Biased estimates of the first and second moments:

$$m_k={\beta }_1{m}_{k-1}+(1-{\beta }_1)g_k,v_k={\beta }_2{v}_{k-1}+(1-{\beta }_2)g_k^2,$$

the parameter values ${\beta }_1=0.9$, ${\beta }_2=0.999$ are taken by default.

3. Bias correction:

$${\widehat{m}}_k={\displaystyle \frac{m_k}{1-{\beta }_1^k}},{\widehat{v}}_k={\displaystyle \frac{v_k}{1-{\beta }_2^k}}.$$

4. Weight update:

$$w_k=w_{k-1}-\eta {\displaystyle \frac{{\widehat{m}}_k}{\sqrt{{\widehat{v}}_k}+ϵ}},$$

where $\eta$ is the initial learning rate, $ϵ$ is a small constant to avoid division by zero (usually ${10}^{-8}$).

One training cycle (epoch) consists of one forward and one backward pass over all collocation points, as well as an optimizer step (Adam), i.e., updating the weights to reduce the loss function. For each training epoch, after updating the weights, it is necessary to zero the gradients to prevent their accumulation.

As a result, after a given number of training epochs, solution (8) is output with an optimal set of weight coefficients $w_n$. The architecture of the network with a shifted basis is shown in Figure 1.

Here, the input x is transformed into the shifted variable $\xi =x-x_0$. Then the basis functions ${\psi }_n\left(\xi \right)={\xi }^2{L}_n\left(\xi \right)$ are multiplied by the trainable weights $w_n$. The term $y_0+y_0^{\prime }(x-x_0)$ ensures exact satisfaction of the initial conditions. The results are summed and fed to the output $y_{LaNN}\left(x\right)$.

5. Implementation of the LaNN Architecture

The LaNN neural-network architecture was implemented in the Python programming language [43] using the PyTorch library [44].

The main implementation stages are as follows:

1. Initialization of parameters—input of the derivative order $\alpha$, the polynomial degree N, the number of collocation points, training and regularization parameters, choice of the weight function, early-stopping settings, and the coordinate of the initial point $x_0$. The early-stopping condition can be, for example, the absence of improvement by at least 0.1% from the best value over a number of epochs, after which training stops. This can save computational time. The regularization parameter is ${10}^{-5}$ by default.

2. Construction of the neural-network model—the LaguerreNN class, which implements the approximation of the solution using shifted Laguerre polynomials and automatically takes into account the initial conditions (5) at the point $x_0$. The LaguerreNN class implements the fractional derivative with the operator ${}^{GC}D_{0x}^{\alpha \left(x\right)}$. The x0 parameter specifies the basis shift point and the point where the initial conditions are specified. For x0 = 0, the statement of the Cauchy problem is well-posed without additional assumptions. For x0 > 0, the method should be applied either in the manufactured-solution mode (when the prehistory is known from the exact solution) or after replacing the operator with ${}^{GC}D_{x_0,x}^{\alpha \left(x\right)}$.

3. Formation of the training set—a set of points on the interval. Three distributions are provided: uniform using torch.linspace; logarithmic using torch.logspace, which gives more points near zero, where the solution and its derivative change rapidly; and the Chebyshev distribution, which gives clustering of points near the boundaries of the interval. The number of points is specified by the user (by default 10,000).

4. Training of the model—minimization of the special loss function (10), including the residual (9), the weight function, and regularization. The weight function allows one to focus on regions with large errors (for example, near zero). The user can choose an arbitrary expression, for example $\phi \left(x\right)=x$. If the user does not choose weights, this is equivalent to the uniform case, where $\phi \left(x\right)=1$. Before use, the weights are normalized so that their mean value equals 1 (this does not affect the optimum but stabilizes the gradient scale).

5. Optimizer—Adam is used with a given initial learning rate. To automatically reduce the step on the loss-function plateau, the ReduceLROnPlateau scheduler is used with a reduction factor of 0.5 and a patience of 500 epochs. The scheduler controls the learning rate $\eta$, which can decrease if the loss function stops decreasing over a given number of epochs, i.e., the learning rate is adaptive.

6. Evaluation of results—comparison with the exact solution on a test grid, output of a table and graphs. The maximum and mean absolute errors are calculated, and a comparison table at the nodes x is constructed. Visualization includes four graphs: comparison of the exact and approximate solutions; error on a logarithmic scale; loss function; and the function $\alpha \left(x\right)$.

7. Computation modes. When solving the Cauchy problem (5) for fractional differential equations using the neural-network approach, the computational costs can be significant, especially when using a large number of collocation points (up to 200,000), high polynomial degrees (up to 30), and a large number of training epochs (up to 10,000). To accelerate computations in the developed method, parallelization algorithms implemented using the PyTorch library and the hardware capabilities of modern computing systems are used.

Two variants of the forward pass of the network are implemented:

• Loop (forward_loop): sequential summation over indices using Python loops. A slow mode, but more understandable.
• Vectorized (forward_vectorized): all operations are performed using PyTorch tensor operations on the CPU. At the algorithmic level, the key technique here is replacing loops over polynomial degrees (from 0 to degree) with tensor operations. In the forward_vectorized variant, all basis functions, their derivatives, and fractional derivatives are calculated using broadcasting and matrix multiplications. This allows efficient use of CPU SIMD instructions. As shown in Section 5, the vectorized version works on average 30 times faster than the loop version, and the use of vectorization does not require changing the network architecture and is completely equivalent in output values.

Some computational aspects: all tensors and the model are moved to the device specified in the device variable (CPU); for numerical stability, the gamma function is calculated using torch.lgamma and exponentiation to avoid overflow for large argument values; the polynomial degree is limited to 30 to prevent instability and overflow.

6. Simulation Results

For further calculations, a personal computer with the following characteristics was used: operating system Windows 11; processor 12th Gen Intel(R) Core(TM) i9-12900KF; logical cores: 24; physical cores: 16; current frequency: 3200 MHz; maximum frequency: 3200 MHz; current CPU load: 2.3%; RAM: 128 GB; graphics card NVIDIA GeForce RTX 4060 Ti. The calculations were performed without using CPU multithreading or CUDA cores.

Remark 6. It should be noted that the selection of hyperparameters (number of collocation points, number of epochs, and regularization coefficient) was carried out manually based on a series of preliminary experiments. The criteria were: training stability, absence of solution oscillations, absence of overfitting, accuracy on the test grid, and acceptable computation time. The absence of systematic automatic hyperparameter selection is a limitation of the work. A more rigorous approach, for example Bayesian optimization or grid search taking into account training time, can be implemented in the future.

Remark 7. (Remark on the statement of initial conditions.) In all examples of this section where the right-hand side $f\left(x\right)=0$ (Examples 3 and 4), the initial conditions are specified at the point $x_0=0$:

$$y\left(0\right)=y_0,y^{\prime }\left(0\right)=y_0^{\prime }.$$

The computational domain begins at $x_{min}>0$, which allows avoiding the singularity of the term ${\displaystyle \frac{2}{x}}{y}^{\prime }\left(x\right)$ at zero. The shift of the Laguerre basis is performed to the point $x_0=0$, which automatically ensures satisfaction of the initial conditions. In test examples 1 and 2, the initial conditions are specified at the point $x_0>0$ and are taken from the exact solution; the prehistory on $[0,x_0]$ is known in this case, and the statement of the Cauchy problem is well-posed.

Remark 8. (On collocation points and the interval $[0,x_{min}]$.) In all examples, the collocation points are chosen on the interval $[x_{min},x_{max}]$, where $x_{min}>0$. The contribution of the segment $[0,x_{min}]$ to the fractional derivative is taken into account analytically through the global LaNN ansatz: the basis functions ${\psi }_n(x-x_0)$ are defined for all $x\ge 0$, and their fractional derivative is calculated by the formula given in Section 3. The initial conditions at zero fix the linear terms. At the same time, the equation on the interval $[0,x_{min}]$ is not included in the loss function, and the behavior of the solution there is controlled only by the analytic continuation of the ansatz. This may introduce additional error that decreases with decreasing $x_{min}$. The study of the sensitivity of the method to the choice of $x_{min}$ was not carried out in this paper and is a direction for further research.

6.1. Test Example 1

Consider the Cauchy problem where equation (5) is specified with a special right-hand side where $\alpha \left(x\right)=(2-0.03x)x^3$:

$${}^{GC}D_{0x}^{\alpha \left(x\right)}y\left(x\right)+{\displaystyle \frac{2}{x}}{y}^{\prime }\left(x\right)+4(2e^{y\left(x\right)}+e^{y\left(x\right)/2})=f\left(x\right),$$

(12)

$$f\left(x\right)={\displaystyle \frac{12x^{3-\alpha \left(x\right)}}{\Gamma (4-\alpha (x\left)\right)}}-{\displaystyle \frac{0.72x^{4-\alpha \left(x\right)}}{\Gamma (5-\alpha (x\left)\right)}}+12x-0.24x^2+4\left(2e^{\alpha \left(x\right)}+e^{\alpha \left(x\right)/2}\right),$$

Thanks to the basis shift to the point $x_0=0.01$, the initial conditions are satisfied exactly:

$$y\left(0.01\right)=0.0000019997,y^{\prime }\left(0.01\right)=0.0005998799,x\in [0.01,1].$$

The solution of problem (12) is the function

$$y_{ex}\left(x\right)=\alpha \left(x\right)x^3=(2-0.03x)x^3.$$

(13)

The hyperparameters for LaNN training are given in Table 1.

The collocation points are chosen uniformly distributed on $[0.01,1]$. Training was completed in 13.42 s, best Loss: $1.49\times {10}^{-3}$. The LaNN training results are given in Table 2.

At the table nodes, the error varies from $4.547474\times {10}^{-13}$ at $x=0.01$ to $4.823208\times {10}^{-4}$ at $x=0.89$; at $x=1$, the error is $1.192093\times {10}^{-7}$. The largest absolute errors are concentrated in the middle and right parts of the interval. Near zero, the absolute error is small due to the smallness of the solution itself, but the relative error there can be noticeable. This is expected for a singular fractional Lane–Emden equation with the point $x_0=0.01$ and variable order $\alpha \left(x\right)$. The weight $\phi \left(x\right)=x$ shifts the emphasis to the right part. Note that the initial conditions at $x_0=0.01$ are satisfied exactly due to the basis shift.

It can be seen (Table 3) that MAE and RMSE are close, which indicates the absence of large outliers. Figure 2 shows graphs of the exact (13) and neural-network solutions, the loss function, and the absolute error.

Note that Loss decreases from $\sim {10}^3$ to $\sim 1.5\times {10}^{-3}$. After about 8000 epochs, the improvement slows down, so early stopping is possible here to save time.

To verify the stability of the calculations with respect to the initialization of the LaNN weights for Example 1, a computational experiment with 10 independent training runs was conducted. The reproducibility parameters of the computational experiment are given in Table 4.

The experiment was conducted without early stopping, so $\delta$ and patience are not defined. The results of the experiment are given in Table 5.

The relative spread of accuracy metrics (ME, MAE, RMSE) does not exceed 9.2%, and that of training time is 2.5% (Table 6).

The largest spread is observed for the internal metric Best Loss (11.0%), which is expected, since it reflects the optimization trajectory rather than the quality of the final solution. All 10 runs produced a correct solution without oscillations or outliers (Figure 3).

It should also be noted that a comparison was made with a neural network based on the Legendre basis under identical conditions (settings). Five independent runs (seeds) were performed, and the averaged metrics were calculated and are presented in Table 7.

The Laguerre-basis neural network demonstrated superior performance here: it outperformed the others by approximately 2.2 times in terms of BestLoss, 2.9 times in RMSE, 3.0 times in ME, and 2.5 times in MAE; regarding ±Std, the Laguerre network was more stable, showing less variation between runs; and there was virtually no difference in training time.

This example demonstrates that LaNN can yield better results than other polynomial-based neural networks.

Remark 9. Equation (12) has a singularity in the term ${\displaystyle \frac{2}{x}}{y}^{\prime }\left(x\right)$ near the boundary point $x=0$; therefore, the numerical solution obtained on the basis of LaNN can differ strongly due to large values of the derivative, which leads to a large deviation from the true solution and complicates the physical interpretation of the process. Therefore, to improve accuracy, it is necessary to cluster points near the boundary point $x=0$. This can be achieved by choosing the distribution of collocation points according to the Chebyshev law.

6.2. Test Example 2 (Stress Test)

Consider the Cauchy problem (12). Choose the interval $[{10}^{-10},1]$; here $x_0={10}^{-10}$, so $y\left({10}^{-10}\right)\approx 2·{10}^{-30},y^{\prime }\left({10}^{-10}\right)\approx 6·{10}^{-20}$. Also choose the weight $\phi \left(x\right)=1/x$, which will shift the emphasis to the left part of the interval under consideration, i.e., clustering of points near the singularity $x=0$ will occur. The values of the remaining parameters for LaNN training are taken from Table 1, and the collocation points are still uniformly distributed. The training results are given in Table 8.

We see (Table 8) that near the singularity $x=0$, a significant deviation from the exact solution (13) occurs, negative values appear in the numerical solution, and the absolute error accumulates (Figure 4). We observe signs of poor conditioning of the method.

If the collocation points are chosen according to the Chebyshev law, which produces clustering of points near the boundaries of the interval $[{10}^{-10},1]$, we obtain a better result (Table 9).

From Table 9, we see that the error near the boundaries is significantly lower than the error in Table 8, and all values of the numerical solution are positive. Note that the maximum error in this example reaches $\sim 2.7\times {10}^{-2}$ (at $x=1$), which is higher than in Example 1 and is explained by the greater singularity of the problem on the interval $[{10}^{-10},1]$.

The LaNN training times for the vectorized version ($t_{vec}$), built on tensor operations, and the loop version ($t_{cic}$), based on classical loops, are given in Table 10. The hyperparameter values were taken from Table 1. The collocation points were chosen uniformly and according to Chebyshev.

Here one can see that the vectorized version trains on average 30 times faster than the loop version (the ratio $t_{cic}/t_{vec}$ ranges from 27.4 to 30.6). We see that the vectorized implementation in PyTorch based on tensor operations can give a speedup of about 30 times for this example. One can try to further accelerate LaNN training by implementing a parallel algorithm for CPU or GPU, but this is a topic for future research.

6.3. Test Example 3

Consider problem (5) with $f\left(x\right)=0$, $y\left(0\right)=1$, $y^{\prime }\left(0\right)=0$. Parameters for LaNN training: $\alpha \left(x\right)=2-0.03x$; number of epochs 8000; initial learning rate $\eta =0.005$; $L_2$-regularization coefficient $\lambda ={10}^{-5}$; polynomial degree $N=15$; number of collocation points $M=20$, uniformly distributed on $[0.01,1]$; weights for the residual: ${\phi }_1\left(x\right)=1$, ${\phi }_2\left(x\right)=x$, ${\phi }_3\left(x\right)=1/x$.

Since in this case we do not know the exact solution of problem (5) as in Example 1, we need to find a numerical solution. For this purpose, we use an implicit nonlocal finite-difference scheme of the first order of accuracy based on the Grünwald–Letnikov approximation $y_{GL}$ with the Gerasimov–Caputo correction [45]. The calculations were performed on a uniform grid with step $h=(x_{max}-x_{min})/N$, $N=500$. The nonlinear equation at each step was solved by Newton’s method with accuracy ${10}^{-12}$. The scheme has the first order of approximation in h. A rigorous study of the grid convergence of the reference solution $y_{GL}$ was not carried out in this paper; indirect confirmation of its suitability is its agreement with the analytical solution in test examples 1 and 2, where the exact solution is known.

The calculated data are given in Table 11.

Table 11 shows how the choice of the weight function $\phi \left(x\right)$ in the loss function affects the accuracy of the LaNN neural-network solution when solving the fractional Lane–Emden equation with $f\left(x\right)=0$ and variable order $\alpha \left(x\right)=2-0.03x$. The numerical solution $y_{GL}$ obtained by the Grünwald–Letnikov finite-difference scheme with the Gerasimov–Caputo correction is used as a reference. The table compares three weight variants: ${\phi }_1=1$, ${\phi }_2=x$, ${\phi }_3=1/x$. The columns $Er{r}_1,Er{r}_2,Er{r}_3$ are the absolute errors $|y_{GL}-y_{LaNN}|$ for the corresponding weights.

At the left boundary point $x=0.01$, all three LaNN variants give values practically coinciding with the reference $y_{GL}$. The errors are $5.50\times {10}^{-7}$ for ${\phi }_1=1$, $1.09\times {10}^{-5}$ for ${\phi }_2=x$, and $1.64\times {10}^{-7}$ for ${\phi }_3=1/x$. This confirms that correct account of the initial conditions $y\left(0\right)=1$, $y^{\prime }\left(0\right)=0$ does not introduce additional error at the boundary of the computational domain.

The reference solution $y_{GL}$ decreases from $\approx 0.9995$ at $x=0.01$ to $\approx -1.3390$ at $x=1.00$, changing sign in the middle part of the interval. All three neural-network solutions $y_{LaNN}$ qualitatively reproduce this behavior without giving nonphysical oscillations or outliers. The differences between the variants with different $\phi$ are small compared to the amplitude of the solution.

The error is mainly of order ${10}^{-3}$–${10}^{-4}$. The most noticeable deviations are concentrated in the region $x\sim 0.12$–$0.45$, where the errors reach $\left(3--6\right)\times {10}^{-3}$. Closer to the right boundary, the errors usually decrease to $\left(0.8--1.3\right)\times {10}^{-3}$. This indicates acceptable accuracy of the method for this example, where the exact analytical solution is unknown and the reference is obtained numerically. Note that the errors given in Table 11 at individual nodes (up to $6\times {10}^{-3}$) do not contradict the integral metrics $ME=2.37\times {10}^{-3}$, $MAE=1.56\times {10}^{-3}$, $RMSE=1.63\times {10}^{-3}$ for the case ${\phi }_1=1$: the integral metrics were calculated on a denser test grid of 2000 points, while the table contains only 10 characteristic nodes, and in the neighborhood of $x\sim 0.12$–$0.23$ the local error may exceed the mean.

The weight ${\phi }_1=1$ gives the most uniform error distribution. For ${\phi }_1=1$, the error varies from $5.66\times {10}^{-3}$ at $x=0.12$ to $9.27\times {10}^{-4}$ at $x=1.00$, without sharp local spikes. In terms of integral metrics, this variant has the smallest maximum error: $ME=2.37\times {10}^{-3}$ with $MAE=1.56\times {10}^{-3}$ and $RMSE=1.63\times {10}^{-3}$. That is, the uniform weight gives a balanced approximation.

The weight ${\phi }_2=x$ shifts the emphasis to the right but does not give an explicit gain. For ${\phi }_2=x$, the error in the left section $x=0.12$–$0.23$ turns out to be the largest among the three variants: $6.11\times {10}^{-3}$ and $4.57\times {10}^{-3}$. At the same time, the total loss function for this variant is the highest: $Loss=8.92\times {10}^{-3}$, and $ME=2.84\times {10}^{-3}$. This means that concentration on the right part of the interval does not lead to an improvement in global accuracy for this example.

The weight ${\phi }_3=1/x$ improves training and accuracy near the left boundary. For ${\phi }_3=1/x$, the smallest value of the loss function is achieved: $Loss=1.01\times {10}^{-3}$, which indicates faster and more stable training. Near the initial point, the error is minimal ($1.64\times {10}^{-7}$ at $x=0.01$), and at a number of right points, for example $x=0.67$, $0.78$, $1.00$, the error is also smaller than for the other variants. At the same time, the maximum error $ME=2.43\times {10}^{-3}$ only slightly exceeds the value for ${\phi }_1=1$. Thus, the weight $1/x$ is a successful compromise: it strengthens the left part, where the solution and its derivative change faster, and at the same time does not worsen the overall accuracy.

The choice of the weight function controls not so much the overall accuracy as the distribution of the error. The mean metrics for all three cases are close: $MAE\approx 1.5\times {10}^{-3}$, $RMSE\approx 1.6\times {10}^{-3}$. This means that LaNN is stable with respect to the choice of $\phi \left(x\right)$, but the local accuracy in different parts of the interval can be purposefully redistributed. If uniform accuracy is needed, ${\phi }_1=1$ is better; if behavior near the singularity $x\to 0$ is more important, ${\phi }_3=1/x$ is preferable.

Consider an example with a different function $\alpha \left(x\right)$.

6.4. Test Example 4

Consider problem (5) with $f\left(x\right)=0$ and the same initial conditions as in Example 3: $y\left(0\right)=1$, $y^{\prime }\left(0\right)=0$. Parameters for LaNN training: $\alpha \left(x\right)=1+sin\left(x\right)$, $\phi \left(x\right)=1$, the remaining parameters are taken from Example 3.

According to Table 12 and Figure 5, we see an acceptable training result.

In Table 12, it can be seen that for $x\gtrsim 0.52$, the neural-network solution $y_{LaNN}$ changes sign and becomes negative. It is important to emphasize that this is not an artifact of LaNN training: the reference numerical solution $y_{GL}$ behaves similarly, and the error remains small (ME $\approx 5.62\times {10}^{-3}$, MAE $\approx 2.76\times {10}^{-3}$, RMSE $\approx 3.25\times {10}^{-3}$). Consequently, the negative values reflect a property of the mathematical model itself for the chosen variable order of the fractional derivative.

The physical admissibility of negative y depends on the interpretation of the function $y\left(x\right)$. Since the equation contains exponential nonlinearities $e^y$ and $e^{y/2}$, it is most natural to interpret y not as density or temperature itself but as the logarithm of dimensionless density/temperature or as a dimensionless potential. In this case, a negative y does not mean negative density or temperature: the physical quantity $e^y$ remains strictly positive. For example, at $x=0.56$ we have $y\approx -0.103$, whence $e^y\approx 0.902$; at $x=1.0$, $y\approx -1.171$, $e^y\approx 0.310$. This means that density/temperature or potential become smaller than some reference value but do not lose physical meaning.

The metrics: ME: $5.62\times {10}^{-3}$, MAE: $2.76\times {10}^{-3}$, RMSE: $3.25\times {10}^{-3}$ indicate the absence of outliers in the values. The numerical algorithm $y_{GL}$ worked in 0.08 s.

From Figure 5, we see that the LaNN neural network copes well with finding the solution of the problem: this is visible when comparing with the numerical solution. Also, from the loss function, it can be seen that the number of epochs can be reduced to 7000, which will accelerate training.

To confirm the efficiency of the neural-network algorithm, the inference time was measured—the direct run of an already trained model with fixed weight coefficients on a test sample of 2000 points. For this purpose, 10 runs were performed for the same trained model (Table 13). The goal is to evaluate the stability of the forward-pass time.

From Table 13, it can be seen that the inference time is approximately 1.5–1.9 milliseconds per run. The spread between runs is small, which indicates the stability of the computations. The fastest execution is 1.46 ms, and the slowest is 1.87 ms. These values show that the forward pass of the network takes fractions of a second, which is quite acceptable for many practical problems.

Next, the metrics were calculated (Table 14).

The mean inference time for 2000 points is 1.66 ms, and per point it is 0.83 microseconds (Table 14). This is a very high computation speed. The standard deviation of $\sim 0.12$ ms confirms that the execution time varies weakly, i.e., the algorithm works predictably. The minimum and maximum are also close to the mean, which indicates the absence of outliers.

7. Conclusions

In this work, the LaNN (Laguerre Neural Network) neural-network method for the numerical solution of the Cauchy problem for the fractional Lane–Emden equation of variable order has been developed, programmatically implemented, and tested. The method is based on representing the solution as an expansion in orthogonal Laguerre polynomials with a shift to the point $x_0$, which ensures exact satisfaction of the initial conditions of the Cauchy problem at $x_0$. The weight coefficients are tuned during training by minimizing a special loss function that includes the residual of the differential equation and $L_2$-regularization. The conducted computational experiments have shown the following.

In test examples with a known analytical solution (Examples 1 and 2), the maximum absolute error was $5.6\times {10}^{-4}$ for the interval $[0.01,1]$ and $\sim 2.7\times {10}^{-2}$ for the substantially singular interval $[{10}^{-10},1]$ with the Chebyshev distribution. In examples with a numerical reference (Examples 3 and 4), the maximum absolute error did not exceed $6.11\times {10}^{-3}$. These values confirm the efficiency of LaNN for the class of singular fractional differential equations. The basis shift allows exact satisfaction of the initial conditions for any $x_0$; in examples with $f\left(x\right)=0$, the initial conditions were specified at $x_0=0$, and the shift to $x_0>0$ was applied in test examples with a manufactured solution.

The uniform and Chebyshev distributions give comparable results; however, the choice of distribution can be optimized depending on the behavior of the solution. In particular, for the strongly singular interval $[{10}^{-10},1]$, the Chebyshev distribution made it possible to reduce the maximum error by more than a factor of 60 compared with the uniform distribution (from $\sim 1.6$ to $\sim 2.7\times {10}^{-2}$).

The introduction of weights into the loss function allows one to control local accuracy: weights growing toward the boundary increase accuracy near the singularity. In Example 3, the weight ${\phi }_3=1/x$ gave the smallest value of the loss function ($1.01\times {10}^{-3}$) and the minimum error near $x_0$, while the weight ${\phi }_1=1$ gave the most uniform error distribution over the interval.

A change in the exponent $\alpha$ substantially affects the complexity of training. For $\alpha \left(x\right)=1+sin\left(x\right)$, the value of the loss function turned out to be more than an order of magnitude higher than for $\alpha \left(x\right)=2-0.03x$, which indicates the need to adapt the architecture (for example, increase the polynomial degree or the number of collocation points) to preserve accuracy.

The vectorized implementation in PyTorch provides an average speedup of 30 times compared with the loop implementation and acceptable computation time (units to tens of seconds) for practical problems. The possibility of parallelization on a GPU can give a further increase in the speed of the training algorithm.

It has been shown that after training, the LaNN neural network provides fast and stable inference. This is a key advantage over traditional numerical methods, which may require re-solving the system of equations when points or parameters change. The obtained time characteristics (about 1.66 ms for 2000 points, or 0.83 $\mu$s per point) confirm that the method is suitable for practical use, especially in problems where multiple evaluations of the solution at different points with fixed model parameters are required. In combination with the possibility of parallelization on a GPU, this makes LaNN a promising tool for the numerical analysis of variable-order fractional differential equations, especially in cases where classical methods encounter difficulties due to singularity or nonlocality.

The proposed LaNN implementation uses the operator ${}^{GC}D_{0x}^{\alpha \left(x\right)}$. This is natural for problems with $x_0=0$ and for verification by the manufactured-solution method. For practical problems with initial conditions at a point $x_0>0$ and unknown prehistory on $[0,x_0]$, a transition to the operator ${}^{GC}D_{x_0,x}^{\alpha \left(x\right)}$ is required. Modification of the method in this direction seems promising and will be considered in further work.

Compared with classical finite-difference schemes, the LaNN method does not require grid construction, provides an analytical representation of the solution, and has fast inference (about 1.5–1.9 ms for 2000 points). Compared with physics-informed neural networks (PINNs), LaNN uses a single weight layer and an orthogonal Laguerre basis, which simplifies training and interpretation. Compared with spectral methods based on Legendre or Chebyshev polynomials, LaNN is naturally defined on a half-axis and agrees better with the nonlocality of the fractional operator. In contrast to the recently proposed VooPINN [?], which handles variable-order operators within a deep PINN architecture, LaNN uses a single trainable layer with an orthogonal Laguerre basis, which simplifies training and interpretation. Compared with the GA-optimized Laguerre network for fractal-fractional Bagley–Torvik equations [?], the present method additionally employs a shifted basis for exact enforcement of the Cauchy initial conditions and is adapted to the singular Lane–Emden nonlinearity. At the same time, a direct quantitative comparison of LaNN with PINNs, spectral methods, and high-order finite-difference schemes was not carried out in this paper. This is a limitation and a direction for future research.

The limitations of the work include: (a) manual selection of hyperparameters; (b) absence of a systematic study of the grid convergence of the reference solution $y_{GL}$; (c) absence of a quantitative comparison with other numerical and neural-network methods; (d) absence of an analysis of the sensitivity of the method to the choice of $x_{min}$ for $x_0=0$; (e) limitation of the domain of applicability to the cases $x_0=0$ and manufactured solution for $x_0>0$.

Further research may be aimed at adapting the method to multidimensional problems, using other orthogonal bases (Hermite, Legendre), and integrating with active-learning methods for automatic selection of collocation points, as well as improving parallel training algorithms on CPU or GPU [26]. A separate direction is the transition to the operator ${}^{GC}D_{x_0,x}^{\alpha \left(x\right)}$ for practical problems with $x_0>0$ and unknown prehistory.