Algorithms for Solving Ordinary Differential Equations Based on Orthogonal Polynomial Neural Networks

1. Introduction

In recent years, machine learning methods, particularly artificial neural networks (ANNs), have found widespread application beyond traditional data analysis tasks, entering the realm of scientific computing [1]. One promising direction is the use of ANNs for the numerical solution of differential equations, opening new possibilities for modeling complex physical and engineering systems [2]. Compared to classical numerical methods (such as the finite difference or finite element method), neural network approaches offer several potential advantages: they allow obtaining a solution in the form of a continuous and differentiable function over the entire domain, can work efficiently in high-dimensional spaces, and are well-suited for problems with inverse parameters or incomplete data.

Various neural network architectures are used to solve ordinary differential equations (ODEs). The most common are multilayer networks, such as the classic multilayer perceptron (MLP) [3] or specialized Physics-Informed Neural Networks (PINNs), which explicitly incorporate governing laws into the loss function [4]. However, these architectures often require significant computational resources for training due to the large number of tunable parameters and can suffer from the vanishing gradient problem. An alternative is offered by single-layer neural networks based on the functional link principle (Functional Link Artificial Neural Network, FLANN) [6]. In the FLANN model, the hidden layer is replaced by a functional expansion block that increases the dimensionality of the input data through decomposition into orthogonal polynomial bases (e.g., Chebyshev, Legendre). This significantly reduces the number of network parameters and the number of iterations required for convergence, making FLANN a computationally efficient model with a high training speed.

FLANN networks based on Chebyshev (ChNN) and Legendre (LeNN) polynomials have already been successfully applied to solve various ODEs [7,8,9], demonstrating competitive accuracy. However, the literature lacks a systematic investigation and direct comparison of the effectiveness of these two approaches. Moreover, the potential of neural networks based on other classical orthogonal polynomials, particularly Laguerre polynomials (LaNN), for solving initial and boundary value problems for ODEs remains virtually unexplored. Thus, there is a clear knowledge gap: there is no comparative analysis within the entire family of polynomial FLANNs, nor an assessment of their advantages compared to standard deep architectures.

The aim of this work is to fill this gap through a comprehensive study of the performance of the family of single-layer polynomial neural networks ChNN, LeNN, and LaNN. The Cauchy problem and the boundary value problem for a nonlinear second-order ODE of the Lane-Emden type, for which an exact analytical solution is known, are chosen as test examples. This allows for an objective assessment of the accuracy of each method. The paper conducts a comparative analysis of the polynomial networks among themselves, as well as with the reference multilayer neural network MLP, based on two key criteria: the accuracy of the solution approximation and the algorithm execution time. It is demonstrated that properly configured polynomial networks can outperform MLP in computational efficiency while maintaining high accuracy, making them a promising tool for rapid numerical integration of differential equations.

2. Preliminaries

Orthogonal polynomials are widely used in approximation theory, in particular for constructing quadrature formulas for approximate calculation of integrals or function approximation [10].

Let $f(x)$ and $g(x)$ be real functions belonging to the class $L_w^2(a,b)$.

Definition 1. 

Legendre polynomials of degree n for $x\in [-1,1]$ form an orthogonal system and can be computed by the recurrence formula:

$$P_0(x)=1,P_1(x)=x,P_{n+1}(x)={\displaystyle \frac{2n+1}{n+1}}x{P}_n(x)-{\displaystyle \frac{n}{n+1}}{P}_{n-1}(x).$$

(1)

Definition 2. 

Chebyshev polynomials of degree n for $x\in [-1,1]$ form an orthogonal system and can be computed by the recurrence formula:

$$T_0(x)=1,T_1(x)=x,T_{n+1}(x)=2x{T}_n(x)-T_{n-1}(x).$$

(2)

Definition 3. 

Laguerre polynomials of degree n for $x\in [0,\infty )$ form an orthogonal system and can be computed by the recurrence formula:

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

(3)

Recall that orthogonality of systems (1)-(3) means that the scalar product of any two functions of this system is equal to zero. More details about the properties of orthogonal polynomials can be found in [10].

For simplicity, we introduce the notation $O{P}_m(x)$, which defines the family of orthogonal polynomials (1)-(3).

3. Neural Network Model Based on Orthogonal Polynomials

Consider a single-layer neural network based on the functional-link principle (Figure 1).

Figure 1 shows the structure of a single-layer neural network, which consists of one input node, one output layer and a functional expansion block based on Legendre, Chebyshev and Laguerre orthogonal polynomials. Here, by the notation OPNN we will mean one of the three neural networks LeNN, LaNN or ChNN. The hidden layer is eliminated by transforming the input pattern into a higher-dimensional space using polynomials (1)-(3).

As input data, we consider a vector $x=\left\{x_1,x_2,\dots ,x_n\right\}$ of dimension n. The enhanced pattern is obtained using orthogonal polynomials (1)-(3):

$$\begin{array}{cc}\hfill \left[1,O{P}_1\left(x_1\right),\dots ,O{P}_{m-1}\left(x_1\right);1,O{P}_1\left(x_2\right),\dots ,O{P}_{m-1}\left(x_2\right);\right.& \\ \hfill \left.\dots ;1,O{P}_1\left(x_n\right),\dots ,O{P}_{m-1}\left(x_n\right)\right].\end{array}$$

(4)

According to (4), n-dimensional input x is extended into m-dimensional improved polynomials. Then a weighted sum z of the expanded input data is formed, which is written as:

$$z=\sum _{j=1}^m{w}_{j}O{P}_{j-1}(x).$$

(5)

To update the weights $w_j$ of the neural network, the error backpropagation learning algorithm is used. Thus, the gradient of the error function with respect to the tunable parameters $w_j$ is determined. The nonlinear hyperbolic tangent function $tanh(z)$ is considered as the activation function. For training, the gradient descent algorithm is used, and the weights are updated using the negative gradient at each iteration. The weights are initialized randomly and then updated:

$$w_j^{k+1}=w_j^k-\eta {\left[{\displaystyle \frac{\partial E\left(x,\mathbf{w}\right)}{\partial w_j}}\right]}^k,$$

(6)

where $\eta$ is the learning parameter taking values from 0 to 1, k is the iteration step used to update the weights, as usual in neural networks, and $E\left(x,\mathbf{w}\right)$ is the error function, $\mathbf{w}=\{w_1,\dots ,w_m\}$.

The output layer $y_{NN}(x)$ is determined by the input vector x and the tunable parameters $w_j$ of the selected neural network by the formula:

$$y_{NN}(x)=tanh(z).$$

(7)

The activation function (7) can be chosen differently depending on the problem.

4. Problem Statement

Consider the following problem for a second order ODE for $x\in [0,1]$:

$$y^{″}(x)+p(x)y^{\prime }(x)+q(x)y(x)=f(x),y(0)=y_0,y^{\prime }(0)=y_0^{\prime }.$$

(8)

Problem (8) is a Cauchy problem.

Definition 4. 

The trial solution of the Cauchy problem (8) we will call the function:

$$y_t(x)=y_0+y_0^{\prime }x+x^{2}N(x,\mathbf{w}).$$

(9)

Consider the following boundary value problem for a second order ODE for $x\in [a,b]$:

$$y^{″}(x)+p(x)y^{\prime }(x)+q(x)y(x)=f(x),y(a)=y_a,y(b)=y_b.$$

(10)

Definition 5. 

The trial solution of the boundary value problem (10) we will call the function:

$$y_t(x)=y_a+{\displaystyle \frac{x-a}{b-a}}(y_b-y_a)+(x-a)(b-x)N(x,\mathbf{w}).$$

(11)

The error function for recalculating the weights according to formula (6) taking into account representations (9) or (11) is given as a residual of the following form:

$$E\left(x,\mathbf{w}\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left[y_t^{″}(x_i)+p(x_i)y_t^{\prime }(x_i)+q(x_i)y_t(x_i)-f(x_i)\right]}^2.$$

(12)

The calculation of derivatives (gradients) in formula (12) is performed taking into account formulas (9) and (11). In the algorithms for implementing polynomial networks, gradient computation was performed automatically using the Python library [11].

5. Algorithms for Implementing Polynomial Neural Networks

The polynomial neural network algorithms were implemented in the Python programming language [11] in the PyCharm environment [12] in the OPNN computer program. The OPNN program flowchart is shown in Figure 2. It consists of two stages: the first stage defines the input parameters and selects the appropriate polynomial neural network, while the second stage contains the procedures for training it and inferring the obtained results. The main procedures in the algorithms are presented as pseudocode (Algorithms 1-7).

Note that Algorithms 1-5 are the same for the Cauchy problem and the boundary value problem, but the differences lie in the definition of the trial function (9) and (11), i.e., Algorithms 6 and 7.

Algorithm 1 Neural Network Training for ODE Solution

1: procedure MainTraining($network\_type$, $params$) 2:     Generate training data $x\_data\in [0.1,1]$ 3:     Define exact solution: $exact\_solution(x)=-2ln(1+x^2)$ 4: $$ 5:     if $network\_type=1$ then 6:         $model\leftarrow LegendreNN(n\_polynomials)$ 7:     else if $network\_type=2$ then 8:         $model\leftarrow ChebyshevNN(n\_polynomials)$ 9:     else if $network\_type=3$ then 10:         $model\leftarrow LaguerreNN(n\_polynomials)$ 11:     else 12:         $model\leftarrow MLP(hidden\_size,num\_layers)$ 13:     end if 14: $$ 15:     $start\_time\leftarrow \mathrm{current}\mathrm{time}()$ 16:     $(trained\_model,loss\_history)\leftarrow TrainModel(model,x\_data)$ 17:     $training\_time\leftarrow \mathrm{current}\mathrm{time}()-start\_time$ 18: $$ 19:     Visualize loss history and results 20:     Output comparative value table 21:     return $trained\_model,loss\_history,training\_time$ 22: end procedure

Algorithm 2 Model Training Procedure

1: procedure TrainModel($net$, $x\_data$, $epochs$, $\eta$) 2:     Initialize Adam optimizer with learning rate $\eta$ 3:     Initialize ReduceLROnPlateau scheduler 4:     $loss\_history\leftarrow []$ 5: $$ 6:     for $epoch\leftarrow 1$ to $epochs$ do 7:         $optimizer.zero\_grad()$ 8:         $loss\leftarrow ComputeLoss(net,x\_data)$ 9:         $loss.backward()$ 10:         $optimizer.step()$ 11:         $scheduler.step(loss.item())$ 12:         $loss\_history.append(loss.item())$ 13:         if $epochmod100=0$ then 14:            Output current loss and learning rate 15:         end if 16:     end for 17: $$ 18:     return $net,loss\_history$ 19: end procedure

Algorithm 3 Loss Function Computation

1: procedure ComputeLoss($net$, x) 2:     Set $x.requires\_grad\leftarrow True$ 3:     $(y,dy\_dx,d2y\_dx2)\leftarrow ComputeDerivatives(x,net)$ 4: $$ 5:     Compute equation components: 6:     $term1\leftarrow d2y\_dx2$ 7:     $term3\leftarrow 4\times (2\times exp(y)+exp(0.5\times y))$ 8: $$ 9:     if $x>{10}^{-6}$ then 10:         $term2\leftarrow (2/x)\times dy\_dx$ 11:     else 12:         $term2\leftarrow 0$ 13:     end if 14: $$ 15:     $eq\_residual\leftarrow term1+term2+term3$ 16:     $mask\leftarrow (x>{10}^{-6})$ 17:     $loss\leftarrow \mathrm{mean}(eq\_residual{[mask]}^2)$ 18: $$ 19:     return $loss$ 20: end procedure

Algorithm 4 Trial Solution Derivative Computation

1: procedure ComputeDerivatives(x, $net$) 2:     $x\_tensor\leftarrow x$ (with $requires\_grad=True$) 3:     $y\leftarrow TrialSolution(x\_tensor,net)$ 4: $$ 5:     Compute first derivative: 6:     $dy\_dx\leftarrow {\displaystyle \frac{\partial y}{\partial x}}$ (using automatic differentiation) 7: $$ 8:     Compute second derivative: 9:     $d2y\_dx2\leftarrow {\displaystyle \frac{\partial dy\_dx}{\partial x}}$ (using automatic differentiation) 10: $$ 11:     return $y,dy\_dx,d2y\_dx2$ 12: end procedure

Algorithm 5 Polynomial Neural Network Forward Pass

1: procedure Forward(x, $weights$) 2:     $polys\leftarrow PolynomialBasis(x)$ 3:     $output\leftarrow {\sum }_{i}weights[i]\times polys[:,i]$ 4:     return $output$ 5: end procedure

PolynomialBasis(x) is a procedure for calculating a polynomial basis for a given input x. Depending on the type of neural network, it invokes different procedures according to definitions (1)-(3).

A trial solution to the Cauchy problem according to formula (9) is calculated using Algorithm 6.

Algorithm 6 Trial Solution with Boundary Conditions (The Cauchy Problem)

1: procedure TrialSolution(x, $net$) 2:     $net\_output\leftarrow net(x)$ 3:     $y\leftarrow x^2\times net\_output$ 4:     return y 5: end procedure

A trial solution to the boundary value problem according to formula (11) is calculated using Algorithm 7.

Algorithm 7 Trial Solution with Boundary Conditions (Boundary value problem)

1: procedure TrialSolution(x, $net$) 2:     $net\_output\leftarrow net(x)$ 3:     $y\leftarrow y_a(1-x)+y_{b}x+x(1-x)\times net\_output$ 4:     return y 5: end procedure

6. Research Results

Consider the operation of orthogonal polynomial networks on the example of solving the Cauchy problem for the singular nonlinear Lane-Emden differential equation:

$$y^{″}(x)+{\displaystyle \frac{2}{x}}{y}^{\prime }(x)+4(2e^{y(x)}+e^{y(x)/2})=0,y(0)=0,y^{\prime }(0)=0,$$

(13)

as well as for the boundary value problem:

$$y^{″}(x)+{\displaystyle \frac{2}{x}}{y}^{\prime }(x)+4(2e^{y(x)}+e^{y(x)/2})=0,y(0)=0,y(1)=-2ln2.$$

(14)

Remark 1. 

The Lane-Emden equation is often encountered in astrophysics, in the theory of stellar structure, in the theory of thermal behavior of a spherical gas cloud and the theory of thermionic currents [8,13].

Remark 2. 

It is known from works [8,13] that the Cauchy problem (13) has an exact solution of the form:

$$y(x)=-2ln(1+x^2).$$

(15)

Based on Definition 4 by formula (9), the trial solution for the Cauchy problem (13) is written as:

$$y_t(x)=x^{2}N(x,\mathbf{w}),$$

and for the boundary value problem (14), also valid is the exact solution (15) and taking into account representation (11) in the form:

$$y_t(x)=-2ln2·x+x(1-x)N(x,\mathbf{w}).$$

The parameters of the orthogonal polynomial neural networks were set as follows: order of polynomials $n=5$, number of training points 11 from the segment $[0.1,1]$, number of epochs 10000, learning rate $\eta =0.01$. The calculation results for solving the Cauchy problem (13) are given in Table 1.

From Table 1 it can be seen that the polynomial neural networks give approximately the same result. Figure 3 shows the solution graph obtained by the LaNN neural network, and Figure 4 shows the loss graph during training.

The graph in Figure 3 demonstrates the efficiency of the Laguerre NN neural network in solving the Cauchy problem (13). The loss function graph on a logarithmic scale is shown in Figure 4.

From Figure 4 it can be seen that the loss decreases with an increase in the number of epochs, which indicates successful model training. At the beginning of training (first 200 epochs), a rapid decrease in loss is observed, then the rate of decrease slows down, which is typical for the neural network training process.

The loss graph allows adjusting the hyperparameters of the neural network: learning rate, number of epochs, etc.

The solution of the boundary value problem (14) is given in Table 2.

Here it can be noted that the LaNN neural network solved the problem more accurately than the ChNN network, but was less accurate than the LeNN network.

Table 3 shows the work of the multilayer neural network MLP with two hidden layers of 20 and 40 neurons respectively for solving the Cauchy problem (13).

It can be seen from Table 3 that to achieve acceptable accuracy in MLP, it is necessary to increase the number of neurons in each layer. This leads, as we will show below, to an increase in the algorithm execution time.

Table 4 shows the solution of the boundary value problem (14) using the MLP network.

From Table 4 it can be seen that increasing the number of neurons in each hidden layer leads to a deterioration in the accuracy of solving the boundary value problem (14). Let us show how the execution time of algorithms for solving problems (13) and (14) depends on the choice of neural network architecture. Neural network parameters: learning rate $\eta =0.01$, number of epochs – 1000, order of polynomials $n=5$.

Table 5 shows the times spent on executing the algorithm for solving the Cauchy problem (13).

From Table 5 it can be seen that $\mathrm{ChNN}<\mathrm{MLP}(20,2)<\mathrm{LaNN}<\mathrm{LeNN}$. Here it should be noted that the neural network based on Laguerre orthogonal polynomials coped with the task faster than the neural network based on Legendre orthogonal polynomials. Moreover, the LaNN and ChNN neural networks were faster in solving the Cauchy problem (13) than the multilayer MLP network.

The results of solving the boundary value problem (14) are given in Table 6.

From Table 6 it follows that $\mathrm{ChNN}<\mathrm{MLP}(20,2)<\mathrm{LaNN}<\mathrm{LeNN}$, and therefore the trend in the execution times of algorithms for solving the boundary value problem (14) remains as in Table 5. Here, a multilayer neural network MLP with three hidden layers of 20 neurons in each layer was used. From Table 6 it follows that single-layer orthogonal polynomial neural networks can work faster than multilayer neural networks of the MLP type, for example, the LeNN neural network.

It should be noted that increasing the order of the polynomial n can also lead to a deterioration in the result. Let us take the parameters for the LeNN neural network from the previous example (Table 3) and find solutions to problems (13) and (14) for polynomial orders $n=5$ and $n=7$. The research results are given in Table 7 and Table 8.

From Table 7 we see that with an increase in the order of the polynomial, the calculation error increased, although the algorithm execution time decreased from 10.01 s for $n=5$ to 8.92 s for $n=7$.

From Table 8 we see that increasing the order of the Legendre polynomial leads to a deterioration in accuracy, and also to a decrease in the algorithm execution time for Table 8 at $n=5$ – 6.54 s, and at $n=7$ – 8.66 s. Therefore, the order of the polynomial must be chosen optimally.

7. Overfitting Issues

It should be noted that the issue of overfitting remains relevant when using neural networks to solve differential equations. Overfitting here in the context of the problem of solving a differential equation is understood differently. Here, the training sample is the differential equation itself, which essentially defines an infinite set of constraints that the solution function must satisfy. The neural network here tries to approximate this function that satisfies these constraints everywhere in a given domain. It should also be said that there is no division of data into test and training samples here. We are looking for a single solution for the entire domain. The concept of "new data" from the same domain is erased — we want the network to work well on the entire domain at once. Thus, pure "memorization of noise in the training data" is not applicable here, since there are no discrete data initially. The concept of overfitting is understood mainly as overfitting at collocation points. We cannot check the fulfillment of the equation on the entire infinite domain. Therefore, we choose a finite set of points (collocation points) at which we calculate the residual of the equation and minimize it. The neural network between collocation points can behave incorrectly, for example, oscillate. The residual at the training collocation points will be very low (sometimes almost zero), but if we check the solution on another, denser grid of points, the residual will be huge. This is a sign of overfitting. Similarly, the concept of overfitting can arise at boundary (initial) points of the domain.

Here, overfitting can be combated by various methods similar to solving the issue of classical overfitting. For example, known methods are L1/L2 regularization, adaptive weights for the loss function, etc. However, in our opinion, the most effective and common method is dynamic selection of collocation points. Instead of using a fixed set of points, at each iteration (or after several epochs) collocations are selected anew (often randomly). Here it is necessary to remember about the non-optimal choice of collocation points. If the points at which the equation residual is calculated are chosen unsuccessfully (for example, uniformly over the entire domain), the network may learn well in some areas and poorly in others (for example, in areas of high gradients) [14]. This approach does not allow the network to "memorize" their location and forces it to learn to satisfy the equation in the entire domain. We also note article [4], in which the authors study the issues of overfitting of the multilayer PINN caused by gradient imbalance.

8. Conclusions

The work for the first time investigated the approximate solution of the Cauchy problem and the boundary value problem for the Lane-Emden differential equation using three single-layer neural networks based on Legendre, Chebyshev and Laguerre orthogonal polynomials. It is shown that the new neural network based on Laguerre polynomials gives acceptable results, and in some cases can be more accurate than the ChNN neural network and faster than the LeNN neural network. All three polynomial neural networks were also compared with the multilayer MLP neural network. Research has shown that single-layer polynomial neural networks can work better than the multilayer MLP neural network. Key parameters for polynomial neural networks, in addition to choosing the number of training epochs and learning rate, is the order of the polynomial. However, this order must be chosen optimally; its further increase or decrease can lead to a deterioration in results. Note that in multilayer neural networks it is necessary to adjust the number of hidden layers and the number of neurons in each of them. The latter means that single-layer polynomial neural networks have a simpler architecture than multilayer neural networks and, accordingly, they are easier to implement on a computer. The choice of neural network architecture depends on the specific task for an ordinary differential equation.

In the future, it is planned to apply polynomial neural networks to solve fractional differential equations and their systems [15,16,17,18,19], as well as to partial differential equations by analogy with work [20].