Closed-form Solutions: A New Perspective on Solving Differential Equations
Authors: Shu Wei
Yanjie Li
Lina Yu
Weijun Li
Min Wu
Linjun Sun
Jingyi Liu
Hong Qin
Yusong Deng
Jufeng Han
Yan Pang
Venue: ICML 2025 poster (2025)
Keywords: Reinforcement learning
Closed-form solution
Differential equations
Policy gradient.
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: A novel solver for differential equations based on reinforcement learning

Abstract: The quest for analytical solutions to differential equations has traditionally been constrained by the need for extensive mathematical expertise.
Machine learning methods like genetic algorithms have shown promise in this domain, but are hindered by significant computational time and the complexity of their derived solutions. 
This paper introduces **SSDE** (Symbolic Solver for Differential Equations), a novel reinforcement learning-based approach that derives symbolic closed-form solutions for various differential equations. 
Evaluations across a diverse set of ordinary and partial differential equations demonstrate that SSDE outperforms existing machine learning methods, delivering superior accuracy and efficiency in obtaining analytical solutions.

Closed-form Solutions: A New Perspective on Solving Differential Equations
Shu Wei 12
Yanjie Li
Lina Yu 12 Weijun Li4 Min Wu Linjum Sun
Jingyi Liu' Hong Qin"
Yusong Deng 12 Jufeng Han 13
Yan Pang 14
Abstract
computational techniques (Beck, 2012; Friz et al., 2020; Aderyani et al., 2022). For example, solving linear partial differential equations (PDEs) via the superposition principle involves determining Green's functions, G(x, x'), which represent solutions to PDEs with a point source at x'. This process demands deep insights into the properties of Green's functions, especially under complex boundary conditions or in multidimensional domains (Duffy, 2015). Nonlinear PDEs, however, pose greater challenges due to their inher- ent complexity, rendering analytic solutions intractable and necessitating numerical methods. Traditional numerical methods, such as the finite volume method, are prone to discretization errors and rely heavily on meshing, which can lead to convergence issues, especially with finely granulated meshes (Moukalled et al., 2016; Li et al., 2023b). Machine learning (ML) methods have emerged as powerful tools to address these limitations, significantly enhancing the ability to model and solve complex PDEs.
The quest for analytical solutions to differential equations has traditionally been constrained by the need for extensive mathematical expertise. Machine learning methods like genetic algorithms have shown promise in this domain, but are hin- dered by significant computational time and the complexity of their derived solutions. This paper introduces SSDE (Symbolic Solver for Differen- tial Equations), a novel reinforcement learning- based approach that derives symbolic closed-form solutions for various differential equations. Eval- uations across a diverse set of ordinary and par- tial differential equations demonstrate that SSDE outperforms existing machine learning methods, delivering superior accuracy and efficiency in ob- taining analytical solutions.
Chhyics-Inffrman Neural Networks (PINNs) (Raissi et al., 2019; Cuomo et al., 2022), grounded in the universal ap- proximation theorem (Hornik et al., 1989), effectively ap- proximate solutions to diverse physical systems. However, due to the nature of neural networks, these methods often ex- hibit low computational efficiency and struggle to accurately satisfy physical constraints, resulting in compromised con- vergence accuracy, particularly for stiff or high- fighffreuuency problems (Wang et al., 2021; Rao et al., 2023). The advent of neural operator methods, such as Deep Operator Network (DeepONet) (Lu al., 2021a) and the Fourier Neural Opera- tor (FNO) (Li et al., 2023b), which learn mappings between functions, offers promising avenue to tackle these issues. Yet, these neural operator algorithms typically require ex- tensive labeled data for training. Moreover, both traditional numerical methods and deep learning approaches sacrifice the interpretability of analytical solutions. This lack of inter- pretability poses challenges in accurately grasping system dynamics and limits the ability to extrapolate beyond the computational domain, thereby constraining their applica- bility in scientific research and engineering applications.
1. Introduction
Differential equations (DEs) are foundational to mathemat- ies, physics. and the natural sciences, providing abstract models for diverse physical phenomena. For example, Pois- son's equation governs the distribution of electrostatic poten- tial, while the heat equation describes temperature variations within an object (Evans, 2022). Solving these DEs is critical for understanding and predicting the dynamics of complex physical systems (Alkhadhr& Almekkawy, 2021; Savović al., 2023).
Deriving analytical solutions requires rigorous analysis of their existence, uniqueness, stability, and behavioral proper- ties, often leveraging advanced mathematical theories and
AAnnLab, .ab, Institute of Semiconductors, Chinese Academy of Sciences, Beijing, China 'College of Materials Science and Opto- Electronic Technology, University Chinese Academy of Sci- ences, Beijing. China School of Integrated Circuits, University of Chinese Academy of Sciences, Beijing. China "School of Industry-education Integration, University of Chinese Academy of Sciences, Beijing, China. Correspondence to: Lina Yu <yulina@ scmi.ac.en. en>>, Min Wu <<ummin (emiiii@emmmac.. ac.
Recent advances in symbolic regression have enabled the derivation of interpretable solutions to DEs without requir- ing expert mathematical knowledge. through data-driven approaches. Genetic programming (GP)-based approaches
Proceedings of the 4nnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Closed-form Solutions: A New Perspective Solving Differential Equations
employ physics-regularized fitness functions to evolve sym- bolic solutions. Tsoulos et al. (2006) propose a GP method that evolves solution expression genomes modeled using Formal Grammars (HOPCROFT, 1979). Oh et al. (2023) introduce a customizable GP approach (Randall et al., 2022) that directly evolves expression trees as solutions for this task. Cao et al. (2023; 2024) enhance GP algorithms with pruning techniques and transfer learning. Similarly, Ka- mali et al. (2015) and Boudouaoui et al. (2020) develop population-based methods that emulate heuristic search be- haviors inspired by resource-seeking populations. How- ever, these approaches often yield complex solutions with limited interpretability and incur high computational costs, particularly in high-dimensional problems, where finding accurate solutions remains challenging. In parallel, alter- native paradigms have emerged. Liang and Yang propose the FEX method (2022). which uses reinforcement learning to search for weighted symbolic networks to solve differen- tial equations. Liu's Kolmogorov-Armold Networks (2025) can be used to derive symbolic solutions through network symbolization. While both methods improve fitting perfor- mance by introducing additional parameters, they frequently lack interpretability. Majumdar et al. propose an approach utilizing trained PINNs to generate datasets, subsequently applying symbolic regression to extract expressions (2023). This method falls short of adhering to physical constraints, primarily due to inherent limitations in numerical precision and error propagation. The core challenge lies in balancing three critical aspects: (1) maintaining strict adherence to physical laws, (2) ensuring computational efficiency, and (3) preserving bumen-interppettble symbolic forms. Cur- rent symbolic regression approaches either sacrifice inter- pretability for accuracy or face scalability issues in complex problems. Alternative paradigms, while offering improved fitting capabilities, fundamentally differ from our goal of deriving precise, physics-compliant symbolic solutions. To overcome these limitations, we propose a novel reinforce- ment learning (RL)-based method for deriving closed-form symbolic solutions to DEs. Our contributions are as follows:
An innovative exploration technique that enables SSDE recursively solve high-dimensional PDEs, achiev- ing superior performance across tested DE types and outperforming mainstream comparative methods.
2. Related Work
Numerical methods based on ML PINNs leverage the nonlinear representation power of neural networks to solve PDEs by embedding physical constraints directly into the training process as inductive biases (Raissi et al., 2019). This approach enables unsupervised learning relying on PDE residuals while faces several challenges. For exam- ple, their dependency physical laws may not thoroughly integrate exact physical priors and they often exhibit. inefficiency when handling with stiff problems, characterized by rapidly changing solu- tions (Wang et al., 2021; Rao et al., 2023). To address these limitations, enhancements have been developed, including the hard-coding boundary conditions (Lu et al., 2021c), and employing temporal representations (Rao al., 2023). In contrast, operator learning methods (Lu et al., 2021a; Cai et al., 2021; Li et al., 2023c;b) utilize neural networks to ap- proximate nonlinear operators, providing a supervised learn- ing approach to solving PDEs that depends on paired input- output data. While these methods offer significant potential, their substantial requirement for labeled data in scientific contexts poses significant challenge. Both PINNs and neu- ral operator learning methods yield approximate solutions to PDEs. These solutions frequently lack interpretability with respect to the underlying physical phenomena and demon- strate limited generalization beyond the specific scenarios used for training.
soft constraints for embedding
Symbolic regression methods Symbolic regression meth- ods can be broadly categorized into two main categories: search-based and supervised learning-based approaches. Search- based methods explore the solution space to iden- tify optimal symbolic expressions, encompassing traditional techniques such as genetic programming (GP) and heuristic search strategies (Forrest, 1993; Koza, 1994; Schmidt & Lip- son, 2009; Staelens et al., 2013; Arnaldo et al., 2015; Iwo & Krawiec, 2019; Randall et al., 2022; Jiang & Xue, 2023). These conventional approaches, however, often grapple with an expansive search space, leading to computationally in- tensive tasks. Moreover, they tend to produce increasingly complex symbolic expressions without corresponding per- formance gains. To mitigate these challenges, recent ad- vancements have integrated neural networks to guide and constrain the search process, enhancing efficiency (Sahoo et al., 2018; Zhang et al., 2023; Dong et al., 2024). Addi- tionally, reinforcement learning techniques, including policy gradient methods and Monte Carlo tree searches, have been employed for symbolic regression (Petersen et al., 2020;
The introduction of the Symbolic Solver for Differen- tial Equations (SSDE), an RL-based paradigm that directly derives closed-form symbolic solutions to dif- ferential equations with high interpretabiility
The development of the Risk-Seeking Constant Opti- mization (RSCO) algorithm, which accelerates con- stant optimization while preserving accuracy and con- vergence.
A novel form of symbolic solutions expressed as para- metric expressions, enhancing the applicability of SSDE for discovering symbolic solutions to DEs in any single dimension.
2

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Mundhenk et al., 2021; Sun al., 2023; Li et al., 2024). Al- though these methods exhibit slower inference times since they necessitate iterative searches during inference and the time-consuming evaluation of regressed symbolic expres- sions, they demonstrate promising performance in SRBench Black-box dataset (La Cava et al., 2021). In contrast, super- vised learning-basee -based approaches leverage large-scale pre- training on synthetic datasets to directly generate symbolic expressions (Biggio et al., 2021; Kamienny al., 2022; Li et al., 2023a). By producing expressions in single forward pass, these methods achieve significantly faster inference speeds compared to search-based techniques. A notable hybrid approach, TPSR (Shojaee et al., 2024), combines elements of search and supervised learning, yielding impres- sive performance. However, the effectiveness of supervised methods diminishes when the training data distribution di- verges from that of the target physical system, leading to reduced model accuracy.
the resulting symbolic expressions may not fully adhere to physical constraints, introducing inaccuracies due to error propagation.
3. Preliminary Studies
Differential equations equations related to unknown univariate or multivariate functions and their derivatives or partial derivatives. Con- sider bounded domain S C R", with point x (I1, ( ,.,,..,, Tn) E Ω. Let =(x) denote the unknown func- tion to be determined. For a positive integer k, Duu E R** represents all k-tth order partial derivatives of u. Given function F : Ω x R x I" x ... x R= X gn* + R. the general form of a k-th order partial differential equation (PDE) be expressed as in Eq. (I). When n. 1, the PDE reduces to an ordinary differential equation (ODE).
Differential equations (DEs) are
Flx, u(x), Du u(x), D u(x) = 0 (1)
ML for solving DEs analytically Symbolic regression's ability to uncover expressive and interpretable formulations makes it a promising tool for deriving understandable sym- bolic solutions to DEs. Genetic algorithm- based symbolic regression methods have been employed to solve DEs (Tsou- los & Lagaris, 2006; Oh et al., 2023; Cao et al., 2023). Oh et al. introduce Bingo (Randall et al., 2022), using auto- matic differentiation to assess physical residuals in DEs as a fitness function to optimize symbolic expressions. Heuris- tic algorithms inspired by swarm intelligence, such as ant colony programming (Kamali et al., 2015) and artificial bee colony programming (Boudouaoui et al., 2020) emulate resource- resource-seeking behaviors search for solutions. However, these methods are computationally intensive, particularly for high-dimensional problems, where the expansive search space poses significant challenges. To address this, trans- fer learning has been explored to reduce the search space by leveraging solutions from single-dimensional instances of nigh-dimensional PDEs (HD-PDEs) (Cao et al., 2024). Yet, this approach is limited, as most physical systems gov- erned by PDEs lack known anglee-imeesioonl equations. Liang and Yang (2022) propose reinforcement learning pipeline to guide a weighted symbolic network toward ap- proximate solution of DEs. The introduction of weights reduces interpretability and complicates the discovery of minimalist analytical expressions. Kolmogorov-Arnold Net- works (KANs) (Liu et al., 2025) offer a novel approach for deriving symbolic solutions to DEs by incorporating a physics-regularized loss function. However, converting KANs' learned representations into symbolic expressions often degrades accuracy, resulting in solutions that fail to precisely satisfy the underlying equations. Majumdar et al. (2023) perform symbolic regression on numerical so- lutions from PINNs to identify interpretable expressions. Since deriede solutions are inherently approximate,
Prominent examples include Poisson's 's equation:
Pu = f(x), - -Au(x) Σ dx7
(2)
defined on the spatial domain 12, and the heat equation:
au(x,)) t) a² t)) fxx, t),
O,
(3)
Ōt defined on the spatiotemporal domain (z, t) E EΩxT, where a is constant and T > 0. To ensure a well-posed prob- lem for analytical solutions, a PDE must be accompanied by appropriate boundary conditions (BCs) and, for time- dependent PDEs, initial conditions (ICs). BCs are typi- cally specified as B(u; x E O22 0 on the boundary 02, while ICs for dependennd problems are given as I(u; t 0, x E 22) O. Together, the differential equation, BCs, and ICs define a complete mathematical problem.
Closed-form solution The closed-form solution ū(x) to differential equation (DE) is composed of a finite com- bination of known functions, such as elementary func- tions (e-g.. polynomials, exponentials, trigonometric func- tions), special functions (e.g., Bessel hypergeometric functions), or, in some contexts, neural operators. Within a bounded domain, the solution must satisfy the DE along with its associated boundary conditions (BCs) and, for time- dependent problems, initial conditions (ICs). For example, á cosh((1)) sin(22), involving the constant e, opera- tors {x, cosh, sin} and variables 21 and x2, is candidate closed-form solution to the Laplace equation - -Au(x) 0.
4. Methodology
SSDE identifies closed-form solutions using an unsuper- vised reinforcement learning framework, as depicted in Fig.
3

---

Closed-form Solutions: A New Perspective Solving Differential Equations
bols. The likelihood of the entire sampled expression is computed as p(ū]0) = I P(Ta|Ti:((G 0), where |r] is the length of the pre-order traversal. The RNN inputs in- clude representations of the parent and sibling nodes of the token being sampled, reinforcing the model's understanding of the expression's 's tree structure.
(a) Sample closed-form solution ä from RNN
Expression Tree T
(b (bJSmplle points from spatio-temporal domain Q T
Physics-regularizes local constant optimization(L.CM To ensure the expression skeleton complies with physical constraints, we employ PyTorch's automatic differentiation method (Paszke et al., 2017) to compute the differential term associated with the constraints. The constants in á are optimized by minimizing the discrepancy between predicted values and the given conditions. For instance, in spatiotem- poral dynamical system satisfying Dirichlet conditions, data points are sampled from the spatiotemporal domain Q x T. We use nonlinear optimization algorithms such as the BFGS algorithm (Fletcher, 2000), to optimize the constants by minimize the sum of mean squared errors La. (Eq. (4), where 's-t' denotes the spatio-temporal system). The loss terms for the collocation points corresponding to the differ- ential term F (defined in Eq. (1)), boundary conditions B, and initial conditions I, where (x),, ( ) (x+, 4/ vt " Nx
Forward computation with gradients
üg xx f) net
رn (x, = = + quzyt
(x,, t0 + ogt
(e) Reward
Physical constraints
R,_1(0 (
PDE
ICs/BC's
(d)
Lx-t
Vzu. a
lann
ū(x, 0)
Local Constants Optimization
Figure 1. Algorithm overview (using spatiotemporal dynamical systems as example). (a) The RNN generates skeletons for candidate solutions. (b) Sampled points are fed into candidate skeletons to construct computational graphs with constants parameters. (c) Physical constraints are built via automatic differ- entiation. (d) Constants are optimized minimize Ces. (e) The evaluator computes rewards based on physical constraints to train RNN. The process iterates until valid solution is found.
{x). 1 represent the respective collocation points, and X denotes the relative weight of each loss term.
Cat = AOMSE + AIMSE + A2MSEE.
MSE N₂ EN |F(x),,t))1
(4)
MSEB = No
1
E, No lū(zf, (t)) -
MSFI NI Σ
Nu |M(ab 0) ub/
1. To accelerate convergence while maintaining precision, we propose a risk-seeking constants optimization method. Additionally, for solving high-dimensional partial differen- tial equations (HD-PDEs), we develop recursive explo- ration technique. Pseudocode for the proposed algorithms is included in the appendix A.
Solution Evaluator We conceptualize the task of discov- ering closed-form solutions Markov Decision Process (MDP) comprising the elements (S, A, P, R). The siblings and parent nodes of the current symbolic node act as the observation S. The policy is defined by the distribution of tokens produced by the RNN, p(T:()). from which the action (token Ti) is sampled, leading to a transition to a new expression state. Each generation of the expression is guided by terminal and undiscounted reward functions, where each episode reflects one cycle. To evaluate the gen- erated expression under physics regularization, we apply squashing function to the average of the root mean square errors (RMSE) physical constraints across different sys- tems. The reward for the expression skeleton á in a spatio- temporal dynamical system is computed as:
4.1. Mhysice-Regularizer Reinforcement Learning
Expression skeleton generator We leverage symbolic ex- pression trees, where internal nodes denote mathematical op- erators, and terminal nodes represent input variables or con- stants. The pre-order traversal of the corresponding expres- sion tree represents the expression skeleton of à. Each con- stant within the skeleton a column vector of length |Caml- Therefore, when (caml + 1, a single expression skeleton can yield multiple outputs for identical inputs, as detailed in Sec- tion 4.3. Inspired by automated machine learning (AutoML) works (Bello et al., 2017; Cuomo al., 2022), employ Recurrent Neural Network (RNN) to generate categorical distributions over tokens from given symbolic library L. For each token Ti (e.g., +, X, 1 .... ), we sample with likelihood P(T,|Tia((-13 6) in an autoregressive manner to ensure contextual dependencies on previously chosen sym-
R Average(VMSE. VMSE
MSE
Training RNN using policy gradients Standard policy gradient methods typically focus on optimizing the average

---

Closed-form Solutions: A New Perspective Solving Differential Equations
performance of the policy, which deviates from the objective of finding an optimal closed-form solution. To maximize the best-case performance, we employ risk-seeking policy gradients with genetic programming introduced in (Petersen et al., 2020; Mundhenk et al., 2021), with the learning ob- jective nna (0; e) parameterized by as:
space. The pre-order traversal length |r]. and the size of the symbol library [ result in a combinatorial search space of size |C|IT), The computational complexity of solving HD-PDEs also grows with dimensionality, due to the in- creased number of variables in the physical constraints. To mitigate this, we propose a recursive exploration scheme that decomposes the solution of multidimensional differen- tial equations into sequential, dependent rounds, with each round focusing on a single variable.
Jrikk (B: e) ss Ea-p(ajo) ( 2 (0) (6)
where (f) is the reward distribution's (1 -)-uunttiee under the current policy. To foster exploration, we integrate the hierarchical entropy of the sampled expressions into the reward term, weighted by (Landajjul et al., 2021):
Parametric expression We focus on finding the expres- sion skeleton for the single variable {m=))_1 The terms formulated using other variables x- are treated as fluc- tuating parameters a(x_;) )) E E., where m denotes the number of such parameters. The skeleton û(x) is trans- formed into a parametric expression ū(x;; ā(x=:)). For the first l observed variable, we can get the solution a paramet- ric form: á ū(I1, "2,'" , T. (x (x.((......)))) Here, x (1. represents all variables except 22, Zg.
F
H(00- AxBBrp((i)))
Σ
(7)
4.2. Risk-Seeking Constant Optimization (RSCO)
Optimizing constants within the expression skeleton to sat- isfy physical constraints, particularly those related to differ- ential equations (e.g., MSEy), ), incurs considerable compu- tational costs. Building on the risk-seeking policy gradients introduced in (Petersen et al., 2020), the training objective focuses solely on the top fraction of episodes. Note that the skeleton must satisfy deterministic conditions meet physical constraints. As i approaches these deterministic conditions, the likelihood of satisfying the differential equa- tion's 's physical constraints increases. Therefore, we propose a seek-nc constant optimization method.
Recursive call When A = 0, the close-form solution is viewed as a parameter û = uy = au. For each round k >> O, we update the expression for the k-th observed variable zg. The last round's parametric expressions á are viewed as new parametric expressions ün dependent 2 and the new parameters d% 5x (x__._.)))- Namely, ύ (zx; (x which forms typical re- cursive relation.
Recursive case We ascertain the parametric skeletons üe on the single dimension, as well as the parameter vector ag by the reinforcement learning method proposed above. In the parametric expression, each parameter is vector of length [caml. corresponding to the varying values calculated from data points collected within the x_; space. Here, [Caim] represents the number of points. These parameters are opti- mized in accordance with the specified physical constraints, and the resulting parametric expressions are evaluated. If the variance of the |caml dimensional parameters falls below predetermined threshold. we surmise that this parameter does not vary with x-; and it degenerates constant. This approach allows us to distinguish between fixed constants and di parameters.
In each iteration, the constants of a batch of sampled expres- sion skeletons are optimized based on deterministic condi- tions by minimizing the modified loss C., (defined in Eq. (8)). The reward R for each expression is still computed by Eq. (5). The top fraction of samples incorporates all phys- ical constraints as defined in Eq. (4), while the constants are refined and the rewards are recalculated to yield more accurate evaluations. The accurate rewards of the precise top e fraction are then used calculate the policy gradient VoJJisk (f,c) c) via Eq. (9), which will guide the RNN. By en- forcing boundary or initial conditions before optimization, this strategy eliminates unnecessary constant optimization based on differential relationships for expression skeletons that already violate physical constraints, echoing the philos- ophy of hard constraint embedding. (Lu et al., 2021c).
Base case When the identified parametric skeleton has no remaining parameters, ã is considered a symbol expression independent of the other variables. In the final round, the parameter expression from the previous round, Æ be- comes a symbolic expression dependent only on the last variable Iz. Further, the closed-form solution to a d-dim DE is given as:
C. > M + MrMS
VoJriss (0,c) -E E.g - (1)) (R(()) (0))-
V (() 2 (0)
(9)
4.3. Recursion-Based Exploration
u (z11, T2,"" ... rkk āi (x (1.2.....))) ã (I (225............................ ( 7))
As the solution dimensionality increases, the symbol library expands, leading to an exponential growth in the search
5

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Ultimately, by substituting all parameters with their corre- sponding symbolic expressions regarding the other variables and further optimizing the constants within the expression skeleton based on physical constraints, we can derive closed-form solution for the multidimensional differential equation.
completely recover the solution in 20 independent experiments. This metric quantifies the symbolic and numerical differences between the regressed expres- sion and the true solution.
5.2. Baselines
The recursive process enables SSDE learn in focused manner, addressing one variable dimension at time.
We evaluate the performance of the proposed recursive ex- ploration method by comparing it with the following meth- ods used for seeking the closed-form solutions:
5. Experimental Settings
PINN+DSR: Performs symbolic regression with DSR (Petersen et al., 2020) on the numerical solutions ob- tained from PINNs (Raissi et al., 2019).
We present experimental settings used to evaluate SSDE by answering the following research questions (RQs):
How does SSDE compare to other mainstream methods in finding closed-form solutions to differential equations?
RQ1:
Kolmogorov-Aroold Networks (KAN): A recently proposed method designed replace Multi-Layer Per- ceptrons (MLP) (Liu et al., 2025). Compared to MLP, KAN performs better in symbolic formula represen- tation tasks (Yu et al., 2024) and shows promise for discovering closed-form solutions to HD-PDEs.
RQ2:
Can RSCO and the recursive exploration strategies improve SSDE's ability to find closed-form solutions?
Uses genetic programming to search for a function that satisfies the given PDE and boundary conditions (Oh et al., 2023).
PR-GPSR:
Does SSDE produce more accurate solutions than performing symbolic regression methods based on numeri- cal solutions?
RQ3:
All details regarding datasets, computing infrastructures, baselines and their parameter settings, along with additional comparative experiments and results, reported in ap- pendix.
5.3. Benchmark Problem Sets
To evaluate the performance of our method across a diverse range of differential equations, we selected the following equations for testing:
5.1. Metries
2D and 3D Poisson's Equations: As a classic exam- ple of -nee-nndependent partial differential equation, Poisson's equation widely used fields such as elec- trostatics, fluid dynamics, and mechanical engineering.
The evaluation of our methodology and baselines is ducted from two distinct perspectives:
The mean Root Mean Squared Error (RMSE) between the predicted closed foreddformmm solution and physical con- straints, denoted by LPHY. Taking the spatiotemporal dynamical system as an example and referring to Eq. (4), it is defined as:
2D and 3D Heat Equations: These equations are canonical examples of spatio-temporal dynamic sys- tems, describing the phenomena related to heat transfer.
2D and 3D Nonlinear Wave Equations: Nonlinear wave equations are fundamental models for describing the propagation of waves in various media, where the wave amplitude depends on both space and time.
CPHY = Average / MSEE + /MSEg- + V MSEE)
(10)
This measure assesses whether the identified solutions are consistent with physical constraints.
The equations and their deterministic conditions are pro- vided in Table
The complete recovery rate of flosed-form solution expressions. For given closed-form solution u, the skeleton of the found expression å after simplification is exactly the same equivalent to u, and after the final optimization of constants, the mean square error between the expression and the true solution is less than x 10 8 the expression is considered to have completely recovered the true solution. PRE represents the proportion of fosed-form solution expressions that
6. Results on Benchmarks
(RQ1) Solution accuracy on benchmarks Table 2 presents comparative analysis of solutions from SSDE against other methods used to derive closed-form solutions. SSDE consistently outperforms the other methods across all benchmarks, discovering solutions with significantly smaller
6

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Table Summary of tested differential equations with their deterministic conditions
NAME
PDE
DETERMINISTIC CONDITIONS
POISSON2D
d"y 6²v -
dul Bug
POISSON3D
d'v 6dv 80
duj Bug
Bug
HEAT2D
Bu
Bt
ozf
pu 82
0xg
= 7..02 [,
-
du
0
8 u
8²u
- 30z + 7..352 -277
HEAT3D
0t Ou?
0ug
02g
8'u
WAVE2D
01 def 0mg
e-t f te 3r 1.64 sin (2>) (in_)))))))((()xxxe__e (x)))
ph (gy d'u d'u 0'y
) ((mr 11 +2.75)
=(X.))=exp(oi+ 000 - ,<(00
WAVE3D
de ozf dzg
+ - u.oc 4 u +22g ccos(22)) (xx)
Table Comparison of average Ceuy (i) and recovery rates among SSDE and baselines for solving differential equations. 95% confidence intervals are obtained from the standard error between mean LruY (t) each problem set.
SSDE
PINN+DSR
KAN
PR-GPSR
NAME
CeHY 1
PRE 1
CruY +
PRE f
PRE f
LruY
PRE 1
POISSON2D
7.20E-5+2.57E-6
100% 5.71E-1#7.88 88E-2
0%
7.20E+0=1 146++0
0%
1.90E 222.322-3
0%
POISSON3D
5.70E-5#1.33E-6
100% 2.60E+0=1.11E--
0%
6.15E+0#1.52E+0
0%
1.80E-1#1.768
0%
HEAT2D
1.47E-6*7.64B-7
100%
3.21E+0=1.35E-1
0%
1.16E+1=1.24B2
0%
2.15E-1#1.28E:
0%
HEAT3D
1.05
100%
2.
0%
1.15=+121.72:+40
0%
8.87E-142.96F-1
0%
WAVE2D
4.25
100%
1.24
2.04T+0+1142mm-
0%
4.24m-143.69F:
0%
WAVE3D
8.161-621.56F--6
100%
1.34
0%
1.36F4 1+2.381+0
0%
9.87R-1+2.46f.
0%
navigate the exponentially expanded search space to find ac- curate closed-forn solutions. In contrast, SSDE's recursive exploration strategy, which employs reinforcement learning independently for each dimension, exhibits a remarkable advantage in solving higher-dimensional PDEs. Specifically, SSDE reliably identifies correct closed-form solutions for both linear and nonlinear differential equations. We also provide detailed illustration of how SSDE progressively isolates each variable and recursively constructs the correct solution through stepwise regression in the appendix
1.0
0.8
Recovery rate
0.6
0.4
(RQ2) Ablation studys We conducted series of abla- tion studies to evaluate the individual contributions of key components, including RSCO and the recursive exploration strategy (RecurExp). Figure 2 illustrates the complete re- covery rate of SSDE on the PDE benchmarks under the same execution time for each ablated configuration. The results highlight the critical role of these components in enhancing overall performance. Specifically, the absence of RSCO significantly impairs SSDE's ability to identify cor- rect closed-foma solutions within the same time frame. Simi- larly, disabling the recursive exploration strategy drastically diminishes SSDE's capability to solve high-dimensional partial differential equations, preventing it from discovering accurate closed-for form solutions.
0.2
0.0
Ppisson2D
Poisson3D
HeatD)
Heat3D
Wave3D
Benchmarks
SSDE
w.o RSCO
w// RecurExp)
Figure Ablation study of SSDE on benchmarks of PDE.
LPHY values. The increased dimensionality poses a substan- tial challenge for the other methods, as they struggle to

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Table 3. Identified Solutions by SSDE on DEs
NAME
GROUND TRUTH
IDENTIFIED SOLUTION
POISSON2D
x2. -1..21 00.5))
x (2.5000z) -13300 (1000000000000000021001000000000
POISSON3D
2.5 -
2.5000z -1.3000.1 1.1.3000..--04499944222222.....
HEAT2D
0.5r +2.z -
0. 0.5000t +log0e===oo(22200 +
HEAT3D
-1.74+-2.52 -.3x
WAVE2D
exp(2f)sin(z2)e 0.50
exp( exppf-0..000+++<<)ssnn..... 2x01--00000++1.0002++==)0ooooo
WAVE3D
(o0(=f==33) 0.56
Table 4. Comparison of average Lavy( (ū) and recovery rates between SSDE and PINN+DSR for Poisson's 's equations within I' dataset.
95% confidence intervals are obtained from the standard error between mean CHyY (û) each problem set.
SSDE
PINN+DSR
NAME DIFFERENTIAL EQUATION
PRE
CPnY
PRE
T1
=2
8.70E-8+7.221-9
100%
2.931-5+2.36F--
100%
T2
-12r+++6112
3.65E-7+2.62B-7
100%
6.55E-5+5.788
100%
T3
12==00++12
3.988-749.221-88
100%
2.521-443.388
100%
T4
8- 400f++00u0 12 +2r
7.24E-747.79E-8
100 %
2.47E-347.833
0%
Ts
am = (ran(rrn)--nsss(())))ooso///op/ -1x)
1.39E-642.06E-6
100%
8.66E-5+4.0SE:
100%
sin sin(2 )
To
7.60E-7+8.75B-7
100%
1.26E-5+8 528.588-
100%
1'r
8-u
4rf/ (
4.06E-615.82E-66
100%
1220-321.15B
0%
Ts
-025///1
1.238-6#2.011-
100%
1.98E-5#1 11776-
70%
6x1
To
8" duu 20co((13
1 1r33in(13)) sin(zz)
3.228-743.277
100 %
5.155-3#1.322 .32B-3
0%
T10
aaa
+ 8au
(sin((1) (308
1.272-6211.22m--6
100 %
2.588-323. 11..31
8x
Dez"
In
8a 8"u
= ) 22 1)/u²
1.248-542.29E-5
100 %
1.89E-342.69B.
0%
F12
8x,
8
+ d
8f
=12xi-62:11
1.228-5+2.32E-5
100%
6.68E- 1+8-1
0%
7. (RQ3) Compared with SR on Numerical Solutions
SSDE bypasses the reliance on numerical surrogates entirely and consistently recovers accurate closed-form expressions across all test cases.
To further demonstrate the superiority SSDE over sym- bolic regression methods that operate on numerical solu- tions, we conduct series of benchmarks based on the Nguyen dataset, adapted for Poisson's 's equations with known source terms. Fach differential equation is designed so that its analytical solution corresponds to one of the expressions in the Nguyen dataset (Uy et al., 2011). The baseline method first computes numerical solutions for these equations us- ing PINNs and then applies DSR directly to the numerical outputs.
8. Conclusion
We propose SSDE, reinforcement Saaning-based frame- work for discovering closed-fome symbolic solutions to differential equations. SSDE employs a recurrent neural network to generate symbolic candidates, guided by an evaluator grounded in physical constraints. To improve learning efficiency and convergence, we introduce risk- seeking constant optimization technique. By formulating symbolic solutions parametric expressions, SSDE can de- compose high-dimensiona PDEs into recursively solvable, dingle-dimensiona components. This recursive formulation allows SSDE to efficiently recover closed-form solutions to complex differential equations. Extensive experiments on various types of differential equations demonstrate that SSDE discover accurate symbolic solutions without
Table 4 reports the results of SSDE and the baseline method in terms of LPHY and recovery rate PRE- Due to the ap- proximation errors inherent in numerical solutions and the compounding of these errors during the symbolic regres- sion process, the baseline methods occasionally recover the correct solution skeletons, but frequently yield signifi- cant inacouracies during constant optimization. In contrast,
8

---

Closed-form Solutions: A New Perspective Solving Differential Equations
prior mathematical knowledge, offering a promising new direction for analytical DE solving.
Biggio, L., Bendinelli, T., Neitz, A., Lucchi, A., and Paras- candolo, G. Neural symbolic regression that scales. In International Conference Machine Learning, pp. 936- 945. PMLR, 2021.
Software and Data
Boudouaoui, Y.. Habbi, H., Ozturk, C., and Karaboga, D. Solving differential equations with artificial bee colony programming. Soft Computing. 24: 14:17991-18088 1-80077 2020.
To facilitate reproducibility and further research, we provide the complete implementation of SSDE, including training scripts, benchmark configurations, and symbolic expression evaluation tools. The source code is publicly available at:
Cai, S., Wang, Z., Lu, L.,, Zaki, T. A., and Karniadakis, G. E. Deepm&mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks. Journal of Computational Physics, 436: 110296, 2021.
https: / //gttttt.bbb//g/ com http://innonennnnSD.10
In addition, all benchmark datasets constructed for this study, including the -ighdimmnsional PDEs and Poisson problems derived from Nguyen expressions, are included in the repository. The code is implemented in Python and relies on Py Torchh and standard scientific computing libraries.
Cao, L., Zheng, Z., Ding, C., Cai, J., and Jiang. M. Genetic programming symbolic regression with simplification- pruning operator for solving differential equations. In International Conference on Neural Information Process- ing. pp. 287-298. Springer, 2023.
Acknowledgement
This work was supported by National Natural Science Foun- dation of China (No.92370117) and the CAS Project for Young Scientists in Basic Research (No. Ko.YYBBR-09))
Cao, L., Liu, Y., Wang, Z., Xu, D., Ye, K., Tan, K. C., and Jiang, M. An interpretable approach to the solu- tions of hign-dimensional partial differential equations. In Proceedings of the AAAI Conference on Artificial In- telligence, volume 38, pp. 2064002000000 20648, 2024.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Cuomo, S., Di Cola, V. S., Giampaolo, F,, Rozza, G., Raissi, M., and Piccialli, F. Scientific machine learning through Poyysiss-iffrmen neural networks: Where we are and what's ss next. Journal of Scientific Computing, 92(3):88, 2022.
References
Dong, J., Zhong, J., Liu, W.-L.,, and Zhang, J. Evolving equation learner for symbolic regression. IEEE Transac- tions on Evolutionary Computation, pp. 1-1, 2024. doi: 10. 10.1199/TEVC.2024444444
Aderyani, S. R., Saadati, R., O'Regan, D., and Alshammari, S.S. Existence, uniqueness and stability analysis with the multiple exp function method for npdes. Mathematics, 10(21):4151, 2022.
Duffy, D. G. Green's 's functions with applications. Chapman and Hall/CRC, 2015.
Alkhadhr, S. and Almekkawy, M. A combination of deep neural networks and physics to solve the inverse problem of burger's 's equation. In 2021 43rd Annual International Conference of the IEFE Engineering in Medicine & Biol- ogy Society (EMBC), pp. 4465-4468. IEEE, 2021.
Evans, L. C. Partial differential equations, volume 19,
American Mathematical Society, 2022.
Fletcher, R. Practical methods of optimization. John Wiley & Sons, 2000.
Amaldo, I., O'Reilly, U.-M., and Veeramachaneni, K. Build- ing predictive models via feature synthesis. In Proceed- ings of the 2015 annual conference on genetic and evolu- tionary computation, pp. 983-990, 2015.
Forrest, S. Genetic algorithms: principles of natural selec- tion applied to computation. Science, (51233)) 1993.
Beck, M. A brief introduction to stability theory for lin- ear pdes. In SIAM conference Nonlinear Waves and Coherent Structures, 2012.
Friz, P. K., Nilssen, T., and Stannat, W. Existence, unique- ness and stability of semi-linear rough partial differential equations. Journal of Differential Equations, 268(4): 1686-1721,
Bello, I., Zoph, B., Vasudevan, V., and Le, Q. V. Neural optimizer search with reinforcement learning. In Inter- national Conference on Machine Learning, pp. 459-468. PMLR, 2017.
HOPCROFT, J. Introduction to automata theory. Languages, and Computations, 127, 1979.

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Hornik, K., Stinchcombe, M., and White, H. Multilayer feedforward networks universal approximators. Neu- ral networks, 2(5):359-36,, 1989.
Li, Z.., Huang, D. Z., Liu, B., and Anandkumar, A. Fourier neural operator with learned deformations for pdes general geometries. Journal of Machine Learning Re- search, 24(388): 2023b.
Iwo, B. and Krawiec, K. Solving symbolic regression prob- lems with formal constraints. In Proceedings of the Ge- netic and Evolutionary Computation Conference, pp. 977- 984, 2019.
Li, Z., Meidani, K., and Farimani, A. B. Transformer for partial differential equations' operator learning. Transac- tions on Machine Learning Research, 2023c. ISSN 2835- 8856. URL https: Ic-EPP qt 3uERT.
httphhts///ppeeeeeeeeeeee /ttp:////ppeeeeeeeeeeeeee net/ forum?
Jiang, N. and Xue, Y. Symbolic regression via control variable genetic programming. In Koutra, D., Plant, C., Gomez Rodriguez, M., Baralis, E., and Bonchi, (eds.), Machine Learning and Knowledge Discovery in Databases: Research Track, pp. 178-195, Cham, 2023. Springer Nature Switzerland. ISBN 978-3-031-43421-1.
Liang, S. and Yang, H. Finite expression method for solving high-dimensional partial differential equations. arXiv preprint arXiv:2206. 10121, 2022.
Liu, Z., Wang. Y., Vaidya, S., Ruehle, F., Halverson, J., Soljacic, M., Hou, T. Y., and Tegmark, M. KAN: Kolmogorov-arnolo networks. In The Thirteenth In- lernational Conference on Learning Representations, 2025. URL https:: /ttp:////pppeeeeeeeeeeeee net/ forum? i id=0zo7qd5vzl.
Kamali, M., Kumaresan, N., and Ratnavelu, K. Solving dif- ferential equations with ant colony programming. Applied Mathematical Modelling, 39(10)3150-33633 2015. ISSN 0307-904X. doi: https///bbioog//00000/////ppmppp
003. URL http://ctg////w...ccenneeeeeeeee com/ stteece//rttccee/p/////000033300
Lu, L.,, Jin, P., Pang. G., Zhang, Z.,, and Karniadakis, G. E. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218- 2222,,,22 2021a.
Kamienny, P.-A., d' Ascoli, S., Lample, G., and Charton, F. End-to-end symbolic regression with transformers. Advances in Neural Information Processing Systems, 35: 10269-10281, 2022.
Lu, L.,, Meng, X., Mao, Z., and Karniadakis, G. E. Deep- XDE: A deep learning library for solving differential equations. SIAM Review, 63(1):208-228, 2021b. doi: 10.1137/19M1274007.
Koza, J. R. Genetic programming means for program- ming computers by natural selection. Statistics and com- puting, 4:27-112, 1994.
Lu, L., Pestourie, R., Yao, W., Wang, Z., Verdugo, F., and Johnson, S. G. Physics-informed neural networks with hard constraints for inverse design. SIAM Journal on Scientific Computing, B3 2021c
La Cava, W., Burlacu, B., Virgolin, M., Kommenda, M., Orzechowski, P., de França, F. O., Jin, Y., and Moore, J. H. Contemporary symbolic regression methods and their relative performance. Advances in neural information processing systems, 2021 2021.
Majumdar, R., Jadhav, V., Deodhar, A., Karande, S., Vig, L., and Runkana, V. Symbolic regression for pdes using pruned differentiable programs. arXiv preprint N:XY::2003.0709, 2023.
Landajuela, M., Petersen, B. K., Kim, S., Santiago, C. P., Glatt, R., Mundhenk, N., Pettit, J. F., and Faissol, D. Discovering symbolic policies with deep reinforcement learning. In International Conference on Machine Learn- ing. pp. 5979-5989, PMLR, 2021.
Moukalled, F., Mangani, L., Darwish, M., Moukalled, F., Mangani, L., and Darwish, M. The finite volume method. Springer, 2016.
Li, W., Li, W., Sun, L., Wu, M., Yu, L.,, Liu, J., Li, Y., and Tian, S. Transformer-basse model for symbolic re- gression via joint supervised learning. In The Eleventh International Conference on Learning Representations, 2023a. URL https: / wttp///nnwweeweeeee/ee/e wtt.//ouuttttto id=ULZyv9M1 j5.
Mundhenk, T., Landajuela, M., Glatt, R., Santiago, C. P., Petersen, B. K., et al. Symbolic regression via deep reinforcement learning enhanced genetic programming seeding. Advances Neural Information Processing Systems, 34:24912-24925 2021.
Li, W., Li, W., Yu, L., Wu, M., Sun, I.,, Liu, Li, Y., Wei, S., Deng, Y., and Hao, M. A neural-guided dynamic sym- bolic network for exploring mathematical expressions from data. In Proceedings of the 41st International Con- ference on Machine Learning, pp. 28222-28242. 2024.
Oh, H., Amici, R., Bomarito, G., Zhe, S., Kirby, R., and Hochhalter, J. Genetic programming based sym- bolic regression for analytical solutions to differential equations, 2023. URL https: /arxivvvvcc...////.aaaaa/aaa/vva org/ /bs// .330200177
10

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z.,, Lin, Z.., Desmaison, A., Antiga, L., and Lerer, A. Automatic differentiation in pytorch. NeurlPS Au- todiff Workshop. 2017. URL https: htttt:::///peeevieeee net/forum d=BJJs
Takamoto, M., Praditia, T.,, Leiteritz, R., MacKinlay, D., Alesiani, F., Pflüger, D., and Niepert, M. Pdehench: An extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems, 35: 1596-1611, 2022.
Petersen, B. K., Larma, M. L., Mundhenk, T. N., Santi- ago, C. P., Kim, S. K., and Kim, T. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, October 2020.
Tsoulos, I. G. and Lagaris, I. F. Solving differential equa- tions with genetic programming. Genetic Programming and Evolvable Machines, 7:33-54, 2006.
Uy, N. Q., Hoai, N. X., O'Neill, M., McKay, R. I., and Galván-López, E. Semantically-bart -aased crossover in ge- netic programming: application to real-valued symbolic regression. Genetic Programming and Evolvable Ma- chines, 12:91-119, 2011.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. Physics- informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Wang. S., Teng, Y., and Perdikaris, P. Understanding and mitigating gradient flow pathologies in thysics-informe neural networks. SIAM Journal on Scientific Computing, 4365):A3055-A3081 2021.
Randall, D. L., Townsend, T. S., Hochhalter, D., and Bomarito, G. F. Bingo: a customizable framework for symbolic regression with genetic programming. In Pro- ceedings of the Genetic and Evolutionary Computation Conference Companion, pp. 2282-22822222 2288, 2022.
Wolfram Research, Inc. Mathematica, Version 14.2, 2024. URL https: / /www. wollram. com/ ma mathemation Champaign, II., 2024.
Rao. C., Ren, P., Wang, Q.. Buyukozturk, O., Sun, H., and Liu, Y. Encoding physics to learn reaction-diffusion processes. Nature Machine Intelligence, 5(7):765-779, 2023.
Yu, R., Yu, W., and Wang. X. Kan or mip: A fairer com- parison. CoRR, abs/2407.16674, 2024. URL https: //doi. oddii.org/00) A0w448550////a
.66774.
Sahoo, S., Lampert, C., and Martius, G. Learning equations for extrapolation and control. In International Conference on Machine Learning, pp.
Zhang, M., Kim, S., Lu, P. Y., and Soljačić, M. Deep learn- ing and symbolic regression for discovering parametric equations. IEEE Transactions on Neural Networks and Learning Systems, 2023.
PMLR, 2018.
Savović, S., Ivanović, M., and Min, R. A comparative study of the explicit finite difference method and physics- informed neural networks for solving the burgers' equa- tion. Axioms, 12(10))982, 2023.
Schmidt, M. and Lipson, H. Distilling free-form natural laws from experimental data, science, 324(5923):81-88) 2009.
Shojaee, P., Meidani, K., Barati Farimani, A., and Reddy, C. Transformer-bast planning for symbolic regression. Advances in Neural Information Processing Systems, 36, 2024.
Staelens, N., Deschrijver, D., Vladislavleva, E.,, Vermeulen, B., Dhaene, T., and Demeester, P. Constructing a no- reference h. 264/ave bitstream- based video quality metric using genetic programming-basse symbolic regression. IEEE Transactions Circuits and Systems for Video Technology, 23(2):1322-1333 2013.
Sun, F., Liu, Y., Wang, J.-X., and Sun, H. Symbolic physics learner: Discovering governing equations via monte carlo tree search. In The Eleventh International Conference on Learning Representations, 2023. URL https:// openrevicn www.:eeee///tttoeem wtt.orrummmmmitttseooo 2.

---

Closed-form Solutions: A New Perspective Solving Differential Equations
A. Pseudocode for SSDE
In this section, we present the pseudocode of SSDE. The overall procedure is described Alg. 1, and the recursive exploration mechanism is shown in Alg. Additionally, we provide the subroutines SampleSkee ekeleeton on (line 5 in Alg. 1), ReplaceP ar arameters (line 12 in Alg. 3), and SingleDimensionale www.eeeniiinnnnsse...... (line in Alg. 3) in Algs. 2, 4, and 5, respectively.
Alg. I describes the overall SSDE framework with the risk-seeking constant optimization (RSCO) strategy. An RNN generates candidate symbolic skeletons, and the constants within each skeleton are optimized using boundary and initial conditions. A reward is then computed based on physical constraints. After collecting a batch of rewards, the (1 - -qquantie computed to select top-performing skeletons. These selected expressions are further refined by re-optimizing constants, and the updated rewards are used compute risk-seeking policy gradients and entropy regularization for RNN parameter updates. This process continues until the termination criterion met. Alg. 2 defines the symbolic skeleton sampling process using the RNN. The SapplyConstaalle
function (line 8) filters infeasible tokens to reduce the search space, while Arity (line 1) returns the argument count for each token. Parent srrentsbbll i ng (line 12) calculates the parent and sibling nodes to guide the next token generation. This top-down sampling process follows a similar mechanism as DSR (Petersen et al., 2020); readers are referred to the original paper for further details.
In Alg. 3, we propose a recursive algorithm based on single-dimension SSDE (described in Alg. 5) to identify the closed-form solutions high-dimensional partial differential equations (HD-PDEs). For an HD-PDE with d variables, the algorithm iteratively represents the symbolic solution as hierarchical parametric traversal. At iteration k; the traversal corresponding the k-th dimension is denoted as (zx) āg) about single variable 2& and containing parameters d where d contains parameters dependent on subsequent variables x (..,,)-- Initially, the entire solution is regarded as
Algorithm SSDE: Reinforcement Learning Framework for Discovering Closed-Form PDE Solutions with RSCO
Input: differential equation F; deterministic conditions (boundary conditions B, initial conditions I)
Parameter: RNN learning rate vy; entropy coefficient Ax; risk factor e; batch size N; sample epochs M; reward function
Output: identified closed-form solution ü' that best satisfies the physical constraints
1: Initialize RNN with parameters 0
2: Randomly sample the domain to build the dataset D = {x}))
ot
No U
(x6
3: for + 1 to M do
4: for < to N do
5: û, - SampleSkeleton(9)
Sample skeleton
6:
(L(x); + BIConnttaintsaall((,,
Calculate the peguaie-regllarized loss concerning B, I > Optimize constants in the skeleton using BFGS Compute the reward
7:
uj < DF
8:
A(a)) < RewardCalc (1))
9:
end for
10:
R (A(1)))
Collect batch rewards
11:
R - -c--uant of R
D. Compute reward threshold > Select expressions with rewards exceeding the c-quantiee threshold
12:
(; + ( R(a,)) 2 R.}
13:
forq 1 do
14:
(Lu) PrLoss (ūg. D)
> Calculate the physics-regularized loss > Optimize constants in the skeleton using BFGS
15:
a + BFGS( &g D. (Cr)g)
16:
R(i,) +- RewardCalc (
Compute precise reward
17:
end for
18:
R. + R
D Collect precise batch rewards > Compute risk-seeking policy gradient
19:
ộ + ReduceMean( (R Ro pou((U)) 0))
20:
d2 + ReduceMean (Ax VoH((1)))
> Compute entropy gradient
21:
0 + r(11 + 92)
> Apply gradients
22:
if max R > R(a") then
23: ú'" +- arg max R
D Update best expression
24:
end if
25:
end for
12

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Algorithm 2 Sampling closed-form solution's ss skeleton from the RNN
Input: RNN with parameters 0: tokens library C Output: The closed-fome form solution's 's skeleton from the RNN
1: 7 7
D. Initialize empty pre-order traversal of skeleton Initialize counter for number of unselected nodes Initial RNN input is empty parent and sibling b Initialize RNN cell state to zero
2: counter <-
3: x + +emptyy|eeme /empty
4: Co 0
5: 1+1
6: while counter + 0 do
7: (1h,, c.) t RNN(I, 1)-1;0)
Emit probabilities; update state Adjust probabilities
8: 2 + ApplyConstraints 1/ L, 7)
9: + Categorical (+).
B> Sample next token D Append token to traversal Update number unselected nodes Compute next parent and sibling
10:
r ++ ー| T
11: counter + counter + Arity(n:)
x - ParentSibling(7))
13: i ++1+1
14: end while
15: û 1-77
Assemble the pre-order traversal 7 into an expression skeleton
16: return i
Algorithm Recursive Exploration based on single-dimensional SSDE
Input: Differential equation F; deterministic conditions (boundary conditions B, initial conditions Z).
Parameter: Number of variables differential equations d.
Output: Identified closed-form solution that best satisfies the physical constraints.
1: Randomly sample the domain build the dataset D {x j U (x], i U (x),, /
2: Treat the solution as parameter ag, initialize with B, I
3: ((0) (2to: (o)) + 00
Initialize empty regressed expression's traversal with on
4: |ăn! <-1
> Initialize number of parameters f(0) as I Initialize pre-order traversal of regressed skeleton with +(0)
5: + ==(0)) (roi ao)
6: for kerk + to do
7: +(K) (FA; -
8: for (m²) (zz; ax)) SingeeDDiSSDDEDDDdiaa-
to | do
9:
T
Seek parametric expressions for each dg-1 Append new parametric expression of kth variable
Fx (x + Fn (,); [-m) (2x;
10:
11:
end for
12:
T + Repalce (aableePaaammeeee////-------111 Tk) Replace parameters in with new sub-expressions discovered for
the kth variable
13:
if l: + d then (A)
14:
PattaalConttrainCCClle///aaii
Calculate constraints across the regressed dimensions Optimize parameters and constants the skeleton with BFGS
f +- DFG
else
L. PiyxicllCOnstmimnnCCCa(c/aaa.a.
Calculate the physics-regularized loss Optimize constants in the skeleton with BFGS
18: T" + BFGS(T, D, Cxr)
19: end if
20: end for
21: 0* <- 7
Assemble the pre-order traversal 7 into an expression skeleton
22:
return i**
single parameter dependent on all variables X1:d- We formally denote this as a trivial traversal f(0) (x0, ao), where xp is merely placeholder without practical meaning, simplifying to f(0) = ag- At each subsequent iteration k. Alg. 5 generates single-dimensions traversals i (zx; ax) for parameters ax-1 identified in the previous step. Alg. then substitutes the previous parameters σ within the current traversal with these newly obtained expressions, producing the updated
13

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Algorithm 4 Replace parameters in regressed skeleton with kth variable's traversals
Input: Regressed skeleton =(No:k-15 d new identified traversals F of k-th variable.
Output: New traversal т of symbolic skeleton variables.
1: 7=
Initialize empty traversal of ls variables
2: for 4--11 to [f(xxx- 1:00 -g--))))))) 1]] do
3:
for m + 1 to [ do
4:
if f = 0%-1 then
5:
T 7-4-lo =
Replace parameter or-1 with the corresponding traversal of the kth variable
6:
else
7:
T T [-
b Append regressed traversal
8: end if
9: end for
10: end for
11: return 7
traversal F involving variables. To mitigate error propagation during the recursion, constants in the traversal are refined by minimizing c(k), which incorporates partial physical constraints. At the final iteration (k = d), when the traversal involves the last variable za without further parameters, we optimize the traversal T against the complete physical constraints Cy- thereby obtaining the final precise symbolic expression.
The underlying principle of Alg. 5 closely parallels that of Alg. 1, yet specifically targets discovering symbolic traversals corresponding to linlssndi variable skeletons. Here, we employ an RNN to generate candidate traversals [Fx]] associated with the k-th variable during each training batch. Each constant token (c) within these traversals is represented as a vector [Puly of dimensionality |edim]. Initially, the vectors are optimized by minimizing deterministic constraints sampled from boundary (B) and initial conditions (I).
We denote the h-th token within traversal [Fxly |pil) thus |adim| |p If the variance of given vector [px]] falls below predefined threshold d, it is treated as a constant; otherwise, it remains as a parameter. The parameters retained in [pk]y subsequently become new parameters for traversal [Fely- As in Alg. we utilize the RSCO strategy to select the top-e traversals within the batch, computing rewards accordingly. Finally, the RNN parameters are updated using risk-seeking policy gradients augmented by entropy regularization.
As an illustrative example, consider discovering the glosed-form solution:
u = 2.5x 1.3r. + 0.5x.
Initially, we represent the symbolic solution simply û ag. The expression skeleton corresponding the first variable Ij is identified using Alg. 5, yielding a parametric traversal:
ū(z (1) q +a
where the parameter a1 depends on subsequent variables x2 and æs. Next, Alg. 5 identifies the skeleton for the second variable Nz:
0n (12; 02) cιaz + 02
Using Alg. 4, we combine the skeletons derived from the first two dimensions to form the intermediate traversal:
à 2 o2z2 02
where 0g remains dependent on Z3. Subsequently, the traversal for the final variable 23 is obtained as
(2(23) = c
Combining all traversals results in the complete symbolic expression
c + canl + cgm2
Finally, optimizing the constants C1, C22 C3 against the full set of physical constraints precisely recovers the exact closed-form solution.
14

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Algorithm 5 Single-dimensional SSDE
Input: Differential equation F; deterministic conditions (boundary conditions B, initial conditions Z); the mith last regressed variable's paramter d the order dth of the variable currently focus Parameter: RNN learning rate Y/i entropy coefficient Ax; risk factor e; batch size N; sample epochs M; reward function R; early stop threshold R5; constant judgement threshold
Output: Pre-order traversal Fi of the parametric expression of k-th variable that best satisfies the physical con- straints
1: Randomly sample the domain build the dataset D {x
M U (x 1i U (xi= uf
2: Initialize generator RNN with parameters 0
3: for < 1 to M do
4: for j < to N do
5: [rl] < SampleSkeleton(0)
D. Sample k-th variable's expression skeleton
6:
[C,,] +- SingeeDinBBCConstaans(((.... 0h [mm]
Calculate single-dimensional constraints B, I
7:
[fr], (Tx, [m..1)) < BFGS((F)); D, [mm]m 1
Optimize parameters in the skeleton using BFGS
8:
for h + 1 to | | do
9:
if ((
then
10:
[ + average((pt]))) end end for āe
>> Fix parameter to average value if its variance falls below threshold 8
11:
(p l is not const.}
Extract variable parameters from optimized constants Compute the reward
14:
R((Fel)) < RewardCalc((Fal))
15:
end for
16:
+ [(((5r]))}}-
Compute batch rewards
R. (1- e) -quantile of R
Compute reward threshold
18:
T: [:([xlg)+ [[xll R((F()) R.}
Select expressions with rewards exceeding the e-quantile threshold
19:
1 N do
20:
[Call + SingddDinnlyxiccCCoossmmmmttt/ D, [m)) (TK, [pkka) <- BFGS((F))g: D, a [m)] (g) )
>> Calculate the physics-regularized loss
Optimize constants in the skeleton using BFGS
22:
for h 1 to [g do
1
23:
if variance (1% < then
j < average [p))))
Fix parameter average value if its variance falls below threshold
25:
end if end for
26:
27:
ai (P [plly is not )onst.}
> Extract variable parameters from optimized constants Compute the reward
28:
[(FFl))) + (ewa)
29:
end for
30:
R [ Jkm
Compute precise batch rewards
31:
δ1 + ReduceMcan (R V ((T())
> Compute risk-seeking policy gradient & Compute entropy gradient
32:
82 + ReduceMean AH ((H(()(0))
33:
0 - + r()1 + (2)
> Apply gradients
34:
if max R > [( then
[7,] + [ jargmax,
b Update best expression
36: end if
37:
end for
B. Computing Infrastructure and Hyperparameter Settings
In this section, we provide additional experimental details, including the computing infrastructure and the specific hyperpa- rameter configurations for SSDE.
Computing infrastructure All experiments reported in this work were conducted on an Intel()) Xeon(R) Gold 6138 CPU 2.00GHz. The algorithm implementation can also leverage GPU acceleration for improved computational efficiency.
15

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Table 5. Hyperparameters for SSDE.
HYPERPARAMETER
SYMBOL
VALUE
LEARNING RATE
0.0010
ENTROPY WEIGHT
A
0.07
ENTROPY GAMMA
7
0.7
RNN CELL SIZE
32
RNN CELL LAYERS
1
RISK FACTOR
0.05
MAX EPOCH
M
200
BATCH SIZE
N
1000
PDE CONSTRAINT WEIGHT
Ao
1
BC CONSTRAINT WEIGHT
À1
1
IC CONSTRAINT WEIGHT
A2
1
Table Convergent relative error of baselines high-dimensional dataset
NAME
SSDE
MATHEMATICA
PR-GPSR
KAN
FEX
DSR
PINN
POISSON2D
1.55E-05
X
1,49E-02
3.24E+00
2.00E-02
9.85E-02
6.24p-04
POISSON3D
3.58E-06
2.01E-02
2.57E+00
9.13E-02
2.15E-01
3.93n-03
HEAT2D
4.45E-06
4.57E-02
6.54E+00
4.70F-02
2.19E-01
6.98E-03
HEAT3D
9.34E-06
1.39E+00
5,58E+00
3.63E-01
3.20E-01
1.15E-02
WAVE2D
3.11E-05
X
-180E-01
1.18E+00
2.34E-01
7.40E-02
1,34E-02
WAVE3D
7.62E-06
4.56E-01
7.26E+00
3.97E-01
1.95E+00
3.69E-02
The hyperparameters of our method mainly include the hyperparameters for the reinforcement learning policy gradient optimization, and weighting coefficients of the physical constraints defined the loss function Lyt
Hyperparameter settings
C.1 = AgMSEg + M,ES + MrMSE,
(11)
For the policy optimizer, we considered the following hyperparameter search spaces: learning rate E {0.0001, 0.0005, 0.0010), entropy weight λx E 0.7.4,0.000000 0.0, and hierarchical entropy's coefficient y 0.7,, 0.9). Hyperparameter tuning was performed via grid search on benchmark problems Is and F10 from the selected by minimizing the average LrHY. The optimal hyperparameters
T dataset, and the best hyperparameters identified through this tuning process summarized in Table 5, and these configurations were consistently used across all benchmark experiments.
C. Additional Experiments and Results
C.1. Experiments with Additional Baselines
To facilitate fair comparison across different methods, we adopt the relative L2 error (defined by Eq. (12)) as unified evaluation metric in our supplementary experiments.
Yime,
y 22
Mean Relative L,2 Loss
(12)
72
| 1
In the main text, we compare SSDE with the baseline method PINN+DSR. which performs symbolic regression on numerical solutions of differential equations. Additionally, we evaluate both methods suite of PDE benchmarks. To further enrich the comparison, we introduce two additional baselines. The first is Mathematica (Wolfram Research, Inc., 2024), a commercial software capable of computing analytical solutions to differential equations. The second is FEX (Liang Yang, 2022), reinforcement learning- based framework for solving differential equations. Detailed descriptions of these two baselines are provided in Appendix D. Table 6 summarizes the performance of all methods across the PDE benchmarks. SSDE consistently identifies correct closed-form solutions and outperforms all baseline methods in terms of accuracy.
16

---

Closed-form Solutions: A New Perspective Solving Differential Equations
(a) 2500
Average of 55DE
Average w./0 RSCO
2000
1500
1000
336
23
85
34
300
53
(b) 200000
SSDE
200000
150000
100000
50000
Figure 3. Ablation study RSCO efficiency (results averaged 20 benchmark runs) (a) Time efficiency comparison: RSCO achieves 4.65x speedup (ju = 934s vs. 200s) in solution discovery time. (b) Search space comparison: While requiring 22% more expression evaluations, RSCO maintains superior time efficiency.
C.2. Complexity Analysis
We conduct a complexity analysis to address two research questions (RQs):
RQ1 Does the proposed risk-seeking constant optimization (RSCO) algorithm accelerate the discovery of closed-form solutions?
RQ2 How does the size of the symbolic library affect the computational complexity of the algorithm?
To answer RQ1, perform 20 independent trials each of the T1-['s benchmarks, comparing the average convergence time and number of iterations required with and without RSCO. As shown in Figure RSCO consistently reduces the total solving time across all benchmarks, despite only marginally improving the average number of iterations. Interestingly, under the same random seed, RSCO sometimes achieves convergence fewer iterations, highlighting the benefit of prioritizing boundary condition satisfaction early in the optimization process.
To investigate RQ2, we evaluate the impact of symbolic library size on solving complexity using the Ii-I', benchmarks. As shown in Figure functions such sin, cos, exp, and log. We observe that enlarging the symbolic library significantly increases the average time and number of iterations required to recover the correct solution. Moreover, the addition of higher-order operators like exp introduces a much larger increase in complexity than trigonometric functions. This is likely due to the presence of high-order polynomial terms in the target solutions, which can be easily confused with exponential or logarithmic expressions during symbolic search, especially given the close derivative relationships among e", In x, and certain polynomial approximations.
start with a minimal operator set +, -, x, + and progressively expand it by adding transcendental
17

---

Closed-form Solutions: A New Perspective Solving Differential Equations
(a)
Library
Time
(b)
Library
b.:K//iin.cos.omm +:uxulsinccss t,.K,
250000
200000
Exprassion
150000
100000
50000
Figure Computational complexity scaling with library size. Average time and (b) explored expression count for SSDE in P1-4 experiments across 20 independent runs.
D. Details of Baselines Algorithms
In this section, we introduce the mainstream methods used for identifying the closed-form symbolic solution to high- dimensional PDEs.
Kolmogorov-Arnold Networks (KAN) newly proposed method designed to replace Multi-Layer P) incorporates capabilities for automatic symbolic regression(Liu al., 2025). Compared with MLP, performs better on symbolic formula representation tasks(Yu al., 2024) and is promising for discovering closed-form solutions HD-PDEs. Moreover, compared to EQL-((Shhoo et al., 2018) 8)), which is a fully-connected work where elementary functions are used as activation functions, this method does not require additional sparsification for symbolic regression and only requires a small number of nodes to complete the task of fitting the data. In this work, we use the open-source code and implementation provided by the authors!. The open source code has been implemented for the 2D Poisson's 's equation merely. We have also implemented the algorithm for the 3D Poisson's 's equation, imperative note that due the automatic symbolic regression process, expressions involving basic operations such as be directly derived from KAN's network structure. Therefore, the
and 3D heat equations based on the open source code.
addition, subtraction, multiplication, and division symbol library used in the KAN method for automatic regression does not include these operators. Additionally, cos token can also be directly obtained from s sin by adding a bias (i.e., sin(x + )) cos(r)). therefore, the cos token is also absent from the symbol library. To ensure a fair comparison, apart from the symbol library employed in our algorithm, we also introduced {o?.03, ) as additional symbols for use in the KAN method. In addition, we utilized the same hyperparameters and network architecture settings as those employed the KAN repository's source code on partial differential equations.
Physics-regularizeed Genetic Programming Based Symbolic Regression(PR-GGRS physics-regularized genetic programming(GT) based symbolic regression method for analytical solutions to differential equations. They use the mean-squared error with physical constraints the fitness of population measure how well
PR-GPSR((Oh et al., 2023)) is a
Phttps://gtt thub. com www.iamiin//ppec.....
18

---

Closed-form Solutions: A New Perspective Solving Differential Equations
the found expression ā satisfies the differential equations. In this work, use the open-source code and implementation provided by the authors?. The open source code has been implemented for the 2D Poisson's equation. We have also implemented the PR-GPSR algorithm for the 3D Poisson's 's equation, 2D and 3D heat equations based on the open source code. To ensure fair comparison, we use the hyperparameters settings reported in the original paper. In addition, we use the same symbol library as our algorithm and add the CONST token to give it the ability to generate constants.
To demonstrate that SSDE has better solution accuracy than the natural paradigm of directly applying symbolic regression algorithms numerical solutions, we designed this baseline. It consists of two parts: the PINN algorithm and the DSR algorithm. The PINN algorithm is used to generate numerical solutions for each equation in the T dataset, and the DSR algorithm directly regresses the corresponding symbolic expressions based on the numerical solutions. The following will introduce these two algorithms and their implementation in detail.
PINN+DSR
Physics-Informen Neural Network (PINN) A neural network method proposed by (Raissi et al., 2019), which is widely used to solve the partial differential equations numerically. PINN directly uses a fully connected neural network to approximate the solution function of the partial differential equation. Its numerical solution has high accuracy when solving wwwdimeesional PDEs with uncomplicated computational domains and boundary conditions. In this work, we use the open-source code and implementation by DeepXDE> (Lu et al., 2021b). To ensure fair and effective comparison, we train the network with PDEs and deterministic conditions defined in I' dataset until the network convergent. And sample the points in the computational domain identically to our approach used for subsequent symbolic regression. The convergent mean relative L2 error on the I' dataset shown in Table 7.
Deep Symbolic Regression (DSR) DSR((Petersen et al., 2020)) is symbolic regression method based on reinforcement learning. The risk-seeking policy gradient is used to guide the RNN to generate the symbolic expressions corresponding to the given observational data. In this work, we use the open-source code and implementation provided by the authors". To ensure fair comparison, we use the same library as our algorithm and add the CONST token to give it the ability to generate constants. In addition, we used the hyperparameter combination reported in table for the policy optimizer.
Finite Expression Method (FEX) FEX is a reinforcement learning- based method for solving high high-dimensional partial differential equations by discovering solutions expressed as finite mathematical expressions (Liang & Yang. 2022). Instead of parameterizing the solution with neural networks, FEX formulates the PDE-solving process as a combinatorial optimization problem over the space of mathematical expressions, represented as binary trees with both discrete operators and continuous parameters. While FEX theoretically designed to recover compact, interpretable expressions that approximate the true solution, in practice, its symbolic structures often coarse approximations, and the accuracy is largely achieved by tuning the continuous parameters. As a result, the discovered expressions tend to resemble parameterized symbolic networks, wherein most of the fitting power comes from weight optimization rather than symbolic structure. This compromises interpretability and makes it difficult to recover exact closed-form solutions, especially when the symbolic search fails to capture the core structure of the ground truth. Nonetheless, FEX still presents a promising approach for symbolic approximation of PDEs with lower memory cost and offers partial interpretability compared to black-box neural solvers. In this work, we use the official open-source implementation of FEX provided by the authors The repository includes support for 2D and 3D Poisson equations. To enable fair and comprehensive comparison, we further extend the implementation to support 2D and 3D heat equations and wave equations.
Mathematica Mathematica is a commercial symbolic computation software developed by Wolfram Research, widely used for solving mathematical problems involving algebraic manipulation, symbolic integration, and differential equations. provides a built-in solver for obtaining analytical solutions to wide range of ordinary and partial differential equations, leveraging rule-based rewriting systems and symbolic inference engines. In contrast data-driven approaches, Mathematica relies on curated mathematical knowledge and a comprehensive symbolic engine to derive closed-form solutions. However, Mathematica's 's performance is strongly tied to the internal heuristics and rule libraries it employs, which are generally optimized for classical PDEs with known solution templates. For high-dimensional PDEs, equations with nonstandard
Pttpps: petttps://gttin ://gtthbb...///// tenttps//gitthe wwtp///wwlxxxii///epdde. 1 i c-opt zaticn
"https://ge thub. htt.hbb.ccom//oooo////eepee symbol
hhttps://gt thub. watt.//////eenngaaaaaaaa /F iniie-ieeeeie tetteee-xxessseeee pwmppeesiio-metteet. -o-mmethod
19

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Table 7. Convergent relative La of PINN I dataset
NAME
EPOCH
MEAN RELATIVE La ERROR
Di
30000
2.25E-05
T2
30000
6.32E-06
T3
30000
7.04E-05
Ta 30000
5.28E-05
I's 30000
2.57E-05
30000 T6
2.
L'T 30000
1.08E-05
Ts
30000
1.05F-05
T9
50000
2.59F-03
T10
50000
5.91F-04
Pn 50000
3.61E-04
P12
50000
2.88E-03
Table 8. Configuration of high-dimensional PDEs dataset.
NAME
DOMAIN
LIBRARY L
GROunD TRUTH
POISSON2D
[-1, 1J2
Co U (12)
u(() = 2.5x| 1.3r] + 0.5r. 1.7x2
POISSON3D
[_,
Lo [32
u(X) 2.5x
1.3r
0.5ug
HEAT2D
0, 1-1,
Lo {x2,t)
u(X, t) =2.5x 1.3x
+.5t
HEAT3D
0,
1,
Co 12,20,,
u(X, t) -2.5x -13r%
0.52g 1.7t
WAVE2D
0, X
Co {x2, t}
u(x, t) exp(xf) (in(2))) 0.5A
WAVE3D
00.1
Co U (22, 133,
u(X, L) exp(2f + 12) cos(2)) 0.5A
Owing to the structure of the KAN network, the symbol library employed in the KAN baseline encompasses variables and constants, as well as the operators (o²,o",e sin, log, exp).
boundary conditions, synthetic benchmarks designed outside its built-in library scope, Mathematica often fails to return exact symbolic solutions may produce no result at all. In this work, we use Mathematica (Version 14.2) as a baseline to assess its symbolic solving capabilities on our PDE benchmarks.
E. Details of Datasets
E.1. Dataset High-dimensional PDEs
Despite the availability of existing datasets for validating numerical PDE solvers based deep learning methods (Takamoto al., 2022), these datasets typically do not provide explicitly defined symbolic solutions. Consequently, evaluating symbolic-solution discovery algorithms on such datasets is not only computationally expensive and inefficient but also does not guarantee the correctness or interpretability of the resulting symbolic expressions.
To comprehensively assess the effectiveness of our algorithm discovering diverse symbolic solutions, we constructed a new benchmark dataset consisting of various types differential equations (DEs). This dataset enables rigorous performance evaluation of our proposed method and other state-of-the-- symbolic-solution discovery algorithms. Specifically, each benchmark example (see Table includes the differential equation itself, the computational domain, associated boundary conditions, and symbol library defining permissible functional forms. For clarity and brevity, define symbol libraries relative to base library Lo. given by Lo = (+, x, +, sin, cos, exp, log, #1, const.), ). The notation [-1, 1] denotes spatial domain represented by a d-dimensional hypercube, while the notation [0, 1] x [-1, 1]d explicitly differentiates between the temporal domain (the first interval) and the spatial domain (the second interval). The benchmark differential equations encompass a range of scenarios, including both stationary and spatiotemporal dynamics, as well as linear and nonlinear PDEs, thereby providing comprehensive coverage for validating symbolic solution discovery methods.
20

---

Closed-form Solutions: A New Perspective Solving Differential Equations
Table Configuration of T dataset Poisson's Equation
NAME
DOMAIN
LIBRARY L
Ground TRUTH
Tí
1
Co
u(x) 2
2j
T2
[-1,0 1
Co
u(x) zt + a
+
T3
[1,
Co
u(x) 2i
+ =
-1,1
Lo
u(x) =xi +
+3
Ts
-1.1
Co
u(x) sin(2))cos(211 -1
T6
1
Co
u(x) sin(xr) sin(21 +
17
1.5
Co
u(x) log(æ1 + + 1) log
Is
0.5, 155
Lo
u(x) /21
0.5, 1.5
Co U [x2)
u(x)
T10
1.5
Co U [x2)
u(x) 2sin(z:)
In
0.5, 1.5
Co [
u(x) x²7
T12
[],
Co U {x2}
u(x) = 21 21
E.2. Benchmarks on Poisson's 's Equation
To compare SSDE with methods that apply symbolic regression (SR) directly to numerical solutions of differential equations, we constructed dataset inspired by the symbolic regression task the Nguyen benchmark (Uy et al., 2011). Specifically, we selected 12 target expressions from the Nguyen benchmark and used them as closed-fone solutions to derive corresponding Poisson equations. The complete configuration of this dataset is summarized in Table It is worth noting that expressions involving only a single variable result ordinary differential equations (ODEs), while those involving two variables give rise partial differential equations (PDEs). This setup allows to evaluate the performance of SSDE and other baselines both ODE and PDE settings in a controlled, ground-truthh-wee SSDE successfully recovers the Nguyen- 12 expressions as exact closed-form solutions, whereas conventional symbolic regression methods like DSR (Petersen et al., 2020) and DSO (Mundhenk et al., 2021) are unable to fully reconstruct the target forms.
environment. Leveraging its recursive exploration strategy,
F. Details of Experimental Results
SSDE successfully recovered the symbolic solutions for all tested partial differential equations, with the intermediate regressed expressions presented in Table 10. In contrast, the solutions obtained by PINN+DSR often omitted essential variables. We hypothesize that DSR's failure to identify correct solutions stems from two key limitations: (1) the increased dimensionality significantly expands the search space, making exploration more difficult; and (2) the policy gradients corresponding to different variables may interfere with each other during training, hindering convergence to the correct expression. Although the KAN method can approximate numerical solutions of PDEs with high accuracy through its network architecture, its symbolic regression capability remains limited. The resulting expressions tend be dense and lack interpretability due insufficient sparsity. In terms of efficiency, PR-GPSR performs poorly: in experiments, it required approximately day to complete 1500 generations. To ensure fair comparison across methods, capped the number of generations for PR-GPSR at 1500. Even within this limit, PR-GPSR consistently failed to recover the correct symbolic solutions across all benchmark equations.
21

---

Closed-form Solutions: A New Perspective Solving Differential Equations
expl -0 -0. 0.49992)
expl explas + 1.0010a7
exp( exp( L001007 0 0.49999)) (x4999())xxp(ee) -expp paatt 1.0010l (cos (22)))
1.0010ln (cos (22)))
exp(x3 exp(x3+1.ud +1 +1.0031
100012f
exp exp[1.0 1.00100 -
0. 0.49999) (.4999t)-xxplee eepp(a3
0.5000f) (#2)
WaVE3D
0. 0.4999t)
sin ( 5002
L88588os(L1 (1)-
1.0000sin( (x2))
a1e 0.30022
e - 70f-0.5002 5002e
r f -0-00.00 5.0000
WAVE2D
0gcos(
Table The stepwise regressed solutions of SSDE used to solve PDEs
100000011.300000 1.3000e
0. 0.4998 1.7000t+
2. 2.5000e 1.00001-1.3000 -122-1.30 1.30000
7177099172 22.500000
2000022 +0.490993
an 1.70000-1
at 1. 70000++++
as 1. 700004
2. 2.5000a1 -1
HEAT3D
2.5000e 1
1. 300000
0.0002 1
0.0006
log (exp(2 (og(exp(2.449656) 4995<1 22 300000
log (expe (e0eepp225000ff 500001 a (1.300022 8. 8.1811e-T))
log (exp (o2 2. 4995xx1)
0.5001
HEAT2D
0. 0.5001e*
0.500112 +
0.0008))
0. 0.5000r² +
2. 2.5000 + 0.65522ez
1.650222
2. 2.5000 1.300322
2. 2.5000 11.00:4-113000000 1.3000x +
0 4999992 +2.3 32.3763--5
-0.0003
1399993+2.7700
2.5000z1+0 +ai
POISSON3D
1.300322
0 0.500023
2f 500000 1332) +
of (250000 +
11 500000 1+3e))++
PoIsSON2D
0. 414997722 (a22
1.4005) 22
0. 1149999cc (302
1.4022) - 222
On
STEP 1
STEP 2
STEP 3
STEP 5
STEP
STEP
22


Reviews:
[{"review_id": "LXupf7BIln", "reviewer": "Reviewer_ateq", "summary": "This paper proposes SSDE, a reinforcement learning-based framework for deriving closed-form symbolic solutions to differential equations. The authors introduce a risk-seeking constant optimization technique and recursive exploration strategy to enhance the method's efficiency. Experiments are conducted on various ordinary and partial differential equations to demonstrate the approach's effectiveness.", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the feedback and suggestions, we will add clarification where needed and include suggestions as space permits.  Extra results: [here](https://anonymous.4open.science/r/SSDE-A47B/SSDE_icml2025_rebuttal.pdf). All figures/tables are from this link. \n\n**Q1: How does SSDE compare to other reinforcement learning approaches for symbolic regression [2-5]?**\n\nAdded comparisons with [2] (`FEX` in Table 3), showing SSDE’s superiority. [2] will be cited in related work. We performed parameter tuning for the PRGP baseline, which improved its performance across benchmarks. Although the GP method can generate highly complex expressions to achieve good performance on the Van der Pol benchmark, it fails to recover closed-form solutions in high-dimensional settings. [3-5] focus on differential equation discovery via symbolic regression, addressing different problem settings and thus not directly comparable.\n\n**Q2:What is the computational efficiency compared to traditional methods(Spline approximation)? Can SSDE handle nonlinear PDEs?**\n  \nSSDE targets closed-form solutions (unattainable by numerical methods like [1]), offering interpretability and mesh-free advantages. Traditional methods trade efficiency for precision via mesh refinement (non-convergence risks). \n\nWe supplement the experiments on nonlinear wave equations. The experimental results is shown in Table 2 and Table 3, which suggests that our method is competitive. Beside, we add Mathematica 14.1 (state-of-the-art symbolic solver) as the extra baseline (In Table 3). \n\n\n**Q3: Lacks convergence/complexity analysis. How does scalability depend on problem size and accuracy?**\n\n1. Convergence follows risk-seeking policy gradients[1].  \n2. The search complexity increases exponentially with the size of the library, but the RecurExp method significantly reduces the complexity caused by the increase in variables(Section 4.3). More detailed analysis will be added in the final submission. Accuracy depends on symbolic expressivity, bypassing numerical discretization limits.\n\n[1] Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. \n\n\n**Q4:How was the training dataset acquired?Can SSDE be parallelized on CPUs and GPUs?**\n\nSSDE uses a **self-supervised paradigm** without external datasets. Collocation points are sampled randomly during training. RNN is used to generate candidate expressions evaluated via PDE residual loss and BC/ICs penalties. Policy gradients are updated to ensure identified solutions are constrained by the governing equations. It supports CPU/GPU parallelization, but experiments were conducted on CPUs due to the lack of access to GPU resources.\n\n**Q5:What are the limitations of SSDE?**\n\nSSDE currently lacks efficiency in solving systems of DEs due to their inherent coupled-solution nature. Exploring the full solution space requires multiple RNNs, and evaluating reward-guided gradients for each solution necessitates novel methods. We aim to address this challenging issue in future work.\n\n**Other suggestions**\n\n1. The evaluation criteria are limited to RMSE and recovery rate, which do not provide a comprehensive assessment of the method's capabilities compared to existing approaches.\n\n   Added L2 relative error comparisons(Table 3). While recovery rate and MRMSE are tailored to our problem's unique requirements, we clarify their rationale:\n\n   Recovery rate: Assesses closed-form solutions by requiring both symbolic and numerical equivalence (<1e-8 MSE threshold).\n\n   MRMSE: Quantifies equation satisfaction when symbolic ground truth is unavailable, complementing traditional numerical metrics.\n\n2. Address mathematical notation inconsistencies and errors.\n   \n   We will thoroughly check and fix grammatical errors in the final submission. The F definition on line 153 refers to: Evans, L. C. Partial differential equations.\n\n\n3. Limited parameter settings are explored.\n   \n   We provide a full scan of the parameters in Appendix B.\n\n4. Key related works in neural operators, traditional numerical methods, and other symbolic regression approaches are not sufficiently discussed or compared. Several important references are missing.\n\n   We have discussed neural operators, traditional numerical methods, and symbolic regression approaches in the introduction/related work, emphasizing their inability to achieve closed-form solutions for DEs.  While neural operator references are included (citing their published works), we acknowledge the need for deeper engagement with additional literature. We will supplement references in the final version.\n\nSSDE targets discovering closed-form solutions for DEs—a goal distinct from numerical approximation[1] or symbolic equation discovery (e.g., [3-5]). Given clarifications on SSDE's focus and added experiments highlighting contributions, we kindly request reconsideration of your score per these revisions.", "claims_and_evidence": "The paper claims that SSDE can effectively find closed-form solutions for differential equations without prior mathematical background. While the idea of combining symbolic learning with neural networks to solve PDEs is promising, the evidence provided is insufficient to fully support this claim. The experiments are limited to linear PDEs and ODEs, which can be solved using traditional methods.", "methods_and_evaluation_criteria": "The methodology combines reinforcement learning with symbolic regression and introduces novel optimization techniques. However, the evaluation criteria are limited to RMSE and recovery rate, which do not provide a comprehensive assessment of the method's capabilities compared to existing approaches.", "theoretical_claims": "The theoretical claims about the method's ability to handle high-dimensional PDEs are not sufficiently supported. The paper lacks a rigorous theoretical analysis of the algorithm's convergence properties and computational complexity.", "experimental_designs_or_analyses": "The experimental design has several limitations:\n\n1. Only linear PDEs and ODEs are tested\n2. No comparison with traditional numerical methods\n3. Limited parameter settings are explored\n4. No analysis of computational efficiency\n5. Dataset acquisition process is not explained", "relation_to_broader_scientific_literature": "The paper does not adequately situate itself within the broader literature. Key related works in neural operators, traditional numerical methods, and other symbolic regression approaches are not sufficiently discussed or compared.", "essential_references_not_discussed": "Several important references are missing, including:\n\n1. Works on neural operators (e.g., DeepONet [1], FNO [2]) (The algorithm presented in this study necessitates a dataset for both training and testing purposes. The nature of this task aligns with the fundamental concepts explored in neural operator research.)\n2. Traditional numerical methods for PDEs [3]\n3. Other reinforcement learning approaches for symbolic regression [4-9]\n\n\n[1] Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators\n[2] Fourier neural operator for parametric partial differential equations\n[3] Spline approximation, part 1: Basic methodology\n[4] Finite Expression Method for Solving High-Dimensional Partial Differential Equations\n[5] Deep Learning and Symbolic Regression for Discovering Parametric Equations\n[6] Deep Reinforcement Learning-Based Symbolic Regression for PDE Discovery Using Spatio-Temporal Rewards\n[7] Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE)\n[8] Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws.\n[9] Symbolic regression via neural-guided genetic programming population seeding", "other_strengths_and_weaknesses": "Strengths:\n\n1. Novel combination of reinforcement learning and symbolic regression\n2. Clear motivation and problem formulation\n3. Introduction of risk-seeking constant optimization and recursive exploration techniques\n\nWeaknesses:\n\n1. Limited experimental validation (only linear PDEs and ODEs)\n2. Lack of comparison with existing methods\n3. No analysis of computational efficiency\n4. Mathematical notation inconsistencies and errors\n5. Limited discussion of practical applications\n6. Missing supplementary material and source code", "other_comments_or_suggestions": "1. Address mathematical notation inconsistencies and errors (e.g., extra multiplication sign on line 153, potential issue with F definition on line 153, missing 's' in \"ordinary differential equation\" on line 156, inconsistent use of subscripts for x on line 131)\n2. Expand experimental validation to include nonlinear PDEs (e.g., NS equations, Gross-Pitaevskii equations)\n3. Compare with traditional numerical methods and other machine learning approaches\n4. Analyze computational efficiency and compare with traditional methods (e.g., spline interpolation)\n5. Provide source code and implementation details\n6. Clarify dataset acquisition process\n7. Discuss the method's limitations and practical applications more thoroughly\n8. Provide GPU details in the appendix (line 297)", "questions_for_authors": "1. How does SSDE compare to other reinforcement learning approaches for symbolic regression [2-5]?\n2. Can SSDE handle nonlinear PDEs like NS, Gross-Pitaevskii equations?\n3. What is the computational efficiency compared to traditional methods [1]?\n4. How was the training dataset acquired?\n5. What are the limitations of the proposed approach?\n6. How does the method scale with problem size and desired accuracy?\n7. Can the algorithm be parallelized on CPUs and GPUs?\n\n[1] Spline approximation, part 1: Basic methodology\n[2] Finite Expression Method for Solving High-Dimensional Partial Differential Equations\n[3] Deep Learning and Symbolic Regression for Discovering Parametric Equations\n[4] Deep Reinforcement Learning-Based Symbolic Regression for PDE Discovery Using Spatio-Temporal Rewards\n[5] Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE)\n[6] Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws.\n[7] Symbolic regression via neural-guided genetic programming population seeding", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {"ethical_review_concerns": "none"}}, {"review_id": "rXlP5pw8Jl", "reviewer": "Reviewer_vQDg", "summary": "The paper proposes a deep learning approach to obtain closed-form solutions for PDEs. The authors exploit this task through a Markov decision process and introduce an RL-based methodology. They also address acceleration and multi-dimensional problems, presenting an ablation study to support the proposed methods. While the proposed approach presents promising results, its applicability to general PDEs remains somewhat uncertain.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the feedback and suggestions, we will add clarification where needed and include suggestions as space permits.  Extra results: [here](https://anonymous.4open.science/r/SSDE-A47B/SSDE_icml2025_rebuttal.pdf). All figures/tables are from this link.\n\n**Q1: Advantages of Formulating PDE Solutions as an RL Problem.**\n\nThe RL framework is critical for symbolic expression discovery because:\n\n* Symbolic Expressions Lack Gradient Backpropagation: Unlike PINN’s numerical solutions (which rely on gradient backpropagation), symbolic expressions require gradient-free optimization. Reinforcement learning (RL) provides policy gradients guided by reward signals, enabling direct exploration of interpretable expressions.\n* Superior Trade-off Between Complexity and Interpretability: Genetic programming (GP) methods optimize solutions via fitness functions but often overfit by increasing expression length and adding excessive constants. In contrast, our RL framework balances accuracy and simplicity by directly rewarding concise, interpretable expressions (e.g., avoiding unnecessary terms).\n\nRegarding RSCO’s impact:\n\nFigure 1 shows that RSCO reduces solution time despite a slight increase in explored expressions. Interestingly, under identical seeds, RSCO often finds solutions faster (e.g., mirroring hPINN’s boundary-condition prioritization for convergence acceleration).\n\n**Q2: Extending to Weak/Distributional Solutions**\n\nYou raise an important point about generalizing to weak/distributional solutions. Potential extensions include:\n\n1. Enriching the Symbol Library: Introducing neural operators or special functions (e.g., Bessel functions) to capture non-smooth solutions.\n2. Combining Multiple Solutions: Treating SSDE’s outputs as particular solutions and aggregating results across multiple runs (analogous to classical DE solving techniques).\n\nWe acknowledge this as a promising future direction and will explore it further.\n\n**Additional Clarifications**\n\n*  While we are unable to conduct a full run-time and memory-usage comparison due to time constraints, we can provide quantitative context for SSDE’s run-time: For the 3D heat equation, SSDE required ~2,700 seconds to converge. RSCO’s design ensures that memory usage remains comparable to DSR , as it avoids storing large amount of computation graphs for automatic differentiation. \n* We will thoroughly check and fix grammatical and clerical errors in the final submission. \n* The code and supplementary experiments are provided to clarify implementation details.\n\nWe appreciate your feedback and welcome further discussion to address any remaining concerns.", "claims_and_evidence": "The advantages of treating the closed-form solution of the PDE as an RL problem should be articulated more clearly.\nBased on the experimental results, it seems that the proposed approach demonstrates superior performance compared to existing methods, even in the absence of RSCO or recursive exploration. I would appreciate it if you could further elaborate on why you believe the proposed framework is particularly advantageous for handling symbolic expressions. Additionally, since RSCO was introduced with the intention of accelerating convergence, why the lack of RSCO results in a decrease in performance, as shown in the results presented in Figure 2?", "methods_and_evaluation_criteria": "The paper currently focuses on smooth solutions, but most general PDEs do not possess classical solutions. As a result, generalized concepts of solutions, such as weak or distributional solutions, are considered, which are generally not unique. For the proposed approach to be of broader significance, it could be applicable/extendable to such scenarios. Could you elaborate on how the methodology might be extended to handle such cases? I am not requesting the addition of experiments but rather would appreciate a discussion on the potential for extending the proposed methodology to such practical scenarios.", "theoretical_claims": "The paper does not include any theoretical content.", "experimental_designs_or_analyses": "The proposed methodology demonstrates superior performance compared to existing approaches. However, the incorporation of policy gradients, RSCO, and recursive exploration likely increases its complexity. To assess this, I recommend comparing computational time and memory usage between the proposed methodology and baseline methods across 1D, 2D, and 3D cases.", "relation_to_broader_scientific_literature": "If the closed form of the solution can be obtained, it would not only enhance the interpretability of solutions predicted by deep learning models but also provide valuable insights into the mathematical understanding of PDEs. This task effectively leverages the strengths of deep learning.", "essential_references_not_discussed": "No essential related works are missing.", "other_strengths_and_weaknesses": "Please refer to the comments provided above.", "other_comments_or_suggestions": "There are some minor typos, for example:\n\n* Line 153, first column: '$R\\times\\times\\cdots\\times R^n$' should be '$R\\times R^n\\times\\cdots\\times $'.\n\n* Line 204, second column: 'We also note that the skeleton satisfy deterministic' should be 'We also note that the skeleton satisfies deterministic.'", "questions_for_authors": "Please refer to the comments provided above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "w32fYwDsCM", "reviewer": "Reviewer_WYoN", "summary": "This work presents an approach to discover analytical solutions of PDEs using reinforcement learning. The authors combine several concepts in this work, including an optimization method constrained by IC/BCs, an iterative construction approach which generates solution skeletons and subsequently refines parameters of this expression, and finally an approach to find solutions along each dimension such that the approach can be extended to PDEs in multiple spatial dimensions. The general algorithm, which is visualized nicely in Figure 1, generates an expression from an RNN. Points from the domain are subsequently sampled and loss with respect to the ICs, BCs, and internal domain is calculated. Finally, a reward function is used to update the RNN in a reinforcement learning approach. \n\nThe authors present several experiments on ODEs, elliptic, and parabolic PDEs. Compared to the baselines, the proposed approach shows superior performance in most cases.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the feedback and suggestions, we will add clarification where needed and include suggestions as space permits.  Extra results: [here](https://anonymous.4open.science/r/SSDE-A47B/SSDE_icml2025_rebuttal.pdf). All figures/tables are from this link.\n\n**Q1: Add the experiments to support the efficiency improvement brought by RSCO.**\n\nAdded ablation experiments of RSCO are perfomed on the Gamma benchmarks with single variable as shown in Fig 1. We conducted 20 experiments on each benchmark and counted the average time required to successfully find an analytical solution and the average number of explored expressions. Despite the RSCO method requires a slight increase in the number of exploration expressions, it greatly improves efficiency. \n\n**Q2: \"Novel closed-form solutions\".**\n\nThis refers to SSDE’s parametric expression framework for high-dimensional PDEs (Appendix Table 10). RecurExp recursively constructs solutions, drastically shrinking search space and boosting efficiency.\n\n**Q3: Can SSDE discover solutions for more complex equations? Lacks of limitations section.**\n\nAdded experiments on the nonlinear wave equation benchmarks are shown in Tabel 1 and Tabel 2. SSDE succeeds on nonlinear Wave2D/Wave3D benchmarks with non-polynomial solutions.\n\nWe also acknowledge the limitations of SSDE in solving systems of DEs due to their inherent coupled-solution nature. Exploring the full solution space requires multiple RNNs, and evaluating reward-guided gradients for each solution necessitates novel methods. We aim to address this challenging issue in future work.\n\n\n**Q4: Discussion about baselines and evaluation method.**\n\nThank you for your suggestions on baseline methods and evaluation metrics. First, we have added relative L2 error comparisons for baselines (Table 3). Regarding baselines:\n\n* We did not include PINN or DSR individually as baselines because our method focuses on discovering analytical solutions, which differs fundamentally from numerical approximation (PINN) or symbolic regression (DSR).\n* PINN+DSR represents a natural approach for symbolic solution discovery (solving DEs numerically with PINN followed by symbolic fitting via DSR), hence its inclusion as a baseline. For reference, we report PINN’s convergence error relative to analytical solutions (Table 4) and DSR’s performance on ground-truth data (Table 5).\n* We use the KAN method as a baseline because they did do the work of solving the symbolic solution of the differential equation in the paper, but this task may indeed be beyond the capabilities of KAN. \n* We added Mathematica 14.1 and the Finite Expression Method (FEX) [1] as stronger baselines (Table 3). While KAN attempts symbolic solution discovery via spline-based activation functions, FEX directly uses symbolic neurons. Both still fail to recover exact closed-form solutions for our benchmarks.\n\nWe performed parameter tuning for the PRGP baseline, which improved its performance across benchmarks. Although the GP method can generate highly complex expressions to achieve good performance on the Van der Pol benchmark, it fails to recover closed-form solutions in high-dimensional settings. You mentioned that Poisson's equation could be solved in their paper, but the sample they tested in their original paper was very simple. Although the samples we tested look like  polynomials, the symbolic expressions involved are very complex and DSR cannot effectively identify them too. Our approach still works very well on these high-dimensional problems.\nThe performance of the PRGP method can also refer to [2].\n\n[1] Finite Expression Method for Solving High-Dimensional Partial Differential Equations\n\n[2] An interpretable approach to the solutions of high-dimensional partial differential equations.\n\nAnswer for `Question For Authors`:\n\nQ1: Symbol libraries are predefined, as in RL/GP-based heuristic search algorithms.\n\nQ2: SSDE requires retraining per PDE system. But transferability is an interesting concern, the solution of different DEs may have very large differences.\n\nQ3: Yes, MSE difference stems from truncating constants to significant digits.\n\nQ4: SSDE’s sensitivity to IC/BC terms mirrors Mathematica’s behavior.\n\nQ5: As we replied to the reviewer ateq(Q3), symbol library size increases search complexity. `RecurExp` it is used to solve the problem of increased search space caused by increased variables.\n\nWe will thoroughly check and fix grammatical and clerical errors in the final submission. Thank you for your constructive feedback and for indicating your openness to revising your evaluation. We have thoroughly addressed all your questions and concerns through additional experiments, which are now included in the rebuttal document here and highlighted in Tables 2/3/4/5,Figure1. These results directly address your specific requests and strengthen the robustness of our claims. We kindly ask that you reassess our submission in light of these additions.", "claims_and_evidence": "The authors present an ablation which supports their claim that the RSCO policy improves recovery rate and performance. As shown in Figure 2, the recovery rate is greatly diminished when RSCO is not used, and even moreso when recursive exploration is not used. The authors also claim that this approach will result in great improvements in accuracy and efficiency (4x speed up), but I do not see any experiments which support this claim. I would kindly request the authors to clarify this, or add the experimental results which support this claim.\n\nThe authors also present several tables with the closed form expressions which closely match the true solution; however, the claim that these are *novel* closed-form solutions is questionable. There are many approaches which have already been able to show the ability to learn closed-form solutions of ODEs and PDEs, including \"Solving differential equations with genetic programming\" (Tsoulos et al. 2006), and other approaches in Genetic Programming, Ant Colony Programming, etc. These approaches have even been shown to work for high dimensional (2 and 3 dimensional) PDEs. \n\nFinally, the authors claim the ability to efficiently discover solutions for *complex* equations. However, I would argue that the DEs considered in this work are relatively simple. More challenging baselines, such as compressible Euler or Shallow Water Equations would be necessary to explore the performance of SSDE for more challenging equations. Likewise, the solutions of all DEs are polynomial. I see in the solution to the Heat2D (table 3), that SSDE identifies a solution which takes the log of an exponential. These effectively cancel out, but suggest the ability to find non-polynomial solutions. Solutions with sin and cos terms are also shown in the Gamma dataset.  Nonetheless, this suggests the approach may have challenges in finding more complex solutions. I believe this claim merits additional experiments to be supported.", "methods_and_evaluation_criteria": "The differential equations tested in this work make sense, but as mentioned in my previous comment, more challenging baselines should also be explored.\n\nThe evaluation of the baselines methods seems questionable. In their work, PR-GSPR is able to recover the analytical solutions for Poisson's equation, which is also tested in this work. Nonetheless, the results in Table 2 show that it failed to recover the correct solution across the experiments. Furthermore, the appendix states that hyperparameter exploration/tuning was not used for this approach, while a fairly large sweep was used for SSDE. A sweep should also be performed for the baselines, where applicable. \n\nPINN+DSR also seems to be a weak benchmark, as one first solves the DEs with PINNs, and then fits an analytical solution to this, both of which will introduce some error. Of course, PINNs can be extremely sensitive to initialization and difficult to optimize, also introducing non-physical dynamics which might be very difficult to model with DSR. I believe simply providing the error of the PINNs with respect to the analytical solution would serve to better illustrate the performance of the method. As this work seems very similar to DSR, I believe this should still be considered as a baseline, however, it should be applied to data from the true solution. \n\nThe KAN baseline is interesting, but I believe this may be quite a bit outside of the realm in which KANs are usually used; however, I will admit I am not as familiar with this work. \n\nI believe established analytical solvers, such as Eureqa and Mathetmatica would serve as stronger baselines in these experiments.", "theoretical_claims": "Not applicable, as proofs are not included.", "experimental_designs_or_analyses": "The design of the proposed experiments is sound, but more challenging equations and stronger baselines must be considered, as previously mentioned. \n\nIn terms of evaluation, I would suggest the use of a relative L1 or L2 error, in addition to the physics loss. This can be provided in the appendix, but it would be helpful to understand the performance of the approaches.", "relation_to_broader_scientific_literature": "The key contribution to the literature is the extension of deep symbolic regression (DSR) to PDEs. The approach outlined in this work is fundamentally similar to DSR, using RNNs to generate expression trees, whose terms are taken from a defined library. Likewise, both approaches use reinforcement learning with additional risk seeking behavior to learn how an optimal policy for generating symbolic expressions. \n\nMore broadly, the paper focuses strongly on the performance improvements which come from risk-seeking optimization. Work in this field has recently been explored, referenced in the paper (see also arxiv 2302.09339). This work helps to illustrate the usefulness of risk-seeking optimization in RL.", "essential_references_not_discussed": "The authors discuss genetic programming, providing a reference to work by Oh et al., but the foundational work in this field for DEs was introduced in 2006, \"Solving differential equations with genetic programming,\" by Tsoulos and Lagaris.", "other_strengths_and_weaknesses": "Strengths\n1. the authors design closed-loop extension of DSR which is able to be applied effectively to DEs\n2. several elements of the proposed architecture improve the efficiency of the approach\n3. the approach can be extended to arbitrarily-high-dimensional PDEs\n\nWeaknesses\n1. the paper lacks a limitations section, which is critical for clarity and transparency on the work\n2. more difficult baseline experiments should be used to understand how well this approach can generalize to PDEs.\n3. baseline models are weak, and not fairly optimized for the test cases\n4. claims regarding efficiency are not supported experimentally: run times and number iterations are not provided", "other_comments_or_suggestions": "There are minor grammatical mistakes and typos throughout the work. Fixing these would improve reading. \n\nLine 356, right column, a reference to comparison between SSDE and the numerical solution is said to exist in Appendix B, but this section only contains computing and hyperparameter details.\n\nI believe the relationship of this work to the DSR paper should be discussed in more detail. There is some discussion in the appendix, but this is limited.", "questions_for_authors": "1. Does the symbol library need to be defined a priori? For example, in the Gamma dataset, there are several solutions which contain sin and cos. Is it necessary that this be prescribed before finding the solution to the PDE?\n\n2. Can knowledge be transferred between systems of PDEs, or is it necessary to train the model from scratch for each set of PDEs?\n\n3. Table 11 in the appendix shows that the PINN+DSR and SSDE both obtain the correct analytical solution for $\\Gamma_1$, yet SSDE has an MSE 22 orders of magnitude better than PINN+DSR. Is this due to small numerical errors in the constants of the solution which aren't shown in the provided solution, or is there some other reason?\n\n4. What is the sensitivity to the IC and BC loss terms? Is the approach sensitive to hyperparameter tuning here? \n\n5. The paper mentions exponential growth in the search space for an increase in the number of symbols. Does this limit the application of this approach to find solutions which have many terms? \n\nI would like to thank the authors for their time, and I would gladly raise my evaluation if additional experimental results are provided to address my questions and concerns.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This submission presents a reinforcement learning based approach for solving differential equations symbolically. The proposed method is quite novel and is used in experiments on ODEs, elliptic and parabolic PDEs, and eventually also the 3d wave equation. The improvement over the baselines is clear, although there was some disagreement among reviewers about the sufficiency of the experiments considered. In fact, 2/3 reviewers urged the authors to consider more realistic, challenging PDEs such as those coming from physics; the authors have agreed and presented some more experiments in this direction. The authors addressed the majority of these concerns in the rebuttal period and it is my opinion that the paper should be accepted despite the fact that one reviewer recommends rejection.
Decision: Accept (poster)