Weak-Form Evolutionary Kolmogorov-Arnold Networks for Solving Partial Differential Equations

Bongseok Kim; Jiahao Zhang; Guang Lin

cs.LGFeb 19, 2026

Weak-Form Evolutionary Kolmogorov-Arnold Networks for Solving Partial Differential Equations

Authors: Bongseok Kim, Jiahao Zhang, Guang Lin

Abstract

Partial differential equations (PDEs) form a central component of scientific computing. Among recent advances in deep learning, evolutionary neural networks have been developed to successively capture the temporal dynamics of time-dependent PDEs via parameter evolution. The parameter updates are obtained by solving a linear system derived from the governing equation residuals at each time step. However, strong-form evolutionary approaches can yield ill-conditioned linear systems due to pointwise residual discretization, and their computational cost scales unfavorably with the number of training samples. To address these limitations, we propose a weak-form evolutionary Kolmogorov-Arnold Network (KAN) for the scalable and accurate prediction of PDE solutions. We decouple the linear system size from the number of training samples through the weak formulation, leading to improved scalability compared to strong-form approaches. We also rigorously enforce boundary conditions by constructing the trial space with boundary-constrained KANs to satisfy Dirichlet and periodic conditions, and by incorporating derivative boundary conditions directly into the weak formulation for Neumann conditions. In conclusion, the proposed weak-form evolutionary KAN framework provides a stable and scalable approach for PDEs and contributes to scientific machine learning with potential relevance to future engineering applications.

Figures & tables

Panel of Figure 1 — Figure 1: Comparison between the Kolmogorov–Arnold Network (KAN) and a conventional fully connected DNN. (a) The KAN architecture employs functional mappings on each edge, enabling richer nonlinear representations and adaptive feature extraction. (b) The DNN architecture is composed of fully connected layers with activation functions.

Figure 2: Evolutionary Kolmogorov–Arnold network with Gaussian RBFs : Strong-form formulation. At each time step, the parameters ${W}(t)$ evolve by minimizing the PDE residual functional $\displaystyle\mathcal{J}(\gamma)=\tfrac{1}{2}\!\int_{\Omega}\!\big\|\tfrac{\partial\hat{u}}{\partial{W}}\gamma+\mathcal{N}(\hat{u})\big\|_{2}^{2}\,\mathrm{d}\mathbf{x}$ , where the optimal update direction satisfies $\displaystyle\frac{\partial{W}}{\partial t}=\arg\min_{\gamma}\mathcal{J}(\gamma)$ . After obtaining $\gamma_{\text{opt}}$ , the network parameters are explicitly updated in time using the forward Euler scheme $\displaystyle{W}(t_{n+1})={W}(t_{n})+\gamma_{\text{opt}}\Delta t$ , yielding a discrete-time parameter evolution associated with the strong-form PDE residual minimization.

Figure 3: Weak-form evolutionary Kolmogorov–Arnold network. At each time step, the parameters $W(t)$ are updated by minimizing the weak–form residual functional $\mathcal{J}(\gamma)=\tfrac{1}{2}\sum_{k=1}^{K}\!\left(\int_{\Omega}\big(\tfrac{\partial\hat{u}}{\partial W}\gamma+\mathcal{N}(\hat{u})\big)v_{k}(\mathbf{x})\,\mathrm{d}\mathbf{x}\right)^{2}$ , where $\{v_{k}\}_{k=1}^{K}$ denote the test functions. The optimal update direction is given by $\gamma_{\text{opt}}=\arg\min_{\gamma}\mathcal{J}(\gamma)$ , and the network parameters are advanced in time via the forward Euler step $W(t_{n+1})=W(t_{n})+\gamma_{\text{opt}}\Delta t$ , resulting in a discrete-time parameter evolution associated with the weak-form residual minimization.

	Strong form	Weak form
System to solve	$J^{T\,(P\times N_{\hat{u}})}J^{(N_{\hat{u}}\times P)}\gamma_{\text{opt}}=-J^{T\,(P\times N_{\hat{u}})}N^{(N_{\hat{u}}\times 1)}$	$J^{T\,(P\times K)}J^{(K\times P)}\gamma_{\text{opt}}=-J^{T\,(P\times K)}N^{(K\times 1)}$
Data–system coupling	Coupled : system size scales with $N_{\hat{u}}$	Decoupled : system size fixed by $K$
$J$ evaluation	$\mathcal{O}(N_{\hat{u}}\,P\,C_{\mathrm{AD}})$	$\mathcal{O}(N_{x}\,P\,C_{\mathrm{AD}}+KN_{x}P)$
$N$ evaluation	$\mathcal{O}(N_{\hat{u}}\,C_{\mathrm{PDE}}^{\mathrm{SF}})$	$\mathcal{O}(N_{x}\,C_{\mathrm{PDE}}^{\mathrm{WF}}+KN_{x})$
Least-squares system formation	$\mathcal{O}(N_{\hat{u}}P^{2})$	$\mathcal{O}(KP^{2})$
Total cumulative cost	$\mathcal{O}\!\left(M\left(N_{\hat{u}}P^{2}+N_{\hat{u}}PC_{\mathrm{AD}}+N_{\hat{u}}C_{\mathrm{PDE}}^{\mathrm{SF}}\right)\right)$	$\mathcal{O}\!\left(M\left(KP^{2}+N_{x}PC_{\mathrm{AD}}+N_{x}C_{\mathrm{PDE}}^{\mathrm{WF}}\right.\right.$ $\left.\left.+KN_{x}P+KN_{x}\right)\right)$

Table 1: Computational complexity comparison of strong-form and weak-form evolutionary solvers over

M

time steps. Here

J=\partial\hat{u}/\partial W

and

N=\mathcal{N}(\hat{u})

denote the parameter Jacobian and PDE residual.

N_{\hat{u}}

and

N_{x}

are the numbers of collocation and quadrature points,

K

the number of test functions, and

P

the number of trainable parameters.

C_{\mathrm{AD}}

denotes parameter automatic differentiation cost, while

C_{\mathrm{PDE}}

denotes spatial differentiation cost for residual evaluation. Due to integration by parts,

C_{\mathrm{PDE}}^{\mathrm{WF}}

is lower than

C_{\mathrm{PDE}}^{\mathrm{SF}}

	EvoKAN-WF	EvoKAN-SF	PINN-SF
Hidden layers	[3, 3, 3, 3]	[3, 3, 3, 3]	[15, 15, 15]
Activation functions	Gaussian RBFs/SiLU	Gaussian RBFs/SiLU	Tanh
Grid points number of activation functions	4	4	-
Number of trainable parameters	162	195	526
Optimizer	Adam	Adam	Adam/L-BFGS-B
Timestep	1e-07	1e-07	-

Table 2: Training configuration for the 1D Allen-Cahn equation (Eq. ( 41 )).

Figure 4: Comparison of predicted solutions for the one-dimensional Allen–Cahn equation (Eq. ( 41 )) at

t=2\times 10^{-5}

, trained with 100 data points. Both EvoKAN–WF and EvoKAN–SF closely match the ground-truth solution at this time, whereas the vanilla PINN does not reproduce the correct profile.

Figure 5: Comparison of predicted solutions for the one-dimensional Allen–Cahn equation (Eq. ( 41 )) at

t=5\times 10^{-5}

, trained with 100 data points. As the solution profile becomes steeper, EvoKAN–WF maintains good accuracy, whereas EvoKAN–SF shows small deviations near the center and boundaries, and the vanilla PINN fails to recover the correct solution.

Panel of Figure 6 — Figure 6: Strong-form (EvoKAN-SF) versus weak-form (EvoKAN-WF) evolutionary KAN solvers for the 1D Allen–Cahn equation (Eq. ( 41 )) with respect to the number of collocation points: (a) $L_{2}$ relative error, (b) condition number, and (c) wall-clock time. The weak-form approach achieves lower errors and substantially better conditioning, while avoiding the rapid growth in computational cost observed in the strong-form solver.

	$L_{2}$ relative error		Wall clock time
Number of collocation points	EvoKAN-SF	EvoKAN-WF	EvoKAN-SF [s]	EvoKAN-WF [s]	Speedup
64	4.2814e-02	7.1881e-03	7.4396 †	11.8412	0.63
100	2.6065e-02	7.1722e-03	8.5188 †	12.0379	0.71
1,000	8.3525e-03	7.1633e-03	16.1555	12.0761	1.34
10,000	8.9389e-03	7.1644e-03	33.5028	12.1654	2.76
50,000	9.0494e-03	6.8493e-03	340.6359	12.1727	27.99

Table 3:

L_{2}

relative error and wall clock time for the one-dimensional Allen–Cahn equation (Eq. ( 41 )) with varying the number of collocation points. Speedup is defined as

T_{\text{SF}}/T_{\text{WF}}

. Entries marked with

\dagger

correspond to EvoKAN-SF solutions whose

L_{2}

relative error remains above

\mathcal{O}(10^{-2})

, indicating insufficient accuracy.

	EvoKAN-WF	EvoKAN-SF	PINN-SF
Hidden layers	[16, 16, 16]	[16, 16, 16]	[25, 25, 25]
Activation functions	Gaussian RBFs/SiLU	Gaussian RBFs/SiLU	Tanh
Grid points number of activation functions	5	5	-
Number of trainable parameters	3392	3472	3670
Optimizer	Adam	Adam	Adam/L-BFGS-B
Timestep	1e-03	1e-03	-

Table 4: Training configuration for the 2D Burgers equation (Eq. ( 47 )).

Figure 7: Comparison of two-dimensional Burgers’ equation (Eq. ( 47 )) solutions at $t=0.1$ . All models are trained using $900$ collocation points. The top row shows the ground-truth FDM solution together with predictions from EvoKAN-WF, EvoKAN-SF, and PINN-SF. The bottom row presents the corresponding absolute error fields with respect to the ground truth. At this early time, EvoKAN-WF yields smaller errors than the other models while providing comparable representations of the overall flow structure.

Figure 8: Comparison of two-dimensional Burgers’ equation (Eq. ( 47 )) solutions at $t=0.4$ . All models are trained using $900$ collocation points. The top row shows the ground-truth FDM solution together with predictions from EvoKAN-WF, EvoKAN-SF, and PINN-SF. The bottom row presents the corresponding absolute error fields with respect to the ground truth. Despite the limited number of collocation samples, EvoKAN-WF provides more accurate predictions than the other models and better captures the flow structure across the domain.

Panel of Figure 9 — Figure 9: Strong-form (EvoKAN-SF) versus weak-form (EvoKAN-WF) evolutionary KAN solvers for the two-dimensional Burgers’ equation (Eq. ( 47 )) with respect to the number of collocation points: (a) $L_{2}$ relative error, (b) condition number, and (c) wall-clock time. The weak-form approach achieves lower errors and improved numerical conditioning, while avoiding the rapid growth in computational cost observed in the strong-form solver.

	$L_{2}$ relative error		Wall clock time
Number of collocation points	EvoKAN-SF	EvoKAN-WF	EvoKAN-SF [s]	EvoKAN-WF [s]	Speedup
900	5.67383e-02	2.27452e-02	10.7565 †	25.3241	0.42
1,600	4.01466e-02	2.25431e-02	28.1469	26.4252	1.07
2,500	4.27184e-02	2.17294e-02	57.7473	26.1492	2.21
10,000	2.53424e-02	2.15498e-02	617.0453	26.1723	23.58
40,000	2.17728e-02	2.16729e-02	2139.5068	29.4292	72.69

Table 5:

L_{2}

relative error and wall clock time for the two-dimensional Burgers’ equation (Eq. ( 47 )) with varying the number of collocation points. Speedup is defined as

T_{\text{SF}}/T_{\text{WF}}

. Entries marked with

\dagger

correspond to EvoKAN-SF solutions whose

L_{2}

relative error exceeds

5\%

, indicating insufficient accuracy.

	EvoKAN-WF	EvoKAN-SF	PINN-SF
Hidden layers	[4, 4, 4]	[4, 4, 4]	[15, 15, 15]
Activation functions	Gaussian RBFs/SiLU	Gaussian RBFs/SiLU	Tanh
Grid points number of activation functions	8	5	-
Number of trainable parameters	388	388	526
Optimizer	Adam	Adam	Adam/L-BFGS-B
Timestep	1e-03	1e-03	-

Table 6: Training configuration for the 2D Heat equation with nonlinear forcing term (Eq. ( 51 )).

Figure 10: Comparison of two-dimensional heat equation (Eq. ( 51 )) solutions at $t=0.20$ . All models are trained using $400$ collocation points. The top row shows the ground-truth solution together with predictions from EvoKAN-WF, EvoKAN-SF, and PINN-SF. The bottom row presents the corresponding absolute error fields with respect to the ground truth. The weak-form solver (EvoKAN-WF) exhibits smoother error distributions and reduced boundary artifacts compared to the strong-form approaches under the same training budget.

Figure 11: Comparison of two-dimensional heat equation solutions at

t=0.50

. All models are trained using

400

collocation points. The top row presents the ground-truth solution along with predictions obtained by EvoKAN-WF, EvoKAN-SF, and PINN-SF. The bottom row displays the corresponding absolute error fields relative to the ground truth. Even at this later time, the weak-form solver (EvoKAN-WF) maintains higher accuracy, exhibiting consistently lower errors and smoother spatial error patterns than the strong-form approaches under the same training budget.

Panel of Figure 12 — Figure 12: Strong-form (EvoKAN-SF) versus weak-form (EvoKAN-WF) evolutionary KAN solvers for the two-dimensional heat equation with respect to the number of collocation points: (a) $L_{2}$ relative error, (b) condition number, and (c) wall-clock time. The weak-form approach achieves lower errors and substantially better conditioning, while avoiding the rapid growth in computational cost observed in the strong-form solver.

Panel of Figure 13 — Figure 13: Gradient error on each boundary for the 2D heat equation (Eq. ( 51 )). The weak form evolutionary KAN progressively reduces the gradient error as the time evolution proceeds and the error approaches zero along all boundaries. The weak formulation provides stable treatment of the Neumann boundary condition and enables monotonic error reduction during network updates. The strong-form evolutionary network enforces a constant-derivative boundary condition, making convergence of the boundary error to zero difficult, while the strong-form PINN exhibits oscillatory boundary errors during the time evolution.

	$L_{2}$ relative error		Wall-clock time
Number of collocation points	EvoKAN-SF	EvoKAN-WF	EvoKAN-SF [s]	EvoKAN-WF [s]	Speedup
400	4.6195e-02	1.7171e-02	24.1326	22.8757	1.06
900	3.0303e-02	1.2976e-02	37.9921	25.1838	1.51
1,600	4.5767e-02	1.5302e-02	74.4028	26.6005	2.80
10,000	1.7737e-02	1.4792e-02	142.7417	31.2263	4.57
40,000	1.5728e-02	1.2892e-02	1111.4212	54.5341	20.39

Table 7:

L_{2}

relative error and wall-clock time for the two-dimensional heat equation (Eq. ( 51 )) with varying numbers of collocation points. Speedup is defined as

T_{\text{SF}}/T_{\text{WF}}

	EvoKAN-WF	EvoKAN-SF	PINN-SF
Hidden layers	[7, 7, 7]	[7, 7, 7]	[15, 15, 15]
Activation functions	Gaussian RBFs/SiLU	Gaussian RBFs/SiLU	Tanh
Grid points number of activation functions	5	5	-
Number of trainable parameters	770	770	526
Optimizer	Adam	Adam	Adam/L-BFGS-B
Timestep	1e-03	1e-03	-

Table 8: Training configuration for the 2D porous medium equation with drift (Eq. ( 57 )).

Figure 14: Solution field comparison for the two-dimensional PME with drift under periodic boundary conditions (Eq. ( 57 )). The solution fields at $t=0$ (top row) and $t=0.4$ (bottom row) are shown for EvoKAN-WF, EvoKAN-SF, and PINN, respectively. All models are trained using $N_{c}=3600$ collocation points, with a uniform time step $\Delta t=10^{-2}$ up to the final time $T=0.4$ . The color map represents the spatial distribution of the solution $u(x_{1},x_{2},t)$ , illustrating the nonlinear diffusion and drift-induced transport dynamics.

Panel of Figure 15 — Figure 15: Strong-form (EvoKAN-SF) versus weak-form (EvoKAN-WF) evolutionary KAN solvers for the two-dimensional porous medium equation (Eq. ( 57 )) with respect to the number of collocation points: (a) condition number and (b) wall-clock time. The weak-form approach achieves substantially better conditioning while avoiding the rapid growth in computational cost observed in the strong-form solver.

Panel of Figure 16 — Figure 16: Energy evolution of strong-form (EvoKAN-SF) and weak-form (EvoKAN-WF) evolutionary KAN solvers for the two-dimensional porous medium equation under different numbers of collocation points. For EvoKAN-SF, monotone energy decay is observed only when a sufficiently large number of collocation points is used. At lower resolutions, EvoKAN-SF exhibits energy increase or incomplete decay over time. In contrast, EvoKAN-WF shows consistent energy decay across all tested resolutions.

	Approximation validity		Wall clock time
Number of collocation points	EvoKAN-SF	EvoKAN-WF	EvoKAN-SF [s]	EvoKAN-WF [s]	Speedup
400	$\times$	$\bigcirc$	21.0145	20.6104	1.02
1,600	$\times$	$\bigcirc$	86.8148	24.9313	3.48
3,600	$\bigcirc$	$\bigcirc$	113.3219	26.1163	4.34
10,000	$\bigcirc$	$\bigcirc$	580.5136	31.6187	18.37
40,000	$\bigcirc$	$\bigcirc$	3599.2407	49.2156	73.12

Table 9: Solution validity and wall clock time for the 2D PME (Eq. ( 57 )) with varying the number of collocation points. The solution validity indicator reports whether the numerical solution exhibits a dissipative energy evolution consistent with the governing equation. Entries marked with

\times

correspond to cases where the computed energy shows non-decreasing behavior or incomplete decay over the simulated time interval, indicating that the solution does not provide an admissible approximation under the given discretization. Speedup is defined as

T_{\text{SF}}/T_{\text{WF}}

Explore similar work

Mar 24, 2026

Weak-PDE-Net: Discovering Open-Form PDEs via Differentiable Symbolic Networks and Weak Formulation

Discovering governing Partial Differential Equations (PDEs) from sparse and noisy data is a challenging issue in data-driven scientific computing. Conventional sparse regression methods often suffer from two major limitations: (i) the instability of numerical differentiation under sparse and noisy data, and (ii) the restricted flexibility of a pre-defined candidate library. We propose Weak-PDE-Net, an end-to-end differentiable framework that can robustly identify open-form PDEs. Weak-PDE-Net consists of two interconnected modules: a forward response learner and a weak-form PDE generator. The learner embeds learnable Gaussian kernels within a lightweight MLP, serving as a surrogate model that adaptively captures system dynamics from sparse observations. Meanwhile, the generator integrates a symbolic network with an integral module to construct weak-form PDEs, avoiding explicit numerical differentiation and improving robustness to noise. To relax the constraints of the pre-defined library, we leverage Differentiable Neural Architecture Search strategy during training to explore the functional space, which enables the efficient discovery of open-form PDEs. The capability of Weak-PDE-Net in multivariable systems discovery is further enhanced by incorporating Galilean Invariance constraints and symmetry equivariance hypotheses to ensure physical consistency. Experiments on several challenging PDE benchmarks demonstrate that Weak-PDE-Net accurately recovers governing equations, even under highly sparse and noisy observations.

Xinxin Li, Xingyu Cui, Jin Qi +3

School of Mathematical Sciences, East China Normal University, Shang Hai, 200241, China · Institute of Applied Physics and Computational Mathematics, 100094, Beijing, China · Shanghai Zhangjiang Institute of Mathematics, 201203, Shanghai, China +2

PDE Discovery Symbolic Regression

Ask PDF

Apr 16, 2026

Python library supporting Discrete Variational Formulations and training solutions with Collocation-based Robust Variational Physics Informed Neural Networks (DVF-CRVPINN)

We explore the possibility of solving Partial Differential Equations (PDEs) using discrete weak formulations. We propose a programming environment for defining a discrete computational domain, introducing discrete functions defined over a set of points, constructing discrete inner products, and introducing discrete weak formulations employing Kronecker delta test functions. Building on this setup, we propose a discrete neural network representation, training the solution function defined over a discrete set of points and employing discrete finite difference derivatives in the automatic differentiation procedures. As a challenging computational model example, we focus on Stokes equations in two-dimensions, defined over a discrete set of points. We train the solution using the discrete weak residual and the Adamax algorithm with discrete automatic differentiation of the discrete gradients. Despite introducing the python environment, we also provide a rigorous mathematical formulation based on discrete weak formulations, proving the well-posedness and robustness of the loss function. The solution of the discrete weak formulations is based on neural network training employing a robust loss function that is related to the true error. In this way, we have a robust control of the numerical error during the training of the neural networks. Besides the Stokes formulation, we also explain the functionality of the proposed library using the Laplace problem formulation.

Tomasz Służalec, Marcin Łoś, Askold Vilkha +1

AGH University of Krakow, Poland

PDE Solving Physics-Informed ML

Ask PDF

Jul 22, 2026

PG-KINN: A Physics-Informed Petrov-Galerkin Kolmogorov-Arnold Network for Solving Forward and Inverse PDEs

Physics-informed learning of partial differential equations (PDEs) has been dominated by multilayer perceptrons (MLPs), whose spectral bias and dense parameterization limit both accuracy and interpretability. Kolmogorov Arnold Networks (KANs) mitigate these limitations because their learnable spline activations are structurally aligned with the piecewise-polynomial bases of classical discretizations. However, the way a PDE is cast into a loss functional is as decisive as the choice of approximator: strong-form residual minimization requires high-order derivatives and heavily weighted losses, the energy (Bubnov-Galerkin) form is restricted to self-adjoint operators and, as we show, collapses to a trivial solution for parameter-identification problems, and boundary integral forms require a known fundamental solution. We propose PG-KINN, a physics-informed KAN built on a Petrov-Galerkin formulation in which the trial space is a KAN and the test space is an independent, compactly supported, piecewise-polynomial space evaluated with Gauss-Legendre quadrature. Integration by parts lowers the differentiation order while retaining applicability to general non-self-adjoint, nonlinear, and inverse problems; the localized test functions turn the global residual into a set of element-wise weak residuals with favorable conditioning. On a suite of benchmarks spanning crack singularities, stress concentration, Neo-Hookean hyperelasticity, inverse parameter identification in heterogeneous media, and complex geometries, PG-KINN consistently outperforms legacy MLP baselines and state-of-the-art KAN-based strong/energy/inverse formulations (PIKAN). These results position the Petrov-Galerkin coupling of KAN trial spaces and polynomial test spaces as a robust and accurate route for AI-based computational mechanics.

Amirhossein Sadr, Nima Soltani, Vahideh Moghtadaiee +3

School of Computer Science, Institute for Research in Fundamental Sciences (IPM), Tehran, Iran · Department of Computer Science and Engineering, Shahid Beheshti University, Tehran, Iran · Cyberspace Research Institute, Shahid Beheshti University, Tehran, Iran +1

Neural PDE Solvers Kolmogorov-Arnold Networks Physics-Informed ML Finite Element Method PDE Inverse Problems

Ask PDF