Physical Information Neural Networks for Solving High-index
  Differential-algebraic Equation Systems Based on Radau Methods
As is well known, differential algebraic equations (DAEs), which are able to
describe dynamic changes and underlying constraints, have been widely applied
in engineering fields such as fluid dynamics, multi-body dynamics, mechanical
systems and control theory. In practical physical modeling within these
domains, the systems often generate high-index DAEs. Classical implicit
numerical methods typically result in varying order reduction of numerical
accuracy when solving high-index systems.~Recently, the physics-informed neural
network (PINN) has gained attention for solving DAE systems. However, it faces
challenges like the inability to directly solve high-index systems, lower
predictive accuracy, and weaker generalization capabilities. In this paper, we
propose a PINN computational framework, combined Radau IIA numerical method
with a neural network structure via the attention mechanisms, to directly solve
high-index DAEs. Furthermore, we employ a domain decomposition strategy to
enhance solution accuracy. We conduct numerical experiments with two classical
high-index systems as illustrative examples, investigating how different orders
of the Radau IIA method affect the accuracy of neural network solutions. The
experimental results demonstrate that the PINN based on a 5th-order Radau IIA
method achieves the highest level of system accuracy. Specifically, the
absolute errors for all differential variables remains as low as $10^{-6}$, and
the absolute errors for algebraic variables is maintained at $10^{-5}$,
surpassing the results found in existing literature. Therefore, our method
exhibits excellent computational accuracy and strong generalization
capabilities, providing a feasible approach for the high-precision solution of
larger-scale DAEs with higher indices or challenging high-dimensional partial
differential algebraic equation systems.
Jiasheng Chen, Juan Tang, Ming Yan, Shuai Lai, Kun Liang, Jianguang Lu, Wenqiang Yang
Physical Information Neural Networks for Solving High-index Differential-algebraic Equation Systems Based on Radau Methods

###### Abstract

As is well known, differential algebraic equations (DAEs), which are able to describe dynamic changes and underlying constraints, have been widely applied in engineering fields such as fluid dynamics, multi-body dynamics, mechanical systems and control theory. In practical physical modeling within these domains, the systems often generate high-index DAEs. Classical implicit numerical methods typically result in varying order reduction of numerical accuracy when solving high-index systems. Recently, the physics-informed neural network (PINN) has gained attention for solving DAE systems. However, it faces challenges like the inability to directly solve high-index systems, lower predictive accuracy, and weaker generalization capabilities. In this paper, we propose a PINN computational framework, combined Radau IIA numerical method with a neural network structure via the attention mechanisms, to directly solve high-index DAEs. Furthermore, we employ a domain decomposition strategy to enhance solution accuracy. We conduct numerical experiments with two classical high-index systems as illustrative examples, investigating how different orders of the Radau IIA method affect the accuracy of neural network solutions. The experimental results demonstrate that the PINN based on a 5th-order Radau IIA method achieves the highest level of system accuracy. Specifically, the absolute errors for all differential variables remains as low as \(10^{-6}\), and the absolute errors for algebraic variables is maintained at \(10^{-5}\), surpassing the results found in existing literature. Therefore, our method exhibits excellent computational accuracy and strong generalization capabilities, providing a feasible approach for the high-precision solution of larger-scale DAEs with higher indices or challenging high-dimensional partial differential algebraic equation systems.

keywords: Differential algebraic equation, Radau IIA method, Physics-informed neural network, Domain decomposition +
Footnote †: journal:

## 1 Introduction

The concept of differential algebraic equations (DAEs) was formally proposed by Gear in the study of network analysis and continuous system simulation problems.[1] Petzold made it explicit through his study of numerical methods that DAEs are not ordinary differential equations (ODEs).[2] DAE systems are composed of coupled ODE systems and algebraic equation systems with physical significance. These systems encompass both differential and algebraic variables, and their system form is more generalized compared to traditional ODE systems. DAEs have gained significant attention since their inception, as they can accurately describe systems that some ODEs cannot represent. They have found extensive applications in various fields, including fluid dynamics, multi-body dynamics, electronic circuits, mechanical systems, control theory, and chemical engineering.

In different developmental periods and research fields, DAEs are also known as singular systems, general systems, descriptor systems, or constrained systems, among other names. They often exhibit various structural forms, such as linear DAEs, nonlinear DAEs, semi-explicit DAEs, implicit DAEs, and Hessenberg-type DAEs. Fortunately, in practical physical modeling, most of the system models obtained are either low-index DAEs or high-index (\(>\)=2) Hessenberg-type DAEs.[3] The index of DAEs measures the 'distance' between DAEs and ODEs. Generally, a higher index implies greater difficulty in transforming DAEs into ODEs or in directly solving DAEs using ODE numerical methods. Traditional numerical methods for solving DAE systems include implicit Runge-Kutta methods[4], BDF methods[5], pseudospectral methods[6], adomian decomposition method[7], exponential integrators[8], generalized-\(\alpha\) methods[9], and Lie group methods[10; 11; 12]. It's worth noting that these direct numerical methods can solve DAEs with an index of 1. However, for high-index DAE systems, these methods are only applicable to a certain class of DAEs and may result in varying order reduction of numerical accuracy.

With the rapid advancement of neural network technology and hardware resources, neural networks are demonstrating increasingly powerful capabilities. Compared to traditional nu
merical computing methods, neural networks offer several advantages, including strong generalization, fault tolerance, and the ability for parallel computation. In 1998, Lagaris et al. [13] approximated solutions to ODEs or PDEs problems by constructing parameterized trial functions. These trial functions consist of two parts: one part satisfies initial conditions or boundary conditions which does not contain trainable parameters, while the other part is a simple feed forward neural network with trainable parameters. In 2019, Raissi et al. [14] introduced an important technique known as Physics-informed Neural Network (PINN) for the numerical approximation of partial algebraic equations (PDEs) problems. The PINN loss function includes not only initial or boundary conditions that reflect physical properties but also a residual term at selected points in the time-space domain where the PDEs hold. It's worth noting that PINN is a data-driven approach that doesn't require prior knowledge of the analytical form of the solution; instead, it learns the solution from data. Various variants of PINN have been proposed based on different collocation methods, such as variational hp-VPINN [15] and conservative PINN (CPINN) [16]. Additionally, PINN has been widely applied to solve problems in various fields, including fluid dynamics [17; 18], seismic wave prediction [19], and optical problems [20].

In recent years, many researchers have attempted to construct neural network models from different perspectives to solve various types of DAEs systems influenced by these methods. For Hessenberg-DAEs with control variables and an index of 3, Kozlov and Tiumentsev [21] achieved the implementation of BDFs method using a semi-empirical neural network model. Zhao Yang et al. [22] constructed a single-layer feed-forward neural network (FFNN) to solve Hessenberg-type DAEs systems. They augmented the loss function in their special Euler-Lagrange equation system with penalty terms for algebraic equations to avoid drifting in the results. Experimental results in their paper showed that the FFNN method with Sigmoid activation function provided approximate analytical solutions close to the numerical solutions of corresponding Runge-Kutta methods, but they didn't provide further details about the method's accuracy. For linear DAEs systems, Hongliang Liu et al. [23] selected Jacobi polynomials as activation functions and constructed a single-hidden-layer feed-forward neural network (JNN). They determined the network parameters using the classical ELM algorithm. Through experimental comparisons with other approximation methods such as Pade approximation, ADM method, and Adams methods, they illustrated the feasibility and superiority of the JNN method. It's worth noting that the examples in the paper involve DAEs with an index of 1 or linear DAEs that have been reduced to index 1. For DAEs systems with an index of 1, Moya et al. [24] proposed a neural network architecture called DAE-PINN based on the PINN method for solving DAEs systems. This neural network model is a discrete-time model based on the implicit Runge-Kutta method, which can directly address most index-1 differential-algebraic equation problems. However, it cannot solve high-index DAEs problems and suffers from low accuracy issues. To address the high-accuracy computation challenges in high-index DAEs systems, we have combined the Radau IIA numerical method with a neural network structure based on attention mechanisms. We have proposed a PINN computational framework based on the Radau method. Furthermore, we have improved the efficiency and accuracy of the solution by applying a strategy of domain decomposition.

In section 2, we briefly introduce the fundamental concepts of DAEs systems, the Radau IIA numerical method, and the neural network structure based on attention mechanisms. Building upon this foundation, we provide a detailed construction of the PINN computing framework based on the Radau IIA method. Additionally, we employ a time domain decomposition strategy for neural network. Section 3 use the neural network designed in this paper to solve two high-index DAEs systems, and we analyze the solving accuracy of this neural network. Finally, we discuss and summarize the advantages, challenges, and potential avenues for improvement in the Radau-PINN architecture.

## 2 Scientific Machine Learning Methods

This section first sequentially introduces the basic concepts of DAEs and the classical Radau IIA numerical method. Then, we introduce a neural network structure based on attention mechanisms. Building upon this, we construct a PINN based on the Radau IIA method. Finally, we enhance the efficiency and accuracy of neural network solutions for DAEs systems by utilizing the concept of domain decomposition.

### Radau IIA Method for DAE Systems

This article first provides a brief introduction to DAEs with an index of 2, with the specific form as follows:

\[\begin{cases}y^{\prime}(t)=&f(t,y(t),z(t)),\\ &0=&g(t,y(t)),\end{cases} \tag{1}\]

where \(y(t)\in\mathbb{R}^{n}\) is the differential function variable, \(z(t)\in\mathbb{R}^{m}\) is the algebraic function variable, \(t\in[t_{0},T]\), \(t_{0}\) is the initial time point, and \(y_{0}=y(t_{0})\) is the initial value. Both \(f(t,y,z)\in\mathbb{R}^{n}\) and \(g(t,y)\in\mathbb{R}^{m}\) are sufficiently smooth, and the Jacobian matrix \(g_{y}f_{z}\) is non-singular.

The Radau IIA method is a class of implicit Runge-Kutta methods, typically defined in the following general form:

\[\xi_{i}=y_{n}+h\sum_{j=1}^{v}a_{i,j}f\left(\xi_{j},\xi_{j}\right), \tag{2}\] \[g\left(\xi_{i},\zeta_{i}\right)=0,\] (3) \[y_{n+1}=y_{n}+h\sum_{j=1}^{v}b_{j}f\left(\xi_{j},\xi_{j}\right),\] (4) \[g\left(y_{n+1},z_{n+1}\right)=0, \tag{5}\]

where \(\xi_{i}=y\left(t_{n}+c_{i}h\right)\), \(\zeta_{i}=z\left(t_{n}+c_{i}h\right)\), \(h\) is the step size, \(n\) is the current step number, \(\left\{a_{ij},b_{j},c_{i}\right\}\) are parameters, and \(c_{i}=\sum\limits_{j=1}^{v}a_{ij}\), \(i,j=1,\cdots,v\).

In table 1, different sets of parameters lead to different implicit Runge-Kutta methods, such as commonly used Gauss
method, Radau method, and Lobatto method. These parameters are determined using Gauss polynomials, Radau polynomials, and Lobatto polynomials, respectively. Among them, the Radau IIA method is a high-precision numerical method with excellent numerical stability. Therefore, in this paper, the Radau IIA method is chosen, and the parameters need to satisfy the following conditions:

\[B(2v-1):\sum_{i=1}^{v}b_{i}c_{i}^{k-1}=\frac{1}{k},\ k=1,\cdots,2v-1 \tag{6}\] \[C(v):\sum_{j=1}^{v}a_{ij}c_{j}^{k-1}=\frac{c_{i}^{k}}{k},\ k=1, \cdots,v\] (7) \[D(v-1):\sum_{i=1}^{v}b_{i}c_{i}^{k-1}a_{ij}=\frac{b_{j}}{k}(1-c_{ j}^{k}),\ k=1,\cdots,v-1 \tag{8}\]

and \(c_{v}=1\), \(b_{j}=a_{vj}\), \(i,j=1,2,\ldots,v\).

### Neural Network Structure Based on Attention Mechanism

Building upon the DAE-PINN structure, we employ adaptive activation functions (9) to train a neural network structure based on an attention mechanism. The specifics are as follows:

The improved neural network model architecture based on attention mechanisms is primarily constructed using two Transformer networks, denoted as U and R, to build two stacked layer networks, as illustrated in Figure 1. Both neural networks map the input variable \(X\) (differential function variable \(y\)) to a high-dimensional feature space. Subsequently, each hidden layer forms new residual connections using element-wise multiplication operations, as expressed below:

\[U = \phi(XW^{1}+b^{1}), \tag{9}\] \[R = \phi(XW^{2}+b^{2}),\] (10) \[H^{(1)} = \phi(\eta\cdot l\cdot XW^{o,1}+b^{o,1}),\] (11) \[M^{(k)} = \phi(H^{k}W^{o,k}+b^{o,k}),\] (12) \[H^{(k+1)} = (1-M^{(k)})\odot U+M^{(k)}\odot R,\] (13) \[P_{\theta}(X) = H^{d+1}W+b, \tag{14}\]

where \(X\) represents the input vector of the neural network, \(W^{o,k}\) is the collection of weights for the \(o\)-th neuron in the \(k\)-th layer, \(b^{o,k}\) denotes the set of biases for the \(o\)-th neuron in the \(k\)-th layer, \(\phi\) is the activation function, \(\odot\) represents element-wise multiplication, \(d\) indicates the number of hidden layers (the depth of the neural network), \(P_{\theta}(X)\) is the final output vector of the neural network, \(\eta\) is a predetermined hyper-parameter that ensures the slope is greater than 1, and \(l\) is a parameter that can modify the slope of the activation function.

### PINN Based on Radau IIA Method

In this section, we use the discrete-time model of PINN as the foundation, incorporating a neural network structure based on attention mechanisms. We have constructed a PINN architecture based on the Radau IIA method, as illustrated in Figure 2.

Firstly, construct a neural network with multiple inputs and multiple outputs, where the inputs consist of the collection of differential variables \(y_{n}\), and outputs

\[\xi_{1}^{\theta},\xi_{2}^{\theta},.....,\xi_{v}^{\theta},{y_{n+1} }^{\theta}; \tag{15}\] \[\zeta_{1}^{\theta},\zeta_{2}^{\theta},.....,\zeta_{v}^{\theta},z _{n+1}}^{\theta}. \tag{16}\]

The first \(v\) values of \(\xi_{i}^{\theta}\) represent intermediate differential variables, and the first \(v\) values of \(\zeta_{i}^{\theta}\) represent intermediate algebraic variables, where \(i=1,2,.....,v\).

Secondly, based on the structure of the DAEs system with an index of 2 and the characteristics of the Radau IIA method, we further design neural network structures based on the attention mechanism for both the differential variable part and the algebraic variable part of the system. There are two specific design approaches: one assigns a single neural network to all differential variables and two neural networks to the algebraic variables, and the other assigns a separate neural network to each differential variable while keeping the algebraic variable part unchanged. In the case of the algebraic variable part, one of the neural networks is used to predict the first v values, while the other neural network is used to predict the v+1-th value. Theoretically, the second approach (as shown in Figure 2) constructs a neural network for each individual differential or algebraic variable, thereby improving the overall model's accuracy and generalization. Through further testing, the second approach's training results are more precise than those of the first approach, consistent with the expected results. As a result, all subsequent experiments in this paper are implemented based on the second approach.

Thirdly, based on the designed network structure, this paper constructs the loss function as follows:

\[\mathcal{L}\left(\theta;\mathcal{T}\right)=W_{f}\mathcal{L}_{f}\left(\theta; \mathcal{T}\right)+W_{g}\mathcal{L}_{g}\left(\theta;\mathcal{T}\right)+W_{s} \mathcal{L}_{s}\left(\theta;\mathcal{T}\right), \tag{17}\]
* \(\mathcal{L}_{f}\left(\theta;\mathcal{T}\right)\) is the loss related to the differential network and is expressed as follows: \[\frac{1}{N_{T}\left(v+1\right)}\sum_{k=1}^{N_{T}}\sum_{i=1}^{v+1} \left\|v_{n,k}-y_{n,k}^{i}\left(\theta\right)\right\|_{2}^{2},\] (18) \[y_{n,k}^{i}(\theta)=\xi_{i,k}^{\theta}-h\sum_{j=1}^{v}a_{i,j}f( \xi_{j,k}^{\theta},\xi_{j,k}^{\theta}),\] (19) \[i=1,\cdots,v,\] \[y_{n,k}^{v+1}(\theta)=y_{n+1,k}^{\theta}-h\sum_{i=1}^{v}b_{i}f( \xi_{i,k}^{\theta},\xi_{i,k}^{\theta});\] (20)
* \(\mathcal{L}_{g}\left(\theta;\mathcal{T}\right)\) is the loss associated with the algebraic network

\begin{table}
\begin{tabular}{c|c c c c c} \(c_{1}\) & \(a_{11}\) & \(a_{12}\) & \(a_{13}\) & \(\cdots\) & \(a_{1v}\) \\ \(c_{2}\) & \(a_{21}\) & \(a_{22}\) & \(a_{23}\) & \(\cdots\) & \(a_{2v}\) \\ \(c_{3}\) & \(a_{31}\) & \(a_{32}\) & \(a_{33}\) & \(\cdots\) & \(a_{3v}\) \\ \(\vdots\) & \(\vdots\) & \(\vdots\) & \(\vdots\) & \(\ddots\) & \(\vdots\) \\ \(c_{v}\) & \(a_{v1}\) & \(a_{i2}\) & \(a_{v3}\) & \(\cdots\) & \(a_{vv}\) \\ \hline  & \(b_{1}\) & \(b_{2}\) & \(b_{3}\) & \(\cdots\) & \(b_{v}\) \\ \end{tabular}
\end{table}
Table 1: The parameter table of the \(v\)-stage implicit Runge-Kutta methods.
and can be expressed as follows:

\[\frac{1}{N_{\mathcal{T}}\left(\nu+1\right)}\sum_{k=1}^{N_{\mathcal{T}}} \left(\sum_{i=1}^{v}\left\|g\left(\xi_{i,k}^{\theta},\xi_{i,k}^{\theta}\right) \right\|_{2}^{2}+\left\|g\left(y_{n+1,k}^{\theta},\xi_{n+1,k}^{\theta}\right) \right\|_{2}^{2}\right); \tag{21}\]
* \(\mathcal{L}_{s}\left(\theta;\mathcal{T}\right)\) is the loss related to the last value of the controlled algebraic variable and can be expressed as follows: \[\frac{1}{N_{\mathcal{T}}}\sum_{k=1}^{N_{\mathcal{T}}}\left(\left\|\xi_{n+1,k} ^{\theta}-\zeta_{v,k}^{\theta}\right\|_{2}^{2}\right),\] (22) where \(W_{f}\) represents the loss weight for the differential neural network, \(W_{g}\) is the weight for the algebraic neural network, and \(W_{s}\) signifies the weight for the control of algebraic variable prediction neural network. The parameters \(a_{i,j}\) and \(b_{i}\) are specific to the Radau IIA method. \(\mathcal{T}\) is the total number of samples, \(N_{\mathcal{T}}\) is the number of training samples in the current batch, and \(\theta\) denotes the neural network parameters. Here, \(f\) represents the differential network, and \(g\) represents the algebraic network. \(y_{n,k}\) corresponds to the sample data of the model, \(\xi_{i,k}^{\theta}\) stands for the values of intermediate differential variables, and \(\zeta_{i,k}^{\theta}\) signifies the values of intermediate algebraic variables. Furthermore, \(y_{n,k}^{i}\left(\theta\right)\) represents the output values of the differential neural network. The notation \(\left\|\cdot\right\|_{2}^{2}\) refers to the square of the L2 norm, \(z_{n+1,k}^{\theta}\) represents the final output of the algebraic neural network, and \(\zeta_{n,k}^{\theta}\) denotes the penultimate output of the algebraic neural network.

Finally, we use gradient descent to solve for the weights, biases, and other parameters of the PINN,

\[\theta^{*}=\arg\min_{\theta}\mathcal{L}\left(\theta;\mathcal{T}\right). \tag{23}\]

### Time Domain Decomposition of Neural Networks

In this section, based on an analysis of the existing limitations of the PINN architecture, we adopt a time-domain decomposition strategy using neural networks.

One limitation of the PINN model is that it exhibits relatively low accuracy in predicting solutions. This is because the inherent inaccuracies involved in solving high-dimensional non-convex optimization problems can lead to local minima, making it challenging to achieve absolute errors below \(10^{-5}\). Another evident limitation is the high training cost.[16] Similarly,

Figure 1: Improved attention neural network structure.

Figure 2: The schematic diagram of PINN based on Radau IIA method.
our proposed PINN model based on the Radau IIA method may encounter similar issues. Furthermore, the iterative format of the Radau IIA method does not fully exploit its high-precision advantages during training.

To address these issues, we propose a time-domain decomposition strategy for neural networks, as illustrated in Figure 3. With this approach, we partition the original problem into segments, which not only enhances solution accuracy but also leverages the advantages of iterative training. In other words, the predicted values from the previous time segment can serve as input values for the subsequent segment. This means that knowing the data values at the initial point \(t_{0}\) for the first segment is sufficient to iteratively compute the solutions over the entire time domain. This approach significantly reduces the amount of required data. Specifically, only the data at the initial point \(t_{0}\) for a set of differential variables, denoted as \(y_{0}\), is needed. Using the time-domain decomposition structure, we can iteratively determine the desired values within the range \([t_{0},T]\). This involves information related to \((T-t_{0})/h\cdot v\) data points, which reduces the need for extensive training data. On the other hand, if we can obtain the initial values for each network at every time segment, parallel training of each neural network becomes possible, significantly reducing the model training time.

## 3 Numerical Experiments

In this section, we apply PINN based on the Radau IIA method to solve two high-index DAEs systems separately and further investigate the influence of the order of the Radau IIA method on the solution results. The experiments were conducted on a Windows 10 operating system with an Intel(R) Core(TM) i7-10875H CPU @ 2.30GHz processor. We used Python 3.9 software and coded the neural network architecture using PyTorch 1.12.1, the GPU version. Additionally, this paper involves two formulas to measure the accuracy of the experiments. One is the commonly used Absolute Error (AE) formula, defined as \(AE=|y_{true}-y_{pred}|\), which reflects the magnitude of the deviation between the neural network's predicted solution and the true solution. The other metric is the Mean Absolute Error (MAE) formula, defined as \(MAE=\frac{1}{n}\sum_{i=1}^{n}\|y_{true}^{i}-y_{pred}^{i}|\), used to assess the differences in accuracy among different orders of the Radau IIA method.

### Hessenberg-type DAEs System

In this section, we explore classical Hessenberg-type DAE systems with an index of 2 that possess exact analytical solutions [12], as follows:

\[\begin{cases}y_{1}^{\prime}(t)=&(y_{3}(t)y_{4}(t)+y_{1}(t)y_{2}(t))y_{5}(t),\\ y_{2}^{\prime}(t)=&-y_{3}(t)y_{4}(t)y_{2}(t)^{2}y_{5}(t),\\ y_{3}^{\prime}(t)=&2y_{3}(t)y_{4}(t)y_{1}(t)y_{2}(t),\\ y_{4}^{\prime}(t)=&-y_{3}(t)y_{4}(t)y_{2}(t)^{2},\\ 0=&y_{1}(t)y_{4}(t)-y_{2}(t)y_{3}(t),\end{cases} \tag{24}\]

where \(t\in[0,1]\), and the initial values \(y_{0}=(1,1,1,1,1)\). The functions \(y_{1}(t),y_{2}(t),y_{3}(t),y_{4}(t)\) represent differential variables, while \(y_{5}(t)\) is an algebraic variable. The system's exact solution expressions are \(y_{1}(t)=e^{2t}\), \(y_{2}(t)=e^{-t}\), \(y_{3}(t)=e^{2t}\), \(y_{4}(t)=e^{-t}\), and \(y_{5}(t)=e^{t}\).

Firstly, we consider the impact of different orders of Radau IIA methods, including 3rd, 5th, 9th, and 13th orders (corresponding to \(v=2,3,5,7\)), on the precision of neural network solutions. Secondly, we explore the influence of activation functions on PINN. Common activation functions for hidden layers include Sigmoid, TanH, Sin, and ReLu, among others. When solving smoothly continuous systems, ReLu is generally not chosen; instead, Sigmoid, TanH, or Sin activation functions are preferred. In the experiments, _Sigmoid_ resulted in better approximate solutions. Within this neural network framework, the initial values of the differential variables, namely \(y_{1}(t),y_{2}(t)\), \(y_{3}(t)\), and \(y_{4}(t)\) for each time segment, are used as a dataset for training. The step size \(h\) is 0.05, which means that each time

Figure 3: The time domain decomposition of neural networks.
interval has a length of 0.05. Each network model in every time segment comprises 5 hidden layers, with each hidden layer containing 100 neurons. Sigmoid is used as the activation function, and the Adam optimizer is applied for 100,000 iterations. The experimental results within the time interval of 0 to 1 are presented in Figures 4.

From Figure 4, it is evident that the accuracy of the mean absolute errors for the 3rd and 13th-order Radau IIA methods corresponds to the blue Y-axis, while the accuracy of the average absolute errors for the 5th and 9th-order Radau methods corresponds to the red Y-axis. For all the differential function variables, the 3rd and 13th-order Radau IIA methods exhibit significantly higher average absolute errors compared to the 5th and 9th-order methods. For the algebraic variable \(y_{5}\), the 13th-order Radau IIA method has notably higher average absolute errors than the 3rd, 5th, and 9th-order methods.

Additionally, we further observe that for all differential function variables from red Y-axis, the 9th-order Radau IIA method's overall trend in average absolute errors is significantly higher than the 5th-order method. For the algebraic variable \(y_{5}\), the 9th-order Radau IIA method exhibits notably higher average absolute errors than the 5th-order method. In other words, the 5th-order Radau IIA-based PINN achieves the highest precision in terms of average absolute errors.

The absolute error results obtained using the 5th-order method are shown in Figure 5. The accuracy of the absolute errors for \(y_{1}(t)\), \(y_{3}(t)\), and \(y_{5}(t)\) corresponds to the blue Y-axis, while the accuracy of the absolute errors for \(y_{2}(t)\) and \(y_{4}(t)\) corresponds to the red Y-axis. From the figure, it is evident that the neural network's predicted values for all four differential variables have their lowest precision of absolute errors maintained at the order of \(10^{-6}\), while the lowest precision of absolute errors for the algebraic variable is kept at \(10^{-6}\). The experimental results suggest that the neural network's predicted solutions have reached a high level of accuracy.

For the neural network structure designed in this paper, the predicted values of the differential variables can be used as the initial values for the next time step's network input dataset. The

Figure 4: The mean absolute errors for solving Hessenberg-DAEs systems using PINN based on Radau IIA of order \(v=3,5,9,13\). The blue curves represent the mean absolute errors \(v=3,13\) on the left Y-axis, while the red curves correspond to \(v=5,9\) on the right Y-axis.
precision of the differential variables can affect the results of the next time step's network. In this context, the precision of the differential variables \(y_{1}\) and \(y_{3}\) is already at the order of \(10^{-6}\), and the precision of the differential variables \(y_{2}\) and \(y_{4}\) is at the order of \(10^{-7}\), which will not significantly affect the precision of the next time step.

### DAE System of the Pendulum Model

In this section, we study the classical pendulum DAEs system with an index of 2, as follows:

\[\left\{\begin{array}{l}y_{1}^{\prime}(t)= y_{3}(t),\\ y_{2}^{\prime}(t)= y_{4}(t),\\ y_{3}^{\prime}(t)=- y_{1}(t)y_{5}(t),\\ my_{4}^{\prime}(t)=- y_{2}(t)y_{5}(t)-\lambda,\\ 0= y_{1}(t)y_{3}(t)+y_{2}(t)y_{4}(t),\end{array}\right. \tag{25}\]

where \(t\in[0,1]\), and the parameters \(m\) and \(\lambda\) are variable parameters, both set to 1 in the experiments of this section. The initial values are \(y_{0}=(1,0,0,1,1)\). In this context, \(y_{1}(t)\), \(y_{2}(t)\), \(y_{3}(t)\), and \(y_{4}(t)\) are differential function variables, while \(y_{5}(t)\) is an algebraic function variable. This DAEs system does not have an exact analytical expression. In this paper, we directly solve the reduced inner ODEs of this system using high-precision ODE solvers from the Python scientific computing library _Scipy_ and compare the obtained approximate solution with the predicted values from the neural network. Similarly, we consider the impact of different orders (3, 5, 9, 13, corresponding to \(v=2,3,5,7\)) in the Radau IIA methods on the accuracy of the neural network's solutions. Secondly, we explore the effect of activation functions on PINN. In this experiment, the _Sin_ activation function provides a better approximation. To maintain consistency in the numerical experiments, other network structural information is consistent with the experiments in the previous section. The results obtained are shown in Figures 6.

From Figure 6, we can observe that the accuracy of the mean absolute errors for the 3rd and 13th-order Radau IIA methods corresponds to the blue Y-axis, while the accuracy of the average absolute errors for the 5th and 9th-order Radau methods corresponds to the red Y-axis. For the differential function variables \(y_{1}(t)\), \(y_{3}(t)\), and \(y_{4}(t)\), the 3rd-order Radau IIA method has significantly higher average absolute errors than the 5th, 9th, and 13th-order methods. For the differential function variable \(y_{2}(t)\), the 13th-order Radau IIA method exhibits significantly higher average absolute errors than the 3rd, 5th, and 9th-order methods. For the algebraic variable \(y_{5}(t)\), the 3rd and 13th-order Radau IIA methods have significantly higher average absolute errors compared to the 5th and 9th-order methods.

Additionally, we further observe that for all differential function variables from red Y-axis, the 9th-order Radau IIA method's average absolute error overall trends similarly to the 5th-order method. For the algebraic variable \(y_{5}(t)\), the 9th-order Radau IIA method exhibits significantly higher average absolute errors in the later time regions compared to the 5th-order method. In other words, a PINN based on the 5th-order Radau

Figure 5: The absolute errors of Hessenberg-DAEs system solved by PINN based on 5th-order Radau IIA. The blue curves represent the absolute errors of \(y_{1}\), \(y_{3}\) and \(y_{5}\) on the left Y-axis, while the red curves correspond to \(y_{2}\) and \(y_{4}\) on the right Y-axis.
IIA method achieves the highest precision in terms of average absolute errors.

The absolute error results obtained using the 5th-order method are shown in Figure 7. The accuracy of the absolute errors for \(y_{1}(t)\), \(y_{2}(t)\), \(y_{3}(t)\), and \(y_{4}(t)\) corresponds to the blue Y-axis, while the accuracy of the absolute errors for \(y_{5}(t)\) corresponds to the red Y-axis. From the figure, we can see that the lowest precision of absolute errors for all four differential variables is maintained at \(10^{-7}\), while the lowest precision of absolute errors for the algebraic variable is kept at \(10^{-5}\). The experimental results suggest that the neural network's predicted solutions for the pendulum's DAEs system can also achieve high precision.

## 4 Summary and Conclusions

DAE systems are widely employed in various domains, including fluid dynamics, multibody dynamics, and control theory. In practical physical modeling, most DAE models are either low-index DAEs or high-index Hessenberg-type DAEs. Classical implicit numerical methods are suitable for a certain class of high-index DAEs, but they often lead to varying order reduction of numerical accuracy. Recently, a novel neural network method, DAE-PINN, has been developed for solving low-index DAEs. However, it cannot directly handle high-index systems. Therefore, this paper proposes a PINN-based approach using the Radau method to solve high-index DAEs systems. This method combines the strengths of the Radau IIA method with a neural network structure based on attention mechanisms and employs a time-domain decomposition strategy to enhance both efficiency and accuracy in solving these systems.

In this paper, two high-index systems, namely Hessenberg-type DAEs and pendulum model DAEs, are studied as examples. The research takes into account the influence of different orders in the Radau IIA methods and the activation functions on the accuracy of neural network solutions. Generally, employing higher-order Radau IIA methods enhances the neural network's generalization capability. However, through compara

Figure 6: The mean absolute errors for solving single Pendulum DAEs systems using PINN based on Radau IIA of order \(v=3,5,9,13\). The blue curves represent the mean absolute errors \(v=3,13\) on the left Y-axis, while the red curves correspond to \(v=5,9\) on the right Y-axis.
tive experiments with two examples, it is found that PINN based on the 5th-order Radau IIA method provide the highest accuracy in solving the systems. This conclusion is consistent with the notion that Radau-5 is a high-precision numerical method [3]. Further experimental results indicate that in high-index systems, the absolute errors for all differential variables maintain a minimum precision of \(10^{-6}\), while the absolute errors for algebraic variables maintain a minimum precision of \(10^{-5}\). This method's numerical accuracy surpasses the corresponding results in the literature [22] and, to some extent, surpasses the accuracy achieved by the DAE-PINN method [24]. This demonstrates that our method can directly and accurately solve high-index DAEs systems, showcasing strong generalization capabilities and offering a viable approach for high-precision solutions to even higher-index DAEs or challenging systems of partial differential algebraic equations. Furthermore, we have maintained the depth and width of the neural networks as in DAE-PINN [24] and have not delved into a detailed study of their impact on the accuracy of our method, which we will need to investigate in our future work.

## Acknowledgements

Project supported by the National Natural Science Foundation of China (Grant No. 12201144), the the Guangdong Basic and Applied Basic Research Foundation of China (Grant No. 2020A1515110554), the Science and Technology Foundation of Guizhou Province (Grant No. QKHJCZK[2021]YB015) of China, and Chongqing Talents Plan Youth Top-notch Project of China (Grant No. 2021000263).
