FoundPDE: A generative foundation model for data-driven symbolic discovery of partial differential equations

Langchen Liu; Junlin Luo; Min Zhu; Zhuoran Yang; Lu Lu

physics.comp-phOct 18, 2026

FoundPDE: A generative foundation model for data-driven symbolic discovery of partial differential equations

Authors: Langchen Liu, Junlin Luo, Min Zhu, Zhuoran Yang, Lu Lu

Organizations: Department of Statistics and Data Science, Yale University, New Haven, CT 06511, USA · Department of Mechanical Engineering, Stanford University, Stanford, CA 94305, USA · Department of Chemical and Environmental Engineering, Yale University, New Haven, CT 06511, USA

Abstract

Discovering governing partial differential equations (PDEs) from observed data is a central inverse problem in computational science and engineering: the goal is not only to predict trajectories, but to identify interpretable symbolic laws that can be analyzed and reused. Existing PDE discovery methods often depend on hand-designed libraries, accurate numerical derivatives, or per-equation symbolic search, making them fragile under sparse observations, noisy data, broad equation families, and unseen symbolic compositions. Here, we introduce FoundPDE, a pretrained data-to-symbol foundation model for PDE discovery. FoundPDE maps discretized spatio-temporal trajectories to symbolic equations with a transformer, uses Fourier and first-order spatial-derivative augmentation to expose complementary trajectory information, refines numerical coefficients, and supports downstream extension with improved generalization through few-shot fine-tuning and hierarchical discovery. On a clean in-domain benchmark spanning 20 PDE families, our model reaches 97.50% type recognition. Under coefficient shifts beyond the training support, type recognition remains 90.53%; under multiplicative Gaussian noise on the model input, it is 94.43% at a 2% noise level but falls to 52.09% at 10%. For unseen composite equations, few-shot fine-tuning reaches over 90% exact recovery on four held-out families, and we further develop a hierarchical discovery algorithm that recovers two-term composites without fine-tuning FoundPDE. By amortizing symbolic discovery across equation families, FoundPDE shifts PDE discovery from per-dataset search toward a reusable pretrained model that extends to new physics with limited additional data, a step toward foundation models that uncover interpretable governing laws directly from scientific data.

Figures & tables

Figure 1: FoundPDE performs end-to-end symbolic PDE discovery from discretized solution data. Stage A constructs the symbolic training corpus and tokenizes each governing equation. Stage B solves the sampled equations under randomized initial conditions to produce spatio-temporal observations. Stage C optionally augments each trajectory with first-order spatial-derivative and Fourier channels. Stage D encodes these observations and autoregressively decodes a symbolic PDE in Polish notation. Stage E returns the predicted expression. The optional BFGS branch updates only its floating-point coefficients, whereas hierarchical discovery reuses the pretrained model to recover selected composite operators sequentially.

Figure 2: The hierarchical discovery algorithm sequentially recovers a composite PDE from familiar operators. The first FoundPDE pass identifies the dominant operator $u_{t}=A(u)$ from a trajectory of the composite equation $u_{t}=A(u)+B(u)$ . The residual-computation panel depicts one update at time step $n$ : the neural operator supplies $\hat{u}^{\,n+1}$ , the dominant-equation solve supplies $u_{A}^{\,n+1}$ , and their difference updates $\tilde{u}^{\,n+1}$ . Repeating this calculation over all time steps constructs the residual dataset $\tilde{\mathcal{D}}$ , as formalized in Algorithm 1 . A second FoundPDE pass identifies $B(u)$ , after which the two operators are combined and calibrated with BFGS.

Figure 3: Effects of feature augmentation across symbolic and numerical metrics. All four variants are trained on the same 512,000 samples for 100 epochs. ( A ) Projected training time for the four model variants, from the measured time per optimization step (Appendix B ). ( B ) The augmented models improve type recognition rate (TRR) relative to the barebone baseline. The Fourier-and-derivative model gives the highest structural accuracy. ( C ) The derivative model gives the lowest mean relative coefficient squared error (MRCSE). The evaluable MRCSE counts are 4,929, 4,959, 4,980, and 4,992 for the barebone, Fourier, derivative, and Fourier-and-derivative models, respectively. ( D ) The barebone model gives the lowest mean relative

L^{2}

error; the evaluable count is 5,115 for every variant. ( E ) The training text-loss curves indicate that the final training loss is not by itself a reliable proxy for downstream symbolic recovery.

Figure 4: In-domain recovery varies across PDE families and evaluation criteria. ( A ) Per-type TRR. Nine families, including heat, advection, porous medium, KdV, wave, Sine–Gordon, and Fokker–Planck, are recognized for every test sample, whereas the diffusion-reaction families have the lowest TRR. ( B ) Mean relative

L^{2}

error by family. The family ordering differs from that in panel A, showing that structural recognition and relative

L^{2}

error need not rank PDE families identically. ( C ) Global relative

L^{2}

error distribution over 5,115 reconstructable predictions. ( D ) Global MRCSE distribution over 4,992 type-matched predictions with extractable coefficients; its long tail shows that a subset of evaluable predictions remains difficult to calibrate.

Panel of Figure 5 — Figure 5: Representative in-domain examples E1–E5 show accurate symbolic recovery on transport, diffusion, and higher-order families. This group covers Burgers, heat, advection, Cahn–Hilliard, and conservation with cubic flux. All five predictions have the correct symbolic type, and the reconstructed fields closely match the ground truth. E4 recovers the Cahn–Hilliard structure with a visibly different fourth-order coefficient while retaining a small relative $L^{2}$ error.

Panel of Figure 6 — Figure 6: Representative in-domain examples E6–E10 separate structural recovery from numerical calibration. This set includes linear-flux conservation, sine-flux conservation, logistic reaction-diffusion, square-logistic reaction-diffusion, and KdV. All five predictions have the correct symbolic type. E7–E10 agree well with the ground truth in field space, whereas small errors in the transport and diffusion coefficients of E6 produce a visible phase shift.

Panel of Figure 7 — Figure 7: Representative in-domain examples E11–E15 show that FoundPDE handles porous-medium, hyperbolic, and driven-diffusion equations well. This group contains porous medium, wave, Klein–Gordon, Sine–Gordon, and Fokker–Planck examples. The porous-medium, wave, Sine–Gordon, and Fokker–Planck cases show close visual agreement between predicted and ground-truth fields, while Klein–Gordon illustrates a regime where the symbolic type is correct but coefficient sensitivity is still visible.

Figure 8: Corruption of the augmented input affects TRR, MRCSE, and relative

L^{2}

error at different rates. The augmented tensor contains solution, first-order spatial-derivative, Fourier-real, and Fourier-imaginary channels, and

\eta

is the relative standard deviation of the multiplicative Gaussian noise applied to each entry. ( A ) TRR as

\eta

increases from clean input to

0.1

. ( B ) Relative

L^{2}

error after reconstructing each decoded PDE from the clean initial condition and comparing it with the noise-free trajectory. ( C ) Mean and median MRCSE on a logarithmic scale; their growing separation indicates an increasingly heavy-tailed distribution. ( D ) Prediction categories over the fixed 5,120-example denominator, showing that wrong-PDE-type and no-matched-PDE outcomes become more frequent as corruption increases.

Figure 9: Sensitivity to augmented-input corruption depends on the PDE family. ( A ) Per-type TRR as $\eta$ increases. Fokker–Planck, Sine–Gordon, inviscid Burgers, and KdV retain high TRR at $\eta=0.1$ , whereas heat, wave, advection, Cahn–Hilliard, and several diffusion-reaction families show larger decreases. ( B ) Most frequent incorrect type transitions at $\eta=0.1$ . The $x$ -axis reports $100\,N(y\to\hat{y})/N(y)$ , where $N(y\to\hat{y})$ is the number of samples whose true PDE type $y$ is predicted as $\hat{y}$ , and $N(y)$ is the total number of samples with true type $y$ . The dominant transitions occur among related reaction-diffusion families and from heat and wave toward reaction or nonlinear-source families.

Figure 10: A corrupted heat example shows a spurious reaction term changing both the symbolic type and the reconstructed trajectory. The top row visualizes the solution channel presented to the model at each corruption level. The bottom row reconstructs every decoded PDE from the same clean initial condition, and relative $L^{2}$ is computed against the noise-free trajectory. ( A ) The clean input is recovered as a heat equation with relative $L^{2}=0.001$ . ( B ) At $\eta=0.02$ , the prediction remains a heat equation and gives relative $L^{2}=0.002$ . ( C ) At $\eta=0.05$ , the prediction adds a logistic reaction term, and the reconstruction has relative $L^{2}=0.146$ . ( D ) At $\eta=0.1$ , the logistic term grows and the relative $L^{2}$ error reaches $0.661$ .

Ground truth	FoundPDE	SINDy	Weak SINDy	DISCOVER
$u_{t}=0.000802u_{xx}+0.0472u$	$u_{t}=0.000722u_{xx}+{\color[rgb]{0,0.55,0}0.0492u}$ (19/20)	$u_{t}=0.000799u_{xx}+{\color[rgb]{0,0.55,0}0.0470u}$ (16/20)	$u_{t}=0.000792u_{xx}+{\color[rgb]{0,0.55,0}0.0472u}$ (20/20)	$u_{t}=0.0008u_{xx}+{\color[rgb]{0,0.55,0}0.0486u}$
$u_{t}=0.00228u_{xx}+0.794u(1-u)$	$u_{t}=0.00231u_{xx}+{\color[rgb]{0,0.55,0}0.813u(1-u)}$ (20/20)	$u_{t}=0.00228u_{xx}+{\color[rgb]{1,0,0}0.778u-0.729u^{2}}$ $\quad{\color[rgb]{1,0,0}{}-0.084u^{3}+0.034u^{4}}$ (6/20)	$u_{t}={\color[rgb]{1,0,0}0.962u-1.024u^{2}}$ (3/20)	$u_{t}=0.0034u_{xx}+{\color[rgb]{0,0.55,0}0.0309u(1-u)}$
$u_{t}=0.00496u_{xx}+0.239u^{2}(1-u)$	$u_{t}=0.00447u_{xx}+{\color[rgb]{1,0,0}0.529u^{2}(1-u)^{2}}$ (13/20)	$u_{t}=0.00492u_{xx}+{\color[rgb]{0,0.55,0}0.248u^{2}-0.253u^{3}}$ (16/20)	$u_{t}=0.00489u_{xx}+{\color[rgb]{0,0.55,0}0.241u^{2}-0.242u^{3}}$ (20/20)	$u_{t}=0.001u_{xx}+{\color[rgb]{1,0,0}0.0428u+0.005}$
$u_{t}=0.00344u_{xx}+0.125u^{2}(1-u)^{2}$	$u_{t}=0.00325u_{xx}+{\color[rgb]{0,0.55,0}0.171u^{2}(1-u)^{2}}$ (16/20)	$u_{t}=0.00345u_{xx}+{\color[rgb]{0,0.55,0}0.125u^{2}}$ $\quad{\color[rgb]{0,0.55,0}{}-0.249u^{3}+0.125u^{4}}$ (12/20)	$u_{t}=0.00343u_{xx}+{\color[rgb]{0,0.55,0}0.126u^{2}}$ $\quad{\color[rgb]{0,0.55,0}{}-0.252u^{3}+0.126u^{4}}$ (10/20)	$u_{t}=0.0034u_{xx}+{\color[rgb]{1,0,0}0.0333u(1-u)}$
Correct over 80 trajectories	68	50	53	–

Table 1: Sparse regression is competitive on clean diffusion-reaction data when its library contains the true terms. All test PDEs have the form

u_{t}=\nu u_{xx}+\rho R(u)

and are observed on the

32\times 64

grid. Each entry is the equation recovered from the first initial condition, the instance also used for DISCOVER; parentheses give the number of the 20 initial conditions of the same equation for which the recovered terms match the ground truth. FoundPDE entries are raw decoded equations before BFGS refinement. SINDy and weak SINDy use the library

\{1,u,u^{2},u^{3},u^{4},u_{x},u_{xx}\}

, and their polynomial reaction terms are shown in expanded form. Green entries match the ground-truth structure, while red entries indicate a structural mismatch. DISCOVER was run on the first initial condition only. Appendix Table 12 reports the same comparison under noise.

Ground truth	DISCOVER $32\times 64$	DISCOVER $128\times 128$	FoundPDE $32\times 64$
Burgers $u_{t}=0.00458u_{xx}-0.697uu_{x}$	$u_{t}=0.0023u_{xx}-0.3596(u^{2})_{x}$	$u_{t}=0.0056u_{xx}-0.8332(u^{2})_{x}$	$u_{t}=0.00357u_{xx}-0.533uu_{x}$
Heat $u_{t}=0.00409u_{xx}$	$u_{t}=0.0164u_{xx}$ $/((x^{2})_{xx})^{3}$	$u_{t}=0.0053u_{xx}/(x^{2})_{xx}$	$u_{t}=0.0041u_{xx}$
Advection $u_{t}=-0.562u_{x}$	$u_{t}=-0.092(x^{3})_{xxx}u_{x}$	$u_{t}=-0.0594(x^{3})_{xxx}u_{x}$	$u_{t}=-0.775u_{x}$
Cahn–Hilliard $u_{t}=-2.3875{\times}10^{-6}u_{xxxx}$ $-1.5u_{x}u_{x}-1.5uu_{xx}$	$u_{t}=0.0089u_{xx}-0.0077u(u_{xx}+u_{xx})$ $-2.9985u^{2}+1.7466u$	$u_{t}=-0.0056uu_{xx}+0.0028u_{xx}+0.9546u$ $-0.8303u(u+u)$	$u_{t}=-2.39{\times}10^{-6}u_{xxxx}$ $-1.5u_{x}u_{x}-1.5uu_{xx}$
Conservation (cubic flux) $u_{t}=-0.341u^{2}u_{x}+0.00057u_{xx}$	$u_{t}=0.0009(u^{2})_{xx}$ $-0.1176(uuu)_{x}$	$u_{t}=-0.454(u^{3})_{x}$ $+0.003(u^{3})_{xx}$	$u_{t}=-1.08u^{2}u_{x}+0.000544u_{xx}$
Conservation (linear flux) $u_{t}=-0.98u_{x}+0.00267u_{xx}$	$u_{t}=-0.165u_{x}(x^{3})_{xxx}$	$u_{t}=-0.844u_{x}$	$u_{t}=-1.12u_{x}+0.00257u_{xx}$

Table 2: A family-level representative comparison shows more consistent symbolic-family recovery by FoundPDE than by per-instance DISCOVER. The ground-truth column gives the target of the sparse

32\times 64

DISCOVER run. The

128\times 128

DISCOVER entry is a separately instantiated representative of the same PDE family, and the FoundPDE entry is selected from the no-noise

32\times 64

evaluation set using the same family and the nearest available coefficient setting. The columns therefore compare recovery behavior across shared PDE families, not predictions on identical coefficients and initial conditions, and coefficient values should not be compared directly across columns. DISCOVER is optimized independently for every displayed instance, whereas all FoundPDE entries are produced by the same pretrained model without per-instance optimization. Green indicates recovery of the target symbolic family for the corresponding instance, orange indicates an algebraically related or simplifiable expression with a non-canonical or incomplete form, and red indicates a structural mismatch.

Figure 11: Strict parameter-OOD evaluation reveals a moderate aggregate decrease and concentrated family-level failures. ( A ) Aggregate TRR and mean relative $L^{2}$ error over 5,120 samples whose coefficient supports are disjoint from the in-domain ranges. ( B ) The ten most frequent incorrect type transitions. ( C ) Per-family TRR and normalized relative $L^{2}$ score, computed as $100\,\mathrm{clip}(1-\mathrm{Rel}L^{2}/0.60,0,1)$ for visualization. ( D ) Per-family TRR versus mean relative $L^{2}$ error. Advection and viscous linear conservation illustrate that high structural recognition does not necessarily imply low relative $L^{2}$ error.

Family	Samples	TRR (%) $\uparrow$	Mean Rel. $L^{2}$ $\downarrow$	Median MRCSE $\downarrow$	Most frequent incorrect type
Burgers	272	100.00	0.088	0.008	None observed
Inviscid Burgers	272	100.00	0.122	0.005	None observed
Advection	272	99.63	0.676	2.568	None observed
Diffusion–reaction, logistic	272	99.63	0.021	0.007	None observed
Viscous conservation, linear flux	272	99.63	0.326	0.011	Sine-flux conservation
Cahn–Hilliard	264	99.62	0.000	0.003	None observed

Table 3: Strict parameter-OOD performance is strong for most pretrained operator families but degrades for sine-flux conservation and for large coefficient shifts in heat, Klein–Gordon, and wave. TRR includes empty decoder outputs as failures. Mean relative

L^{2}

error is computed on 5,011 reconstructable predictions. Median MRCSE is computed on 4,630 type-matched predictions with extractable packed coefficients. Both numerical metrics use the raw decoded equations before BFGS refinement. The final column reports the most frequent parsed incorrect type; “none observed” indicates that the only failures, if any, are empty decoder outputs.

Figure 12: Single-family checkpoint adaptation reaches high composite recovery within 20 epochs, while joint four-family scratch training requires a longer horizon. ( A ) Separate pretrained fine-tuning runs (solid curves) versus one jointly trained scratch model (dashed curves). Both protocols use 10 optimizer updates per epoch, but their batch composition differs. By epoch 20, fine-tuning reaches

91.0\%

exact recovery on advection–Burgers,

94.0\%

on advection–cubic-Burgers,

95.0\%

on advection–sine, and

97.0\%

on Burgers–cubic-Burgers. At epoch 100, scratch training remains below these values on three families but reaches

100\%

on Burgers–cubic-Burgers. The curves are single-run trajectories, and epochs 20 and 100 are fixed comparison horizons rather than test-selected checkpoints. ( B ) Retention along the advection–Burgers fine-tuning run: target-family TRR increases rapidly, while the fixed-denominator macro-average TRR across the 20 original families decreases from

96.52\%

88.50\%

. Missing outputs are counted as failures.

Ground truth	Decomposition results
$u_{t}=-1.08\,uu_{x}+0.00541\,u_{xx}$	Step 1 (Dominant):	$u_{t}=-1.14\,uu_{x}$
	Step 2 (Subordinate):	$u_{t}=0.00476\,u_{xx}$
	Combined:	$u_{t}=-1.14\,uu_{x}+0.00476\,u_{xx}$
	After BFGS:	$u_{t}=-1.080\,uu_{x}+0.005407\,u_{xx}$
$u_{t}=-0.902\,u^{2}u_{x}+0.00334\,u_{xx}$	Step 1 (Dominant):	$u_{t}=-0.989\,u^{2}u_{x}$
	Step 2 (Subordinate):	$u_{t}=0.00505\,u_{xx}$

Table 4: The hierarchical discovery algorithm recovers the tested unseen two-term composites through sequential operator identification and BFGS calibration.

Appendix figures & tables9 assets

Supplementary material from the paper’s appendix.

Appendix

Family	Canonical form	Reference setting	In-domain support
Heat / Diffusion	$u_{t}=cu_{xx}$	$c_{0}=3\times 10^{-3}$ , $t_{f}=3$	$c/c_{0}\in[0.1,1.9]$
Porous medium	$u_{t}=(u^{m})_{xx}$	$t_{f}=0.1$	$m\in\{2,3,4\}$
Klein–Gordon	$u_{tt}=c^{2}u_{xx}-m^{2}c^{4}u$	$c_{0}=1$ , $m_{0}=0.1$ , $t_{f}=1$	$c/c_{0},m/m_{0}\in[0.1,1.9]$
Sine–Gordon	$u_{tt}=u_{xx}-c\sin(u)$	$c_{0}=1$ , $t_{f}=1$	$c/c_{0}\in[0.1,1.9]$
Cahn–Hilliard	$u_{t}=-\epsilon^{2}u_{xxxx}-6(uu_{x})_{x}$	$\epsilon_{0}=10^{-2}$ , $t_{f}=0.5$	$\epsilon/\epsilon_{0}\in[0.1,1.9]$
KdV	$u_{t}=-\delta^{2}u_{xxx}-uu_{x}$	$\delta_{0}=0.022$ , $t_{f}=1$	$\delta/\delta_{0}\in[0.1,1.9]$

Appendix

Table 5: The synthetic corpus contains 20 PDE types constructed from the one-dimensional benchmark families of PROSE-PDE [ 13 ] . The table reports each family’s canonical form, reference setting, and in-domain parameter support. Diffusion–reaction and conservation-law families each contain multiple symbolic types, yielding 20 PDE types in total.

Figure 13: The principal model converges with respect to symbolic recovery before its training loss stops decreasing. ( A ) Per-epoch mean training loss of the principal model (512,000 samples, 100 epochs) and of the same model trained on 128,000 samples for 50 epochs. ( B , C ) Teacher-forced cross-entropy of the principal model at each saved checkpoint on 5,120 training samples and on the 5,120 test samples, for equation-structure tokens and for coefficient tokens (sign, mantissa, and exponent). ( D ) Test TRR of each checkpoint.

Setting	128K samples, 50 epochs	512K samples, 100 epochs
Clean in-domain	93.40	97.50
Corruption, $\eta=0.02$	89.86	94.43
Corruption, $\eta=0.1$	49.26	52.09
Strict parameter OOD	89.36	90.53
Heat	98.53	67.28
Klein–Gordon	85.23	71.97

Appendix

Table 6: A larger training budget improves recognition overall but changes the parameter-OOD failure profile. Both models use the Fourier-and-derivative input and the same architecture. Values are TRR (%) with the complete 5,120-sample denominators; the per-family rows are strict parameter-OOD results.

Composite family	Symbolic form
Advection–Burgers	$u_{t}=-\alpha u_{x}-kuu_{x}+\nu u_{xx}$
Advection–cubic Burgers	$u_{t}=-\alpha u_{x}-ku^{2}u_{x}+\nu u_{xx}$
Advection–sine	$u_{t}=-\alpha u_{x}-k\cos(u)u_{x}+\nu u_{xx}$
Burgers–cubic Burgers	$u_{t}=-\alpha uu_{x}-ku^{2}u_{x}+\nu u_{xx}$

Appendix

Table 7: Custom composite PDE families used for few-shot fine-tuning. These equations combine transport, nonlinear flux, and diffusion terms that are individually familiar from pretraining but not present as full composite symbolic forms. The initial conditions for these composite datasets use a sinusoidal superposition generator rather than the Gaussian-distribution or Gaussian-process initial conditions used by some PROSE-PDE families.

Equation type	Training initial condition	In-domain test initial condition
Heat, Diffusion–reaction, Klein–Gordon, Sine–Gordon, Cahn–Hilliard, Viscous conservation, Inviscid conservation, KdV	Sinusoidal superposition with $n_{\max}=2$	Independent draw from the same generator
Advection, wave	Gaussian distribution with $N=2$	Independent draw from the same generator
Fokker–Planck	Gaussian distribution with $N=1$	Independent draw from the same generator
Porous medium	Quadratic function	Independent draw from the same generator

Appendix

Table 8: FoundPDE uses family-specific initial-condition generators adapted from the PROSE-PDE benchmark design [ 13 ] . Most equation families use sinusoidal superpositions, advection, wave, and Fokker–Planck use Gaussian-distribution initial conditions, and porous-medium equations use quadratic initial conditions. Training and in-domain testing draw independently from the same generator within each family.

Equation type	Solver
Heat, Klein–Gordon, Sine–Gordon, Porous medium, Cahn–Hilliard	Method of lines
Diffusion–reaction, Viscous conservation, Inviscid conservation	PDEBench solver pipeline
Advection, Wave	Analytical solution
KdV	Fourier spectral method
Fokker–Planck	Matrix numerical method

Appendix

Table 9: Numerical solvers used in our dataset follow PROSE-PDE [ 13 ] . Diffusion–reaction and conservation-law systems use the PDEBench pipeline [ 39 ] . KdV uses a Fourier spectral solver. Fokker–Planck uses a matrix numerical method. Advection and wave equations use exact solutions defined by the initial condition. The remaining families use a method-of-lines solver.

Setting	Eq. settings $\times$ ICs	$n_{L^{2}}/n_{c}$	TRR (%)	Mean relative $L^{2}$	Mean MRCSE	Median MRCSE
Clean in-domain and corruption control, $\eta=0$	$256\times 20$	$5115/4992$	$97.50\;[96.54,98.30]$	$0.0308\;[0.0241,0.0384]$	$0.185\;[0.052,0.393]$	$0.000705\;[0.000536,0.000855]$
Corruption, $\eta=0.02$	$256\times 20$	$5108/4833$	$94.43\;[92.64,96.02]$	$0.0358\;[0.0293,0.0435]$	$16.87\;[0.235,49.58]$	$0.00129\;[0.000910,0.00176]$
Corruption, $\eta=0.05$	$256\times 20$	$4868/4164$	$81.39\;[78.44,84.10]$	$0.0805\;[0.0700,0.0924]$	$2150\;[10.68,5709]$	$0.00430\;[0.00301,0.00615]$
Corruption, $\eta=0.1$	$256\times 20$	$3821/2640$	$52.09\;[48.20,56.00]$	$0.211\;[0.191,0.233]$	$8.78{\times}10^{7}\;[293,2.79{\times}10^{8}]$	$0.0206\;[0.0130,0.0302]$
Strict parameter OOD	$640\times 8$	$5011/4630$	$90.53\;[89.22,91.80]$	$0.156\;[0.139,0.174]$	$0.228\;[0.140,0.331]$	$0.00458\;[0.00395,0.00544]$

Appendix

Table 10: Clustered uncertainty intervals separate robust aggregate trends from heavy-tailed coefficient errors. Values in brackets are coefficient-setting-level 95% bootstrap intervals. TRR uses the complete 5,120-example denominator; mean relative

L^{2}

error uses reconstructable predictions; and MRCSE uses type-matched predictions with extractable packed coefficients. The column

n_{L^{2}}/n_{c}

reports these two evaluable counts. The widening mean-MRCSE intervals under stronger corruption and strict parameter OOD reflect a small number of extreme coefficient errors, whereas median MRCSE remains substantially more stable.

Family	Reference coefficient(s)	Strict OOD support
Heat	$c_{0}=3\times 10^{-3}$	$c/c_{0}\in[1.91,2.5]$
Advection	$\beta_{0}=1$	$\beta\in[-2.5,-1.91]\cup[1.91,2.5]$
KdV	$\delta_{0}=0.022$	$\delta/\delta_{0}\in[1.91,2.2]$
Fokker–Planck	$\eta_{0}=10^{-3}$	$\eta/\eta_{0}\in[1.91,2.0]$
Diffusion–reaction, logistic / linear	$\nu_{0}=3\times 10^{-3}$ , $\rho_{0}\in\{1,0.1\}$	Either $\nu/\nu_{0}\in[1.91,2.2]$ or $\rho/\rho_{0}\in[1.91,2.3]$
Diffusion–reaction, bistable / square-logistic	$\nu_{0}=3\times 10^{-3}$ , $\rho_{0}=1$	Either $\nu/\nu_{0}\in[1.91,2.0]$ or $\rho/\rho_{0}\in[1.91,2.2]$

Appendix

Table 11: Strictly disjoint coefficient supports define the parameter-OOD benchmark. Each listed interval lies outside the corresponding in-domain support in Table 5 . Multi-parameter equations shift one coefficient at a time unless a coupled restriction is stated. Family-specific upper limits and coupling constraints preserve numerical stability without changing the symbolic PDE family.

$\eta$	FoundPDE (recomputed)	FoundPDE (Sec. 3.3 )	SINDy	SINDy SG 11	SINDy SG 5	Weak SINDy $H=0.25$	Weak SINDy $H=0.375$
0	68	68	50	9	27	53	50
0.01	54	62	0	9	4	17	6
0.02	45	47	0	2	1	6	5
0.05	24	29	0	0	0	5	0
0.1	5	23	0	0	0	1	2

Appendix

Table 12: FoundPDE retains structural recovery at noise levels at which sparse-regression baselines fail. Number of the 80 diffusion-reaction trajectories (four equations with 20 initial conditions each) whose recovered terms match the ground truth. FoundPDE (recomputed) derives its derivative and Fourier channels from the noisy solution; FoundPDE (Sec. 3.3 ) uses noise injected into each channel of the clean augmented input. SG denotes Savitzky–Golay smoothing with the given window length, and

H

is the weak-form test-function half-width.

PDE Discovery

Partial differential equation discovery infers governing equations from observed fields to reveal mechanisms and build predictive physical models.

See all 36 papers

See, Hypothesize, Validate: Multimodal Agentic Framework for Discovering Governing PDEs
Aug 28, 2026Sarang Manoj Pekhale, Amartya Roy, Rajat Sarkar +1PDE Discovery Partial Differential Equations
Department of Applied Mechanics Indian Institute of Technology Delhi, India Yardi School of Artificial Intelligence (ScAI) Indian Institute of Technology Delhi, India
Symbolic Graph Networks for Robust PDE Discovery from Noisy Sparse Data
Mar 23, 2026Xingyu Chen, Junxiu An, Jun Guo +1PDE Discovery Graph Neural Networks
School of Software Engineering, Chengdu University of Information Technology, Chengdu, 610225, China · School of Statistics, Chengdu University of Information Technology, Chengdu, 610103, China · College of Applied Mathematics, Chengdu University of Information Technology, Chengdu, 610225, China +1
LLM-PDESR: Robust PDE Discovery via Subdomain Weighted Residuals and LLM-Guided Symbolic Hypothesis Generation
Jul 12, 2026Jinyang Du, Hao Ma, Xiaohu Shi +4PDE Discovery Large Language Model-Guided Scientific Discovery
College of Computer Science and Technology, Jilin University, Changchun, China · School of Big Data and Artificial Intelligence, Guangdong University of Finance and Economics, Guangzhou, China · Key Laboratory of Symbolic Computation and Knowledge Engineering of Ministry of Education, Jilin University, Changchun, China +2

Aug 28, 2026J/K move · Enter open · S save

See, Hypothesize, Validate: Multimodal Agentic Framework for Discovering Governing PDEs

Sarang Manoj Pekhale, Amartya Roy, Rajat Sarkar, Souvik Chakraborty

Department of Applied Mechanics Indian Institute of Technology Delhi, India Yardi School of Artificial Intelligence (ScAI) Indian Institute of Technology Delhi, India

PDE Discovery Partial Differential Equations AI Agents for Scientific Discovery Large Language Model-Guided Scientific Discovery

Open paper AskPDF arXiv

Discovering governing partial differential equations (PDEs) from observational data remains a core challenge across the sciences. Existing sparse-regression, symbolic-regression, and LLM-based approaches can be constrained by predefined libraries, noise sensitivity, hallucination, or limited iterative refinement. We introduce MAGE (Multimodal Agentic Governing Equation Discovery), an agentic framework that organizes PDE discovery as a confidence governed hypothesis validation loop inspired by the scientific cycle of observation, hypothesis, and falsification. Four role-specialized agents collaborate: a Differential Observer computing derivatives and diagnostic visualizations; a VLM-powered Phenomenology Extractor distilling qualitative cues from multimodal diagnostics; an LLM-driven Governing Law Synthesizer proposing candidates without a predefined library; and an Equation Arbiter fitting coefficients and assigning confidence scores. Discovery iterates until the top candidate clears a user-specified threshold, providing a structured process with an explicit accept-reject protocol. On the evaluated canonical PDE suite, MAGE obtains 8/8 exact structural recovery and the lowest coefficient error among the compared methods on 7/8 systems, with improvements of up to 4 orders of magnitude and a geometric-mean improvement of approximately 3 orders of magnitude. The pipeline also recovers the expected operators in two complex geometries and, on one laboratory sensor record, selects a cubic restoring-force model with held-out

R^2=0.98538

. These results support further study of structured agentic reasoning for library-free governing-law discovery, while broader generalization remains to be evaluated.

Abstract

PDE Discovery

Explore similar work

See, Hypothesize, Validate: Multimodal Agentic Framework for Discovering Governing PDEs

Symbolic Graph Networks for Robust PDE Discovery from Noisy Sparse Data

LLM-PDESR: Robust PDE Discovery via Subdomain Weighted Residuals and LLM-Guided Symbolic Hypothesis Generation