Algorithmically Designed Artificial Neural Networks (ADANNs): Higher
  order deep operator learning for parametric partial differential equations
In this article we propose a new deep learning approach to approximate
operators related to parametric partial differential equations (PDEs). In
particular, we introduce a new strategy to design specific artificial neural
network (ANN) architectures in conjunction with specific ANN initialization
schemes which are tailor-made for the particular approximation problem under
consideration. In the proposed approach we combine efficient classical
numerical approximation techniques with deep operator learning methodologies.
Specifically, we introduce customized adaptions of existing ANN architectures
together with specialized initializations for these ANN architectures so that
at initialization we have that the ANNs closely mimic a chosen efficient
classical numerical algorithm for the considered approximation problem. The
obtained ANN architectures and their initialization schemes are thus strongly
inspired by numerical algorithms as well as by popular deep learning
methodologies from the literature and in that sense we refer to the introduced
ANNs in conjunction with their tailor-made initialization schemes as
Algorithmically Designed Artificial Neural Networks (ADANNs). We numerically
test the proposed ADANN methodology in the case of several parametric PDEs. In
the tested numerical examples the ADANN methodology significantly outperforms
existing traditional approximation algorithms as well as existing deep operator
learning methodologies from the literature.
Arnulf Jentzen, Adrian Riekert, Philippe von Wurstemberger
# Algorithmically Designed Artificial

###### Abstract

In this article we propose a new deep learning approach to solve parametric partial differential equations (PDEs) approximately. In particular, we introduce a new strategy to design specific artificial neural network (ANN) architectures in conjunction with specific ANN initialization schemes which are tailor-made for the particular scientific computing approximation problem under consideration. In the proposed approach we combine efficient classical numerical approximation techniques such as _higher-order Runge-Kutta schemes_ with sophisticated deep (operator) learning methodologies such as the recently introduced _Fourier neural operators (FNOs)_. Specifically, we introduce customized adaptions of existing standard ANN architectures together with specialized initializations for these ANN architectures so that at initialization we have that the ANNs closely mimic a chosen efficient classical numerical algorithm for the considered approximation problem. The obtained ANN architectures and their initialization schemes are thus strongly inspired by numerical algorithms as well as by popular deep learning methodologies from the literature and in that sense we refer to the introduced ANNs in conjunction with their tailor-made initialization schemes as _Algorithmically Designed Artificial Neural Networks_ (ADANNs). We numerically test the proposed ADANN approach in the case of some parametric PDEs. In the tested numerical examples the ADANN approach significantly outperforms existing traditional approximation algorithms as well as existing deep learning methodologies from the literature.
###### Contents

* 1 Introduction
* 2 A rough overview of the ADANN approach
	* 2.1 Base model with highly specialized initializations
	* 2.2 Difference models
	* 2.3 Multiple runs over initialization and training
* 3 Derivation of base models for semilinear heat PDEs
	* 3.1 One-dimensional semilinear heat PDEs with Dirichlet boundary conditions
	* 3.2 Designing algorithms
		* 3.2.1 Spatial finite difference discretization
		* 3.2.2 Temporal linearly implicit Runge-Kutta discretizations
		* 3.2.3 A compact reformulation of the designing algorithms
	* 3.3 Designing the base model
* 4 Numerical simulations
	* 4.1 One-dimensional reaction diffusion type equation
	* 4.2 One-dimensional Sine-Gordon type equation
	* 4.3 Two-dimensional semilinear heat equation
* A Second order linearly implicit Runge-Kutta methods
* A.1 Order conditions for general LIRK methods
* A.2 A family of 2 stage linearly implicit Runge-Kutta methods of order 2
* A.3 The special case of the Crank-Nicolson explicit Euler method

## 1 Introduction

Deep learning approximation methods - usually consisting of deep artificial neural networks (ANN) trained through stochastic gradient descent (SGD) optimization methods - belong nowadays to the most heavily employed approximation methods in the digital world. The striking feature of deep learning methods is that in many situations numerical simulations suggest that the computational effort of such methods seem to grow only at most polynomially in the input dimension \(d\in\mathbb{N}=\{1,2,3,\dots\}\) of the problem under consideration. In contrast, classical numerical methods usually suffer under the so-called _curse of dimensionality_ (cf., e.g., Bellman [4], Novak & Wozniakowski [37, Chapter 1], and Novak & Wozniakowski [38, Chapter 9]) in the sense that the computational effort grows at least exponentially in the dimension.

In the recent years, deep learning technologies have also been intensively used to attack problems from scientific computing such as the numerical solutions of partial differential equations (PDEs). In particular, deep learning approximation methods have been used to approximately solve high-dimensional nonlinear PDEs (see, e.g., [25, 10, 11, 14, 16, 2, 42] and the references mentioned therein) such as high-dimensional nonlinear pricing problems from financial engineering and Hamiltonian-Jacobi-Bellman equations from optimal control. In the context of such high-dimensional nonlinear PDEs, the progress of deep learning approximation methods is obvious as there are - except of in some special cases (see, e.g., [19, 20, 36] and the references therein for Branching type methods and see, e.g., [11, 12, 13, 22] and the references therein for multilevel Picard methods) - essentially no alternative numerical approximation methods which are capable of solving such high-dimensional nonlinear PDEs.

There is nowadays also a huge literature on deep learning approximation methods for low-dimensional PDEs (cf., e.g., [24, 41]). For low-dimensional PDEs, in most cases, there usually
already exist a number of highly efficient traditional (non-deep learning based) approximation methods in the scientific literature (cf., e.g., [23, 43]). Nonetheless, there are several convincing arguments that deep learning approximation methods might have the potential to significantly outperform such efficient traditional approximation methods from classical numerics.

One situation where this strongly applies is in the context of _parametric_ PDE approximation problems. Specifically, in applications one is often not only interested to approximately solve the considered PDE models once but instead there is often the need to approximately solve such models again and again but with different initial values and/or different model parameters and the idea of deep learning approaches in this context is to try to not only solve one fixed PDE model but instead to use deep ANNs to learn the whole solution mapping which maps initial values and model parameters to the corresponding PDE solutions. In particular, even though the original PDE model has often only one to three space-dimensions, the associated parametric approximation problem becomes very high-dimensional due to the high number of parameters to approximately describe the initial value and the model parameters. Deep learning methods seem to be very natural candidates for such kind of problems in a way that the deep ANNs learn the mapping from parametrizations of the initial values and/or the model parameters to the PDE solution based on training data - the deep learning methods in this situation are then often referred to as _operator learning_ approaches (cf., e.g., [29, 30, 32]).

However, even though very remarkable advances have been accomplished in this direction of research, for instance, by means of so-called _Fourier neural operator_ (FNO) approximations (see Li et al. [30]), so far in most situations deep operator learning techniques do not outperform the most efficient higher order classical numerical methods for the considered approximation problem. This is also not entirely surprising due to fundamental lower bounds established in the literature that a wide class of methods, including typical deep learning approximations, can in general not overcome the curse of dimensionality in the \(L^{\infty}\)-norm (cf., e.g., Heinrich & Sindambiwe [18], Heinrich [17], Grohs & Voigtlander [15]).

It is precisely the objective of this work to introduce a new deep operator learning approximation approach which aims to overcome this challenge and outperforms highly efficient higher order classical numerical methods for the considered approximation problems. For this, we introduce a new strategy to design specific ANN architectures in conjunction with specific ANN initialization schemes which are tailor-made for the particular scientific computing approximation problem under consideration. In the proposed approach we combine efficient classical numerical approximation methods such as higher-order Runge-Kutta schemes with sophisticated deep (operator) learning techniques such as FNO approximations. The obtained ANN architectures and their initialization schemes are thus strongly inspired by numerical algorithms as well as by popular deep learning methodologies from the literature and, in that sense, we refer to the introduced ANNs in conjunction with their tailor-made initialization schemes as _Algorithmically Designed Artificial Neural Networks_ (ADANNs). We numerically test the proposed ADANN approach in the case of some parametric PDEs. In the tested numerical examples the ADANN approach significantly outperforms existing traditional approximation algorithms as well as existing deep learning methodologies from the literature.

We now discuss some ideas in the scientific literature which are related to the ADANN approach introduced in this paper. The ADANN technology is partially inspired by the _learning the random variable_ methodology in Becker et al. [3] where so-called Monte Carlo neural networks have been introduced which have the property that at initialization the realizations of those networks correspond to sample realizations of Monte Carlo algorithms. Further approaches related to the here proposed ADANN methodology can, e.g., be found in [1, 8, 34, 39] where certain parameters of classical numerical approximation methods such as finite difference methods and higher-order Runge-Kutta schemes have been improved through a training process.

Next we mention several other promising deep operator learning approaches for learning the solution operators associated to PDEs. One of the most successful methods in practice are the
FNOs introduced in Li et al. [30]. The derivation of FNOs is based on Li et al. [29], an earlier paper by the same authors, where so-called graph kernel networks are employed. In [29] each layer of the network represents a convolution with a kernel computed by a neural network, which is replaced by a multiplication in Fourier space in [30]. The article Li et al. [28] generalizes FNOs to more complicated geometries. In Brandstetter et al. [6] the FNO methodology is extended by using Clifford layers, where calculations are performed in higher-dimensional non-commutative Clifford algebras. Another successful approach is the deep operator network (DeepONet) architecture introduced in Lu et al. [32], which consists of two types of DNNs that take as input the output space points and the input function values, respectively. For a comparison between the DeepONet and FNO methodologies we refer to Lu et al. [33]. Generalizing DeepONets, the work Lanthaler et al. [27] uses a more sophisticated nonlinear DNN architecture for operator learning. In Pham & Warin [40] operators on Wasserstein spaces, for example, mean-field interactions of measures, are learned using networks based on standard DNNs and DeepONets. The article Nelsen & Stuart [35] uses random feature maps associated to operator-valued kernels to approximate operators between Banach spaces. The paper Liu et al. [31] approximates the entire flow map associated to ODEs by training a different DNN in each time-step and combining these architectures with classical Runge-Kutta methods on different time scales. We also refer to [7, 26] for estimates for approximation and generalization errors in network-based operator learning for PDEs. Finally, we refer to [6, Appendix D] for a much more detailed literature overview regarding different variants of FNOs and other neural network architectures for solving PDEs.

The remainder of this article is organized as follows. In Section 2 we introduce the main ideas of the ADANN methodology in an abstract setting. In Section 3 we describe the ADANN methodology in detail in the case of semilinear heat PDEs. In Section 4 we present three numerical simulations comparing the ADANN methodology to existing methods in the literature.

## 2 A rough overview of the ADANN approach

In this section we describe the main ideas of the ADANN methodology for a class of representative approximation problems. Specifically, we consider the problem of numerically approximating a function

\[\mathcal{S}\colon\mathcal{I}\to\mathcal{O} \tag{1}\]

where \(\mathcal{I}\) and \(\mathcal{O}\) are topological spaces1. Roughly speaking, the ADANN methodology proposed in this paper relies on the following three main ingredients:

Footnote 1: For example, we think of the mapping which assigns to all suitable initial conditions of an initial value PDE problem of evolutionary type the corresponding solution of the PDE at the terminal time.

1. _Base model with highly specialized initializations_: Based on a family of classical (higher order) numerical approximation algorithms we design a tailor-made problem-specific ANN type model and a corresponding family of highly specialized initializations for the model and train the model to approximate \(\mathcal{S}\) (cf. Section 2.1).
2. _Difference model_: We employ existing ANN technologies from the literature to approximately learn the difference between realizations of the base model of step (i) and the solution operator \(\mathcal{S}\). Adding the difference model to the base model results in the _full ADANN model_ (cf. Section 2.2).
3. _Multiple runs over initialization and training_: We use an additional suitable optimization approach aiming to minimize the overall approximation error of the full ADANN model to repeat the training of the base model with different (possibly random) highly specialized initializations and the training of the difference model with different random initializations (cf. Section 2.3).
### Base model with highly specialized initializations

To derive the base model we consider a parametric family \((\Phi_{p})_{p\in\mathfrak{P}}\) of classical numerical approximation algorithms2 of the solution operator \(\mathcal{S}\colon\mathcal{I}\to\mathcal{O}\) indexed over a parameter set \(\mathfrak{P}\) and given for every \(p\in\mathfrak{P}\) by a mapping

Footnote 2: For example we think of a family of Runge-Kutta and/or finite element methods

\[\Phi_{p}\colon\mathcal{I}\to\mathcal{O}. \tag{2}\]

We then design an ANN type model \(\mathscr{B}\colon\mathbb{R}^{\mathrm{d}_{\mathrm{Base}}}\times\mathcal{I}\to \mathcal{O}\) with \(\mathbf{d}_{\mathrm{Base}}\in\mathbb{N}\) trainable parameters which can reproduce all the approximation algorithms in (2) as realizations of the model in the sense that for every \(p\in\mathfrak{P}\) there exists a parameter vector \(\mathbf{W}_{p}\in\mathbb{R}^{\mathrm{d}_{\mathrm{Base}}}\) such that

\[\mathscr{B}(\mathbf{W}_{p},\,\cdot\,)=\Phi_{p}\approx\mathcal{S}. \tag{3}\]

We refer to the function \(\mathscr{B}\) as the base model, for every choice of a parameter vector \(W\in\mathbb{R}^{\mathrm{d}_{\mathrm{Base}}}\) we call the function \(\mathscr{B}(W,\,\cdot\,)\) a realization of the base model, and we refer to the algorithms in (2) as the designing algorithms for the base model \(\mathscr{B}\).

Note that (3) shows that we already have parameters for the base model which yield reasonable approximations of \(\mathcal{S}\). In the next step we propose to train the base model with SGD type optimization methods to improve these approximations. In order to generate training data we assume that we have a random input value \(\mathfrak{I}\colon\Omega\to\mathcal{I}\) on a probability space \((\Omega,\mathcal{F},\mathbb{P})\) and an algorithm \(\Psi\colon\mathcal{I}\to\mathcal{O}\) which approximates \(\mathcal{S}\) with a very high accuracy (and potentially very high computational effort3) in the sense that the approximation

Footnote 3: For example, one can think of \(\Psi\) as being of the same class of algorithms as the algorithms employed in (2) but with much finer discretizations.

\[\Psi\approx\mathcal{S} \tag{4}\]

is much more accurate than the approximations in (3). Moreover, in order to describe the training objective we introduce a seminorm \(\left|\!\left|\!\left|\cdot\right|\!\right|\!\right|\colon\mathcal{O}\to[0,\infty)\) and use it to define the loss function \(L_{\mathrm{Base}}\colon\mathbb{R}^{\mathrm{d}_{\mathrm{Base}}}\times\mathcal{ I}\to[0,\infty)\) and the objective function \(\mathbf{L}_{\mathrm{Base}}\colon\mathbb{R}^{\mathrm{d}_{\mathrm{Base}}}\to[0,\infty]\) by imposing for all \(W\in\mathbb{R}^{\mathrm{d}_{\mathrm{Base}}}\), \(i\in\mathcal{I}\) that

\[L_{\mathrm{Base}}(W,i)=\left|\!\left|\!\left|\mathscr{B}(W,i)-\Psi(i)\right|\! \right|\!\right|^{2}\qquad\text{and}\qquad\mathbf{L}_{\mathrm{Base}}(W)= \mathbb{E}[L_{\mathrm{Base}}(W,\mathfrak{I})]. \tag{5}\]

We then propose to minimize \(\mathbf{L}_{\mathrm{Base}}\) by means of SGD type processes. As starting values for the SGD type processes we use some of the (possibly randomly chosen) parameters in (3) which result already at initialization in small expected loss values.

### Difference models

Once satisfactory base model parameters \(\mathfrak{W}\in\mathbb{R}^{\mathrm{d}_{\mathrm{Base}}}\) have been found in the previous step (see Section 2.1) we propose to employ existing deep learning strategies from the literature to approximate a scaled up difference between the reference solution \(\Psi\) and the realization of the base model \(\mathscr{B}(\mathfrak{W},\,\cdot\,)\). More precisely, to scale the error of the base model we use a Monte Carlo method to get an estimate \(\varepsilon\in(0,\infty)\) of the square root of the expected loss

\[\mathbf{L}_{\mathrm{Base}}(\mathfrak{W})=\left(\mathbb{E}\big{[}\!\left|\! \left|\mathscr{B}(\mathfrak{W},\mathfrak{I})-\Psi(\mathfrak{I})\right|\! \right|\!\right|^{2}\big{]}\right)^{\nicefrac{{1}}{{2}}}\approx\varepsilon \tag{6}\]

and then introduce an ANN structure \(\mathscr{D}\colon\mathbb{R}^{\mathrm{d}_{\mathrm{Diff}}}\times\mathcal{I}\to \mathcal{O}\) with \(\mathbf{d}_{\mathrm{Diff}}\in\mathbb{N}\) trainable parameters for which we aim to find a parameter vector \(\theta\in\mathbb{R}^{\mathrm{d}_{\mathrm{Diff}}}\) such that

\[\mathscr{D}(\theta,\,\cdot\,)\approx\tfrac{1}{\varepsilon}(\Psi-\mathscr{B}( \mathfrak{W},\,\cdot\,)). \tag{7}\]
We refer to the function \(\mathscr{D}\) as the difference model. To train this model we define the loss function \(L_{\text{Diff}}\colon\mathbb{R}^{\text{d}_{\text{Base}}}\times\mathbb{R}^{\text{ d}_{\text{Diff}}}\times(0,\infty)\times\mathcal{I}\to[0,\infty)\) and the objective function \(\mathbf{L}_{\text{Diff}}\colon\mathbb{R}^{\text{d}_{\text{Diff}}}\to[0,\infty]\) by requiring for all \(W\in\mathbb{R}^{\text{d}_{\text{Base}}}\), \(\theta\in\mathbb{R}^{\text{d}_{\text{Diff}}}\), \(v\in(0,\infty)\), \(i\in\mathcal{I}\) that

\[L_{\text{Diff}}(W,\theta,v,i)=\big{\|}\big{\|}\mathscr{D}(\theta,i )-\big{(}\tfrac{1}{v}(\Psi(i)-\mathscr{B}(W,i))\big{)}\big{\|}\big{\|}^{2}\quad \text{and}\quad\mathbf{L}_{\text{Diff}}(\theta)=\mathbb{E}\big{[}L_{\text{Diff} }(\mathfrak{W},\theta,\varepsilon,\mathfrak{I})\big{]} \tag{8}\]

and suggest to use SGD type methods to minimize \(\mathbf{L}_{\text{Diff}}\).

Combining the base model of Section 2.1 and the difference model of this subsection, we define the full ADANN model \(\mathscr{A}\colon\mathbb{R}^{\text{d}_{\text{Base}}}\times\mathbb{R}^{\text{ d}_{\text{Diff}}}\times\mathbb{R}\times\mathcal{I}\to\mathcal{O}\) by imposing for all \(W\in\mathbb{R}^{\text{d}_{\text{Base}}}\), \(\theta\in\mathbb{R}^{\text{d}_{\text{Diff}}}\), \(\epsilon\in\mathbb{R}\), \(i\in\mathcal{I}\) that

\[\mathscr{A}(W,\theta,\epsilon,i)=\mathscr{B}(W,i)+\epsilon \mathscr{D}(\theta,i). \tag{9}\]

### Multiple runs over initialization and training

The results of SGD type methods can strongly depend on the initial parameters and (3) suggests many choices of good initial parameters for the base model. In view of this, we propose to run through the trainings of the models described in Section 2.1 and Section 2.2 several times with different initialization parameters for the base and the difference model and use the best run as final approximation of the mapping \(\mathcal{S}\).

To make this procedure more clear, we now describe it in a simplified setting in detail. Specifically, let \(R\in\mathbb{N}\) be the number of runs. **For every run**\(r\in\{1,2,\dots,R\}\) we consider a random designing algorithm parameter \(\mathfrak{p}_{r}\colon\Omega\to\overline{\mathfrak{P}}\) and **train the base model** with an SGD process \(\mathcal{W}^{(r)}=(\mathcal{W}^{(r)}_{n})_{n\in\mathbb{N}_{0}}\colon\mathbb{N }_{0}\times\Omega\to\mathbb{R}^{\text{d}_{\text{Base}}}\) starting at \(\mathcal{W}^{(r)}_{0}=\mathbf{W}_{\mathfrak{p}_{r}}\) and proceeding for all \(n\in\mathbb{N}\) by

\[\mathcal{W}^{(r)}_{n}=\mathcal{W}^{(r)}_{n-1}-\frac{\gamma_{\text{ Base}}}{B_{\text{Base}}}\Big{[}{\sum_{b=1}^{B_{\text{Base}}}}(\nabla_{W}L_{ \text{Base}})(\mathcal{W}^{(r)}_{n-1},\mathfrak{I}^{(r,n,b)}_{\text{Base}}) \Big{]}, \tag{10}\]

where \(\gamma_{\text{Base}}\in(0,\infty)\) is the learning rate of the SGD method, where \(B_{\text{Base}}\in\mathbb{N}\) is the batch size of the SGD method, and where \(\mathfrak{I}^{(r,n,b)}_{\text{Base}}\colon\Omega\to\mathcal{I}\), \(n,b\in\mathbb{N}\), are i.i.d. copies of the random input \(\mathfrak{I}\). We stop these SGD processes after \(N_{\text{Base}}\in\mathbb{N}\) training steps.

Next, for every run \(r\in\{1,2,\dots,R\}\) we use a Monte Carlo method to get an estimate \(\varepsilon^{(r)}\colon\Omega\to(0,\infty)\) of the square root of the expected loss

\[\big{(}\mathbb{E}\big{[}\big{\|}\mathscr{B}(W,\mathfrak{I})- \Psi(\mathfrak{I})\big{\|}^{2}\big{]}\big{)}^{\nicefrac{{1}}{{2}}}\big{|}_{W= \mathcal{W}^{(r)}_{N_{\text{Base}}}}\approx\varepsilon^{(r)} \tag{11}\]

and then **train the difference model** starting with some standard initialization from the literature with the SGD process \(\Theta^{(r)}=(\Theta^{(r)}_{n})_{n\in\mathbb{N}_{0}}\colon\mathbb{N}_{0}\times \Omega\to\mathbb{R}^{\text{d}_{\text{Base}}}\) satisfying for all \(n\in\mathbb{N}\) that

\[\Theta^{(r)}_{n}=\Theta^{(r)}_{n-1}-\frac{\gamma_{\text{Diff}}}{B_ {\text{Diff}}}\Big{[}{\sum_{b=1}^{B_{\text{Diff}}}}(\nabla_{\theta}L_{\text{ Diff}})(\mathcal{W}^{(r)}_{N_{\text{Base}}},\Theta^{(r)}_{n-1},\varepsilon^{(r)}, \mathfrak{I}^{(r,n,b)}_{\text{Diff}})\Big{]}, \tag{12}\]

where \(\gamma_{\text{Diff}}\in(0,\infty)\) is the learning rate of the SGD method, where \(B_{\text{Diff}}\in\mathbb{N}\) is the batch size of the SGD method, and where \(\mathfrak{I}^{(r,n,b)}_{\text{Diff}}\colon\Omega\to\mathcal{I}\), \(n,b\in\mathbb{N}\), are i.i.d. copies of the random input \(\mathfrak{I}\). We stop these SGD processes after \(N_{\text{Diff}}\in\mathbb{N}\) steps.

Finally, we approximately **calculate the optimal run**\(\mathbf{r}\in\{1,2,\dots,R\}\) such that the expected loss

\[\begin{split}&\Big{(}\mathbb{E}\Big{[}\big{\|}\big{\|}\mathscr{A} \big{(}W,\theta,\epsilon,\mathfrak{I}\big{)}-\Psi(\mathfrak{I})\big{\|}^{2} \Big{]}\Big{)}\Big{|}_{(W,\theta,\epsilon)=(\mathcal{W}^{(r)}_{N_{\text{Base}}}, \Theta^{(r)}_{N_{\text{Diff}}},\varepsilon^{(r)})}\\ &=\min_{r\in\{1,2,\dots,R\}}\Big{(}\mathbb{E}\Big{[}\big{\|} \big{\|}\mathscr{A}\big{(}W,\theta,\epsilon,\mathfrak{I}\big{)}-\Psi( \mathfrak{I})\big{\|}^{2}\Big{]}\Big{)}\Big{|}_{(W,\theta,\epsilon)=( \mathcal{W}^{(r)}_{N_{\text{Base}}},\Theta^{(r)}_{N_{\text{Diff}}},\varepsilon^{( r)})}\end{split} \tag{13}\]
of the full ADANN model \(\mathscr{A}\) in (9) is minimal and, thereafter, we propose to use

\[\mathscr{A}\big{(}\mathcal{W}^{(\mathbf{r})}_{N_{\text{Base}}},\Theta^{(\mathbf{r })}_{N_{\text{Diff}}},\varepsilon^{(\mathbf{r})},\,\cdot\,\big{)}=\mathscr{B}( \mathcal{W}^{(\mathbf{r})}_{N_{\text{Base}}},\,\cdot\,)+\varepsilon^{(\mathbf{ r})}\mathscr{D}(\Theta^{(\mathbf{r})}_{N_{\text{Diff}}},\,\cdot\,)\approx \mathcal{S} \tag{14}\]

as an approximation of \(\mathcal{S}\) in (1).

In the procedure described above we made several major simplifications when compared to the methodology used in our numerical simulations in Section 4. First, in our numerical simulations we train the models with sophisticated SGD type methods such as the ADAM optimizer with adaptive learning rates and adaptive batch sizes. However, for simplicity in this section we only described the case of the plain-vanilla SGD method with constant learning rates and constant batch sizes. Second, we did not specify above how to sample the random designing algorithm parameter vectors \(\mathfrak{p}_{r}\colon\Omega\to\mathfrak{P}\), \(r\in\{1,2,\ldots,R\}\), used to initialize the base model. In our numerical simulations we used different approaches to this end such as, e.g., deterministically exploring the parameter set \(\mathfrak{P}\) or drawing random samples from a uniform distribution over the parameter set \(\mathfrak{P}\). Finally, above we restricted ourselves to the case where in each run we train exactly one base network and one corresponding difference network, whereas in some of our numerical simulations, we use different approaches to decide for each run whether to train a new base model and a corresponding difference model or whether to only train a new difference model for an already trained base model.

## 3 Derivation of base models for semilinear heat PDEs

In this section we describe a possible way to design and initialize base models for the problem of approximating the mapping from the initial condition of a semilinear heat PDE to the terminal value. We will use this base model in our numerical simulations in Sections 4.1 and 4.2.

We describe the approximation problem for semilinear heat PDEs in Section 3.1, we introduce the designing algorithms in Section 3.2, and we present the resulting base model and its tailor-made problem specific initializations in Section 3.3.

### One-dimensional semilinear heat PDEs with Dirichlet boundary conditions

We now introduce the setting for the approximation problem considered in this section. Let \(T\in(0,\infty)\), \(f\in C(\mathbb{R},\mathbb{R})\), \(\mathcal{I}\subseteq C([0,1],\mathbb{R})\) and for every \(g\in\mathcal{I}\) let \(u_{g}\in C^{1,2}([0,T]\times[0,1],\mathbb{R})\) satisfy for all \(t\in[0,T]\), \(x\in[0,1]\) that

\[\big{(}\tfrac{\partial}{\partial t}u_{g}\big{)}(t,x)=(\Delta_{x}u_{g})(t,x)+f( u_{g}(t,x)),\quad u_{g}(0,x)=g(x),\quad\text{and}\quad u_{g}(t,0)=u_{g}(t,1)=0. \tag{15}\]

Our goal is to approximate the map \(\mathcal{S}\colon\mathcal{I}\to C([0,1],\mathbb{R})\) given for all \(g\in\mathcal{I}\) by \(\mathcal{S}(g)=u_{g}(T,\cdot)\).

### Designing algorithms

We now derive a family of approximation algorithms for the semilinear heat PDE in (15) which will serve as designing algorithms for the base model derived in this section. The algorithms are based on discretizing the space with the finite difference method and on discretizing the time with a family of second order linearly implicit Runge-Kutta (LIRK) methods.

#### 3.2.1 Spatial finite difference discretization

For the spatial discretization of the PDE in (15) we use \(N\in\mathbb{N}\) equidistant points \(\mathfrak{x}_{1},\mathfrak{x}_{2},\ldots,\mathfrak{x}_{N}\in[0,1]\) given for all \(i\in\{1,2,\ldots,N\}\) by \(\mathfrak{x}_{i}=\frac{i}{N+1}\) and consider the corresponding finite difference
discretization of the Laplace operator with Dirichlet boundary conditions given by

\[A=(1+N)^{2}\begin{pmatrix}-2&1&0&0&\cdots&0&0&0\\ 1&-2&1&0&\cdots&0&0&0\\ 0&1&-2&1&\cdots&0&0&0\\ &&&\ddots&&\\ 0&0&0&0&\cdots&1&-2&1\\ 0&0&0&0&\cdots&0&1&-2\end{pmatrix}\in\mathbb{R}^{N\times N}. \tag{16}\]

Note that for all \(g\in C^{2}([0,1],\mathbb{R})\) with \(g(0)=g(1)=0\) we have that

\[A\begin{pmatrix}g(\mathfrak{x}_{1})\\ g(\mathfrak{x}_{2})\\ \vdots\\ g(\mathfrak{x}_{N})\end{pmatrix}\approx\begin{pmatrix}(\Delta g)(\mathfrak{x} _{1})\\ (\Delta g)(\mathfrak{x}_{2})\\ \vdots\\ (\Delta g)(\mathfrak{x}_{N})\end{pmatrix}. \tag{17}\]

Using this spatial discretization on the PDE in (15) results in an initial value problem. Formally, for every \(\mathfrak{g}\in\mathbb{R}^{N}\) we let \(\mathbf{u}_{\mathfrak{g}}\in C^{1}([0,T],\mathbb{R}^{N})\) satisfy4 for all \(t\in[0,T]\) that

Footnote 4: Throughout this paper for every \(h\colon\mathbb{R}\to\mathbb{R}\), \(n\in\mathbb{N}\), \(x=(x_{1},x_{2},\ldots x_{n})\in\mathbb{R}^{n}\) we denote by \(h(x)\in\mathbb{R}^{n}\) the vector given by \(h(x)=(h(x_{1}),h(x_{2}),\ldots h(x_{n}))\).

\[\big{(}\tfrac{\partial}{\partial t}\mathbf{u}_{\mathfrak{g}}\big{)}(t)=A \mathbf{u}_{\mathfrak{g}}(t)+f(\mathbf{u}_{\mathfrak{g}}(t))\qquad\text{and} \qquad\mathbf{u}_{\mathfrak{g}}(0)=\mathfrak{g}. \tag{18}\]

If \(N\) is chosen sufficiently large we expect for all \(g\in\mathcal{I}\) that

\[\mathbf{u}_{(g(\mathfrak{x}_{1}),g(\mathfrak{x}_{2}),\ldots,g(\mathfrak{x}_{N }))}(T)\approx\begin{pmatrix}u_{g}(T,\mathfrak{x}_{1})\\ u_{g}(T,\mathfrak{x}_{2})\\ \vdots\\ u_{g}(T,\mathfrak{x}_{N})\end{pmatrix}\approx\begin{pmatrix}\mathcal{S}(g)( \mathfrak{x}_{1})\\ \mathcal{S}(g)(\mathfrak{x}_{2})\\ \vdots\\ \mathcal{S}(g)(\mathfrak{x}_{N})\end{pmatrix}. \tag{19}\]

#### 3.2.2 Temporal linearly implicit Runge-Kutta discretizations

In the next step we use a parametric family of second order LIRK methods5 to discretize the ODE in (18). For all parameters \(p=(p_{1},p_{2})\in(0,\infty)^{2}\) and all step sizes \(h\in[0,\infty)\) we let the LRIK time step \(\phi_{p}^{h}\colon\mathbb{R}^{N}\to\mathbb{R}^{N}\) satisfy for all \(U,k_{1},k_{2}\in\mathbb{R}^{N}\) with

Footnote 5: We refer to Section A for a derivation of this family

\[k_{1}=(I_{N}-hp_{2}A)^{-1}(AU+f(U))\qquad\text{and} \tag{20}\]

\[k_{2}=(I_{N}-hp_{2}A)^{-1}\big{(}A(U+h2p_{1}(\tfrac{1}{2}-p_{2})k_{1})+f(U+hp_{ 1}k_{1})\big{)} \tag{21}\]

that

\[\phi_{p}^{h}(U)=U+h\Big{[}(1-\tfrac{1}{2p_{1}})k_{1}+(\tfrac{1}{2p_{1}})k_{2} \Big{]}, \tag{22}\]

where \(I_{N}\in\mathbb{R}^{N\times N}\) is the identity matrix. For every number of time steps \(\mathbf{m}\in\mathbb{N}\) and every choice of parameters \(p\in(0,\infty)^{2}\) the corresponding designing algorithm \(\Phi_{p}^{\mathbf{m}}\colon\mathbb{R}^{N}\to\mathbb{R}^{N}\) is subsequently given by

\[\Phi_{p}^{\mathbf{m}} \tag{23}\]

If \(N\) is chosen sufficiently large we then expect for all parameters \(p\in(0,\infty)^{2}\) and all sufficiently large \(M\in\mathbb{N}\) that for all \(g\in\mathcal{I}\) we have that

\[\Phi_{p}^{M}(g(\mathfrak{x}_{1}),g(\mathfrak{x}_{2}),\ldots,g(\mathfrak{x}_{N }))\approx\mathbf{u}_{(g(\mathfrak{x}_{1}),g(\mathfrak{x}_{2}),\ldots,g( \mathfrak{x}_{N}))}(T)\approx\begin{pmatrix}\mathcal{S}(g)(\mathfrak{x}_{1})\\ \mathcal{S}(g)(\mathfrak{x}_{2})\\ \vdots\\ \mathcal{S}(g)(\mathfrak{x}_{N})\end{pmatrix}. \tag{24}\]
#### 3.2.3 A compact reformulation of the designing algorithms

To make the designing algorithms of Section 3.2.2 amenable to be written as realizations of an ANN type base model we now present a more compact reformulations of the algorithms in (23). For this we fix the number of time steps \(M\in\mathbb{N}\) and the corresponding time step size \(H=T/M\) and for every \(p=(p_{1},p_{2})\in(0,\infty)^{2}\) we let \(\mathbf{W}_{p}=(\mathbf{W}_{p,i})_{i\in\{1,2,\ldots,5\}}\in(\mathbb{R}^{N\times N })^{5}\) satisfy

\[\mathbf{W}_{p,1}=(I_{N}-Hp_{2}A)^{-1}(I_{N}+H(1-p_{2})A)+H^{2}( \tfrac{1}{2}-p_{2})\big{[}(I_{N}-Hp_{2}A)^{-1}A\big{]}^{2}, \tag{25}\] \[\mathbf{W}_{p,2}=H(1-\tfrac{1}{2p_{1}})(I_{N}-Hp_{2}A)^{-1}+H^{2} (\tfrac{1}{2}-p_{2})(I_{N}-Hp_{2}A)^{-1}A(I_{N}-Hp_{2}A)^{-1},\] (26) \[\mathbf{W}_{p,3}=H(\tfrac{1}{2p_{1}})(I_{N}-Hp_{2}A)^{-1},\qquad \mathbf{W}_{p,4}=(I_{N}-Hp_{2}A)^{-1}(I_{N}+H(p_{1}-p_{2})A),\] (27) \[\text{and}\qquad\mathbf{W}_{p,5}=Hp_{1}(I_{N}-Hp_{2}A)^{-1}. \tag{28}\]

Note that for all \(p\in(0,\infty)^{2}\), \(U\in\mathbb{R}^{N}\) we have that

\[\phi_{p}^{H}(U)=\mathbf{W}_{p,1}U+\mathbf{W}_{p,2}f(U)+\mathbf{W} _{p,3}f\big{(}\mathbf{W}_{p,4}U+\mathbf{W}_{p,5}f(U)\big{)}. \tag{29}\]

### Designing the base model

Roughly speaking, we propose to design the base model by considering the matrices in (29) as trainable parameters. Specifically, let \(\mathscr{B}\colon((\mathbb{R}^{N\times N})^{5})^{M}\times\mathcal{I}\to \mathbb{R}^{N}\) satisfy for all \(W=((W_{m,i})_{i\in\{1,2,\ldots,5\}})_{m\in\{1,2,\ldots,M\}}\in((\mathbb{R}^{N \times N})^{5})^{M}\), \(g\in\mathcal{I}\), \(U_{0},U_{1},\ldots,U_{M}\in\mathbb{R}^{N}\) with \(U_{0}=(g(\mathfrak{x}_{1}),g(\mathfrak{x}_{2}),\)\(\ldots,g(\mathfrak{x}_{N}))\) and

\[\forall\,m\in\{1,2,\ldots,M\}\colon\qquad U_{m}=W_{m,1}U_{m-1}+W_ {m,2}f(U_{m-1})\\ +W_{m,3}f\big{(}W_{m,4}U_{m-1}+W_{m,5}f(U_{m-1})\big{)} \tag{30}\]

that

\[\mathscr{B}(W,g)=U_{M}. \tag{31}\]

Note that if \(N\) and \(M\) were chosen sufficiently large (24), (29), (30), and (31) imply that for all \(p\in(0,\infty)^{2}\), \(g\in\mathcal{I}\) it holds that

\[\mathscr{B}\big{(}(\underbrace{\mathbf{W}_{p},\ldots,\mathbf{W}_{ p}}_{\text{$M$-times}}),g\big{)}=\Phi_{p}^{M}\big{(}(g(\mathfrak{x}_{1}),g( \mathfrak{x}_{2}),\ldots,g(\mathfrak{x}_{N}))\big{)}\approx(\mathcal{S}(g)( \mathfrak{x}_{1}),\mathcal{S}(g)(\mathfrak{x}_{2}),\ldots,\mathcal{S}(g)( \mathfrak{x}_{N})). \tag{32}\]

We have thus derived a model together with a family of model parameters which have the property that, under sufficient assumptions, the realization of the model for any of the model parameters in the family is approximating the solution operator at discrete space points.

## 4 Numerical simulations

In this section we test the ADANN methodology numerically in the case of three selected parametric PDE problems.
### One-dimensional reaction diffusion type equation

In this section we apply the ADANN methodology to a reaction diffusion type equation. Roughly speaking, we consider the semilinear heat PDE problem described in (15) when the nonlinearity is taken to be of the reaction diffusion type. Specifically, let \(\mathcal{I}=\{g\in C^{2}([0,1],\mathbb{R})\colon g(0)=g(1)=0\}\) and for every \(g\in\mathcal{I}\) consider the PDE

\[\big{(}\tfrac{\partial}{\partial t}u\big{)}(t,x)=\tfrac{1}{100}( \Delta_{x}u)(t,x)+\frac{1-u(t,x)}{1+(u(t,x))^{2}} \tag{33}\]

with Dirichlet boundary conditions \(u(t,0)=u(t,1)=0\) and initial value \(u(0,x)=g(x)\) for \((t,x)\in[0,\infty)\times[0,1]\). We want to approximate the mapping \(\mathcal{S}\colon\mathcal{I}\to C^{2}([0,1],\mathbb{R})\) from the initial condition of the PDE to the terminal value at time \(T=1\) (cf. Section 3.1). We apply the ADANN methodology as described in Section 2 to this problem with the following choices.

**Base model:** Roughly speaking, we use the base model for semilinear heat PDEs derived in Section 3 with \(N=35\) equidistant space discretization points and \(M=5\) timesteps adapted to the situation when the Laplace operator is multiplied with the factor \(\nicefrac{{1}}{{100}}\) as considered in (33). Concretely, this means that we let \(f\colon\mathbb{R}\to\mathbb{R}\) satisfy for all \(U\in\mathbb{R}\) that \(f(U)=\frac{1-U}{1+U^{2}}\) and we let \(\mathcal{B}\colon((\mathbb{R}^{35\times 35})^{5})^{5}\times\mathcal{I}\to \mathbb{R}^{35}\) satisfy for all \(W=((W_{m,i})_{i\in\{1,2,\ldots,5\}})_{m\in\{1,2,\ldots,5\}}\in((\mathbb{R}^{3 5\times 35})^{5})^{5}\), \(g\in\mathcal{I}\), \(U_{0},U_{1},\ldots,U_{5}\in\mathbb{R}^{35}\) with \(U_{0}=(g(\nicefrac{{1}}{{36}}),g(\nicefrac{{2}}{{36}}),\ldots,g(\nicefrac{{ 35}}{{36}}))\) and

\[\forall\,m\in\{1,2,\ldots,5\}\colon\qquad U_{m}=W_{m,1}U_{m-1}+W_ {m,2}f(U_{m-1})\\ +W_{m,3}f\big{(}W_{m,4}U_{m-1}+W_{m,5}f(U_{m-1})\big{)} \tag{34}\]

that

\[\mathcal{B}(W,g)=U_{5}. \tag{35}\]

To initialize this base model we use for some choices of \(p\in(0,\infty)^{2}\) the model parameters

\[(\underbrace{\mathbf{W}_{p},\ldots,\mathbf{W}_{p}}_{\text{5- times}})\in((\mathbb{R}^{35\times 35})^{5})^{5} \tag{36}\]

as proposed in (32) with the modification that the matrix \(A\) needs to be multiplied with the factor \(\nicefrac{{1}}{{100}}\).

**Difference model:** As difference model we use a classical feed forward ANN with layer dimensions \((35,50,150,35)\) and the GELU activation function. In every training run this difference model is initialized randomly with the standard Glorot uniform initialization.

**Training objective:** As a distribution for the initial value \(\mathfrak{I}\colon\Omega\to\mathcal{I}\) we use a sine expansion with decaying randomly distributed coefficients. Specifically, we have for all \(x\in[0,1]\) that

\[\mathfrak{I}(x)=\sum_{n=1}^{32}\frac{5Z_{n}\sin(\pi nx)}{n^{2}}, \tag{37}\]

where \(Z_{n}\colon\Omega\to\mathbb{R}\), \(n\in\{1,2,\ldots,16\}\), are independent standard normally distributed random variables. In addition, to measure the loss we use the seminorm \(\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left| \!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left| \!\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\! \left|\!\left|\!\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\!\left|\!\!\left|\! \left|\!\left|\!\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\! \left|\!\left|\!\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\! \left|\!\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\!\left|\!\!\left|\! \left|\!\left|\!\left|\!\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\! \left|\!\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\! \left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\! \left|\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\! \left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\!\left|\!\left|\!\! \left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\!\left|\!\left|\!\!\left|\! \left|\!\!\left|\!\left|\!\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\! \left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\!\left|\!\left|\! \left|\!\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\!\!\left|\!\!\left|\!\left|\!\! \left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\left|\!\left|\!\!\left|\!\left|\!\left|\!\!\left|\!\!\
\((1/2,1/2)\) with 300 timesteps for the temporal discretization (this LIRK scheme is sometimes referred to as the Crank-Nicolson explicit midpoint method (cf. Section A.3)), and only keeping the approximation of the terminal value at the space points \(\nicefrac{{1}}{{36}},\nicefrac{{2}}{{36}},\ldots,\nicefrac{{35}}{{36}}\). In summary, the training objective is to find a function \(\mathscr{N}\colon\mathcal{I}\to\mathbb{R}^{35}\) such that the loss

\[\mathbb{E}\big{[}\tfrac{1}{35}\left\|\mathscr{N}(\mathfrak{I})-\Psi(\mathfrak{ I})\right\|^{2}\big{]} \tag{39}\]

is as small as possible. If the loss is suitably small, we expect for suitable \(g\in\mathcal{I}\) that \(\mathscr{N}(g)\approx(\mathcal{S}(g)(\nicefrac{{1}}{{36}}),\mathcal{S}(g)( \nicefrac{{2}}{{36}}),\ldots,\mathcal{S}(g)(\nicefrac{{35}}{{36}}))\) is a good approximation of the terminal value at the space points \(\nicefrac{{1}}{{36}},\nicefrac{{2}}{{36}},\ldots,\nicefrac{{35}}{{36}}\).

**Runs over initializations and trainings:** We perform two numerical simulations with different loops over the initializations and training procedures of the base and difference model. In the first numerical simulation we use a grid-based approach to explore different possible initializations of the base model. Specifically, for every \(p\in\{\frac{1}{10},\frac{2}{10},\ldots,\frac{9}{10}\}\times\{\frac{3}{10}, \frac{4}{10},\ldots,\frac{9}{10}\}\) we

* initialize the base model with the specialized initialization described (32) corresponding to the parameter vector \(p\),
* train the base model with the ADAM optimization method with adaptive learning rates and batch sizes,
* randomly initialize the difference model with the standard Glorot uniform initialization, and
* train the difference model with the ADAM optimization method with adaptive learning rates and batch sizes.

We illustrate the results of this numerical simulation in Figure 1 and approximately summarize the performance in the rows 9-10 in Table 1. In the second numerical simulation we choose the parameters used to initialize the base model by randomly sampling them uniformly distributed on the set \([\frac{1}{10},\frac{9}{10}]\times[\frac{3}{10},\frac{9}{10}]\). We perform 50 runs and in each run we use an optimization based approach to decide whether

* to randomly sample a new parameter vector \(p\in[\frac{1}{10},\frac{9}{10}]\times[\frac{3}{10},\frac{9}{10}]\) and carry out (i)-(iv) above or
* to initialize a new instance of a difference model and train it to approximate the scaled error for one of the already trained base models.

We illustrate the results of this numerical simulation in Figure 2 and approximately summarize the performance in the row 11 in Table 1.

In order to compare the ADANN methodology with existing techniques from the literature we also show the results of existing deep learning methods and classical numerical methods when applied to the optimization problem in (39). Specifically, in row 1 in Table 1 we approximately summarize the performance of classical feedforward ANNs with the architecture \((35,100,220,150,35)\) and the GELU activation function and in row 2 in Table 1 we approximately summarize the performance of FNOs with 4 layers, width 32, and 13 Fourier modes. For both the ANNs and the FNOs we performed 3 training runs using the ADAM optimizer with adaptive learning rates and batch sizes to minimize the objective in (39) and present the resulting error of the best of the 3 runs. In rows 3-8 in Table 1 we approximately present the resulting errors when the approximations are computed by applying the classical Crank-Nicolson explicit midpoint method with finite difference spatial discretization with 35 space discretization steps and 15-20 timesteps to the PDE in (33). Note that this corresponds to the approximation described in (24) resp. (32) with \(M\in\{15,16,\ldots,20\}\), \(N=35\), and \(p=(1/2,1/2)\).

All the simulations in this subsection were run on a remote machine on [https://vast.ai](https://vast.ai) equipped with an NVIDIA GeForce RTX 3090 GPU with 24 GB RAM and an Intel(r)
Xeon(r) E5-2696 v2 CPU with 32 GB of total system RAM. All the ADANN, ANN, and FNO models were trained using \(2^{19}\) training samples, which were computed ahead of the training in 311 seconds.

### One-dimensional Sine-Gordon type equation

In this section we present numerical simulation results for the ADANN methodology applied to a one-dimensional Sine-Gordon type PDE with periodic boundary conditions. Specifically, let \(\mathcal{I}=\{g\in C^{2}([0,1],\mathbb{R})\colon g(0)=g(1)\text{ and }g^{\prime}(0)=g^{\prime}(1)\}\) and for every \(g\in\mathcal{I}\) consider the PDE

\[\big{(}\tfrac{\partial}{\partial t}u\big{)}(t,x)=\tfrac{1}{100}(\Delta_{x}u)(t,x)+\sin(u(t,x)) \tag{40}\]

with periodic boundary conditions \(u(t,0)=u(t,1)\) and \((\tfrac{\partial}{\partial x}u)(t,0)=(\tfrac{\partial}{\partial x}u)(t,1)\) and initial value \(u(0,x)=g(x)\) for \(t\in[0,\infty)\), \(x\in[0,1]\). We want to approximate the mapping \(\mathcal{S}\colon\mathcal{I}\to C^{2}([0,1],\mathbb{R})\)

Figure 1: Error plots in dependence of the initialization of the base model for the reaction diffusion type equation in (33). _Left_: Estimated \(L^{2}\)-errors of the base model prior to training. _Middle_: Estimated \(L^{2}\)-errors of the base model after training. _Right_: Estimated \(L^{2}\)-errors of the full ADANN model after training.

Figure 2: Estimated \(L^{2}\)-errors for the optimization based run for the reaction diffusion type equation in (33). All runs in between two vertical lines correspond to the same base model initialization.
from the initial condition of the PDE to the terminal value at time \(T=2\). We apply the ADANN methodology as described in Section 2 to this problem with the following choices.

**Base model:** Roughly speaking, we use the base model for semilinear heat PDEs derived in Section 3 with \(N=30\) equidistant space discretization points and \(M=15\) timesteps adapted to the situation of periodic boundary conditions and when the Laplace operator is multiplied with the factor \(\nicefrac{{1}}{{100}}\) as considered in (40). Concretely, this means that we let \(\mathscr{B}\colon((\mathbb{R}^{30\times 30})^{5})^{15}\times\mathcal{I} \to\mathbb{R}^{30}\) satisfy for all \(W=((W_{m,i})_{i\in\{1,2,\ldots,5\}})_{m\in\{1,2,\ldots,15\}}\in((\mathbb{R}^{3 0\times 30})^{5})^{15}\), \(g\in\mathcal{I}\), \(U_{0},U_{1},\ldots,U_{15}\in\mathbb{R}^{30}\) with \(U_{0}=(g(\nicefrac{{0}}{{30}}),g(\nicefrac{{1}}{{30}}),\ldots,g(\nicefrac{{ 29}}{{30}}))\) and

\[\forall\,m\in\{1,2,\ldots,15\}\colon\qquad U_{m}=W_{m,1}U_{m-1}+W _{m,2} \text{sin}(U_{m-1})\\ +W_{m,3}\text{sin}\big{(}W_{m,4}U_{m-1}+W_{m,5}\text{sin}(U_{m-1 })\big{)} \tag{41}\]

that

\[\mathscr{B}(W,g)=U_{15}. \tag{42}\]

To initialize this base model we use for some choices of \(p\in(0,\infty)^{2}\) the model parameters

\[(\underbrace{\mathbf{W}_{p},\ldots,\mathbf{W}_{p}}_{\text{15-times}})\in(( \mathbb{R}^{30\times 30})^{5})^{15} \tag{43}\]

as proposed in (32) with the modification that the matrix \(A\) is multiplied with the factor \(\nicefrac{{1}}{{100}}\) and adjusted to periodic boundary conditions.

**Difference model:** As difference model we use an FNO with 3 layers, width 20, and 16 Fourier nodes. In every training run this difference model is initialized randomly as proposed in [30].

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline \multirow{2}{*}{Method} & Estimated & Estimated & Number & Training & Time for 4096 \\  & \(L^{1}\)-error & \(L^{2}\)-error & of trainable parameters & time (in s) & evaluations (in s) \\ \hline ANN & 0.003430 & 0.003950 & 64255 & 367 & 0.007 \\ \hline FNO & 0.002411 & 0.003175 & 61921 & 593 & 0.050 \\ \hline Classical (5 timesteps) & 0.001683 & 0.001785 & 30625 & 0 & 0.023 \\ \hline Classical (6 timesteps) & 0.001301 & 0.001360 & 36750 & 0 & 0.025 \\ \hline Classical (7 timesteps) & 0.001123 & 0.001157 & 42875 & 0 & 0.033 \\ \hline Classical (8 timesteps) & 0.001025 & 0.001047 & 49000 & 0 & 0.040 \\ \hline Classical (9 timesteps) & 0.000964 & 0.000982 & 55125 & 0 & 0.035 \\ \hline Classical (10 timesteps) & 0.000925 & 0.000940 & 61250 & 0 & 0.044 \\ \hline Best trained base model (grid) & 0.000174 & 0.000189 & 30625 & 12270 & 0.020 \\ \hline Best trained full ADANN (grid) & 0.000015 & 0.000019 & 45360 & 21147 & 0.027 \\ \hline Best trained full ADANN (opt) & 0.000062 & 0.000065 & 45360 & 13393 & 0.026 \\ \hline \end{tabular}
\end{table}
Table 1: Numerical simulations for the reaction diffusion type PDE in (33)

Figure 3: Example approximation plots for a randomly chosen initial value for the reaction diffusion type equation in (33). _Left_: Best ANN approximation. _Middle_: Best FNO approximation. _Right_: Best full ADANN approximation of the optimization based run.
**Training objective:** As a distribution for the initial value \(\mathcal{I}\colon\Omega\to\mathcal{I}\) we use a Fourier expansion with decaying randomly distributed coefficients. Specifically, we have for all \(x\in[0,1]\) that

\[\mathfrak{I}(x)=2Z_{0}+\sum_{n=1}^{16}\frac{2Z_{n}\sin(2\pi nx)}{n^{2}}+\frac{2Z _{n}\cos(2\pi nx)}{n^{2}}, \tag{44}\]

where \(Z_{n}\colon\Omega\to\mathbb{R}\), \(n\in\{-32,-31,\ldots,32\}\), are independent standard normally distributed random variables. In addition, to measure the loss we use the seminorm \(\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left| \!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\left|\! \left|\!\left|\!\left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\! \left|\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left| \!\!\left|\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\! \left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\!\left|\!\!\left|\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left| \!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left| \!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\! \!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\!\left|\!\!\left|\!\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\! \left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\left|\!\!\! \left|\!\left|\!\!\left|\!\!\left|\!\!\!\left|\!\!\left|\!\!\!\left|\!\!\left|\!\!\! \left|\!\!
We illustrate the results of this numerical simulation in Figure 5 and approximately summarize the performance in the row 12 in Table 2.

In order to compare the ADANN methodology with existing techniques from the literature we also show the results of deep learning methods and classical numerical methods when applied to the optimization problem in (46). Specifically, in row 1 in Table 2 we approximately summarize the performance of classical feedforward ANNs with the architecture \((30,100,300,160,30)\) and GELU activation function and in row 2 in Table 2 we approximately summarize the performance of FNOs with 4 layers, width 34, and 13 Fourier modes. For both the ANNs and the FNOs we performed 3 training runs using the ADAM optimizer with adaptive learning rates and batch sizes to minimize the objective in (46) and present the resulting error of the best of the 3 runs. In rows 3-8 in Table 2 we approximately present the resulting errors when the approximations are computed by applying the classical Crank-Nicolson explicit midpoint method with finite difference spatial discretization with 30 space discretization steps and 15-20 timesteps to the PDE in (40). Note that this corresponds to the approximation described in (24) resp. (32) with \(M\in\{15,16,\ldots,20\}\), \(N=30\), and \(p=(1/2,1/2)\).

All the simulations in this subsection were run on a remote machine on [https://vast.ai](https://vast.ai) equipped with an NVIDIA GeForce RTX 3090 GPU with 24 GB RAM and an Intel(r) Core(tm) i3-10100 CPU with 32 GB of total system RAM. All the ADANN, ANN, and FNO models were trained using \(2^{19}\) training samples, which were computed ahead of the training in 819 seconds.
### Two-dimensional semilinear heat equation

In this section we present numerical simulation results for the ADANN methodology applied to a two-dimensional semilinear heat PDE with periodic boundary conditions. Specifically, let

\[\mathcal{I}=\left\{g\in C^{2}([0,1]^{2},\mathbb{R})\colon\left(\begin{array}{c}g (z,0)=g(z,1),\\ g(0,z)=g(1,z),\\ g^{\prime}(z,0)=g^{\prime}(z,1),\\ g^{\prime}(0,z)=g^{\prime}(1,z)\end{array}\right)\right\} \tag{47}\]

and for every \(g\in\mathcal{I}\) consider the PDE

\[\big{(}\tfrac{\partial}{\partial t}u\big{)}(t,x)=\tfrac{1}{100}(\Delta_{x}u)( t,x)+\big{(}1+(u(t,x))^{2}\big{)}^{1/2} \tag{48}\]

with periodic boundary conditions \(u(t,(z,0))=u(t,(z,1))\), \(u(t,(0,z))=u(t,(1,z))\), \(\big{(}\tfrac{\partial}{\partial x}u\big{)}(t,(z,0))\)\(=(\tfrac{\partial}{\partial x}u)(t,(z,1))\), and \(\big{(}\tfrac{\partial}{\partial x}u\big{)}(t,(0,z))=(\tfrac{\partial}{ \partial x}u)(t,(1,z))\) and initial value \(u(0,x)=g(x)\) for \(t\in[0,\infty)\), \(z\in[0,1]\), \(x\in[0,1]^{2}\). We want to approximate the mapping \(\mathcal{S}\colon\mathcal{I}\to C^{2}([0,1]^{2},\mathbb{R})\) from the initial condition of the PDE to the terminal value at time \(T=2\). We apply the ADANN methodology as described in Section 2 to this problem with the following choices.

Figure 5: Estimated \(L^{2}\)-errors for the optimization based run for the Sine-Gordon type equation in (40)

Figure 6: Example approximation plots for a randomly chosen initial value for the reaction diffusion type equation in (33). _Left_: Best ANN approximation. _Middle_: Best FNO approximation. _Right_: Best full ADANN approximation of the optimization based run.
**Base model:** Roughly speaking, we adapt the base model for one-dimensional semilinear heat PDEs derived in Section 3 to the two-dimensional situation with periodic boundary conditions considered in (48). For the finite difference spatial discretization we employ a \(40\times 40\) grid of equidistant points and for the second order LIRK temporal discretization we use 2 timesteps. Since we consider \(40^{2}=1600\) space discretization points the trainable base model parameters are matrices of the form \(W=((W_{m,i})_{i\in\{1,2,\ldots,5\}})_{m\in\{1,2\}}\in((\mathbb{R}^{1600\times 1 600})^{5})^{2}\) and they are employed to compute the output of the base model analogous to the description in (30)-(31). To initialize this base model we use initializations corresponding to parameters \(p\in(0,\infty)^{2}\) of the LIRK scheme as described in (32), with the modification that the matrix \(A\) represents the two-dimensional Laplace operator discretized with finite differences on a \(40\times 40\)-grid, multiplied by the factor \(\nicefrac{{1}}{{100}}\) and adjusted to periodic boundary conditions.

**Difference model:** As difference model we use an FNO with 4 layers, width 30, and 15 Fourier modes. In every training run this difference model is initialized randomly as proposed in [30].

**Training objective:** We choose the initial value \(\mathfrak{I}\colon\Omega\to\mathcal{I}\) to be \(\mathcal{N}(0,4(2I-\Delta)^{-2})\)-distributed, where \(\Delta\) denotes the Laplace operator with periodic boundary conditions on \(L^{2}([0,1]^{2})\). In our numerical simulations we approximate this distribution on an \(80\times 80\) equidistant grid. Specifically, let \(\mathbf{N}=80\), let \(X\colon\Omega\to\mathbb{R}^{\mathbf{N}\times\mathbf{N}}\) be a random matrix with independent \(\mathcal{N}(0,\mathbf{N}^{2})\)-distributed entries, let \(\Delta_{\mathbf{N}}\colon\mathbb{R}^{\mathbf{N}\times\mathbf{N}}\to\mathbb{R} ^{\mathbf{N}\times\mathbf{N}}\) be the discrete Laplace operator with periodic boundary conditions, let \(I_{\mathbf{N}}\colon\mathbb{R}^{\mathbf{N}\times\mathbf{N}}\to\mathbb{R}^{ \mathbf{N}\times\mathbf{N}}\) be the identity operator. We employ the approximation

\[\big{(}\mathfrak{I}\big{(}\tfrac{i}{\mathbf{N}},\tfrac{j}{\mathbf{N}}\big{)} \big{)}_{(i,j)\in\{0,1,\ldots,\mathbf{N}-1\}^{2}}\approx 2(2I_{\mathbf{N}}- \Delta_{\mathbf{N}})^{-1}X. \tag{49}\]

In addition, to measure the loss we use the seminorm \(\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|} \!\big{|}\!\big{|}\!\big{|}\!\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\! \big{|}\!\big{|}\!\big{|}\!\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!\big{|}\!
and GELU activation function and in row 2 in Table 3 we approximately summarize the performance of FNOs with 5 layers, width 80, and 20 Fourier modes. For both the ANNs and the FNOs we used the ADAM optimizer with adaptive learning rates. In row 3 in Table 3 we approximately present the resulting errors when the approximations are computed by applying the classical Crank-Nicolson explicit midpoint method with finite difference spatial discretization with \(40\times 40\) space discretization steps and 2 timesteps to the PDE in (48).

All the simulations in this subsection were run on a remote machine on [https://vast.ai](https://vast.ai) equipped with an NVIDIA GeForce RTX 3090 GPU with 24 GB RAM and an AMD EPYC(tm) 7502 CPU with 32 GB of total system RAM. All the ADANN, ANN, and FNO models were trained using \(2^{14}\) training samples, which were computed ahead of the training in 2065 seconds.

## Appendix A Second order linearly implicit Runge-Kutta methods

In this section we present a formal derivation of a well-known family of second order LIRK methods for semilinear ODEs, which are used to construct a base model and a corresponding family of initialization parameters in Section 3 (cf., e.g., Deuflhard & Bornemann [9, Section 6.4] and Hochbruck & Ostermann [21]). We will work in the following setting. Let \(d\in\mathbb{N}\), \(A\in\mathbb{R}^{d\times d}\), \(f\in C(\mathbb{R}^{d},\mathbb{R}^{d})\) and consider the ODE

\[\dot{u}(t)=Au(t)+f(u(t)) \tag{51}\]

for \(t\in(0,\infty)\).

### Order conditions for general LIRK methods

We first introduce the one step increment function of general LIRK methods for the ODE in (51). Specifically, let \(s\in\mathbb{N}\), \(\alpha=(\alpha_{i,j})_{(i,j)\in\{1,2,\ldots,s\}^{2}}\in\mathbb{R}^{s\times s}\), \(\beta=(\beta_{i,j})_{(i,j)\in\{1,2,\ldots,s\}^{2}}\in\mathbb{R}^{s\times s}\)

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|} \hline \multirow{2}{*}{Method} & Estimated & Estimated & Number & Precomputation & Time for 512 \\  & \(L^{1}\)-error & \(L^{2}\)-error & parameters & time (in s) & evaluations (in s) \\ \hline ANN & 0.220621 & 0.240503 & 30731200 & 1472 & 0.76 \\ \hline FNO & 0.013329 & 0.013707 & 25643217 & 2085 & 0.82 \\ \hline Classical (2 timesteps) & 0.666985 & 0.939649 & 25600000 & 0 & 0.79 \\ \hline Best trained base model & 0.012739 & 0.014994 & 25600000 & 10768 & 0.79 \\ \hline Best trained full ADANN & 0.001291 & 0.001543 & 27227937 & 23359 & 1.01 \\ \hline \end{tabular}
\end{table}
Table 3: Numerical simulations for the two-dimensional semilinear heat PDE in (48)

Figure 7: Error plots in dependence of the initialization of the base model for the two-dimensional semilinear heat PDE in (48). _Left_: Estimated \(L^{2}\)-errors of the base model prior to training. _Middle_: Estimated \(L^{2}\)-errors of the base model after training. _Right_: Estimated \(L^{2}\)-errors of the full ADANN model after training.
\(b=(b_{i})_{i\in\{1,2,\ldots,s\}}\in\mathbb{R}^{s}\) and let \(\Phi=(\Phi^{h}(u))_{(h,u)\in[0,\infty)\times\mathbb{R}^{d}}\colon[0,\varepsilon] \times\mathbb{R}^{d}\to\mathbb{R}^{d}\) satisfy for all \(h\in[0,\infty)\), \(U,k_{1},k_{2},\ldots,k_{s}\in\mathbb{R}^{d}\) with

\[\forall\,i\in\{1,2,\ldots,s\}\colon\quad k_{i}=A(U+h\sum_{j=1}^{i}\beta_{i,j} k_{j})+f(U+h\sum_{j=1}^{i-1}\alpha_{i,j}k_{j}) \tag{52}\]

that

\[\Phi^{h}(U)=U+h\sum_{i=1}^{s}b_{j}k_{j}. \tag{53}\]

We refer to the number \(s\) as the number of stages of the LIRK method, we refer to \(\alpha\) as the nonlinear LIRK parameters, we refer to \(\beta\) as the linear LIRK parameters, we refer to \(b\) as the LIRK integration weights, and we refer to \(k_{1},k_{2},\ldots,k_{s}\) as the LIRK stages. Although the LIRK stages are defined implicitly in (52), under suitable conditions they can be computed explicitly. Specifically, under suitable conditions, we have for all \(h\in[0,\infty)\), \(U,k_{1},k_{2},\ldots,k_{s}\in\mathbb{R}^{d}\) with

\[\forall\,i\in\{1,2,\ldots,s\}\colon\quad k_{i}=(I_{d}-h\beta_{i,i}A)^{-1} \Big{(}A(U+h\sum_{j=1}^{i-1}\beta_{i,j}k_{j})+f(U+h\sum_{j=1}^{i-1}\alpha_{i,j }k_{j})\Big{)} \tag{54}\]

that

\[\Phi^{h}(U)=U+h\sum_{i=1}^{s}b_{j}k_{j}. \tag{55}\]

Order conditions for the one step method \(\Phi\) are obtained by formally setting the Taylor expansion of \(\Phi\) equal to the Taylor expansion of the solution of the ODE in (51) for a fixed initial value \(U\in\mathbb{R}^{d}\) up to terms of a certain order. The resulting order conditions for a second order scheme are given by

\[\sum_{i=1}^{s}b_{i}=1\qquad\text{and}\qquad\sum_{i=1}^{s}b_{i}C_{i}=\sum_{i=1} ^{s}b_{i}c_{i}=\frac{1}{2}, \tag{56}\]

where \((C_{i})_{i\in\{1,2,\ldots,s\}},(c_{i})_{i\in\{1,2,\ldots,s\}}\subseteq\mathbb{R}\) satisfy for all \(i\in\{1,2,\ldots,s\}\) that

\[C_{i}=\sum_{j=1}^{i}\beta_{i,j}\qquad\text{and}\qquad c_{i}=\sum_{j=1}^{i-1} \alpha_{i,j}. \tag{57}\]

Under suitable regularity on the nonlinearity \(f\) the conditions in (56) ensure that the ODE integration scheme defined through the one step increment function \(\Phi\) will have global convergence order 2.

### A family of 2 stage linearly implicit Runge-Kutta methods of order 2

In this section we solve the order conditions in (56) in the case of \(s=2\) stages and under the assumption that \(\beta_{1,1}=\beta_{2,2}\). For this let \(p_{1},p_{2}\in(0,\infty)\) and assume that

\[\alpha_{1,2}=p_{1}\qquad\text{and}\qquad\beta_{1,1}=\beta_{2,2}=p_{2}. \tag{58}\]

This and (56) imply that

\[b_{1}=1-\tfrac{1}{2p_{1}},\qquad b_{2}=\tfrac{1}{2p_{1}},\qquad\text{and} \qquad\beta_{1,2}=2p_{1}(\tfrac{1}{2}-p_{2}). \tag{59}\]

This and (54) in turn imply, under suitable conditions, that for all \(h\in[0,\infty)\), \(U,k_{1},k_{2}\in\mathbb{R}^{d}\) with

\[k_{1}=(I_{d}-hp_{2}A)^{-1}(AU+f(U))\qquad\text{and} \tag{60}\]

\[k_{2}=(I_{d}-hp_{2}A)^{-1}\big{(}A(U+h2p_{1}(\tfrac{1}{2}-p_{2})k_{1})+f(U+hp_{ 1}k_{1})\big{)} \tag{61}\]

it holds that

\[\Phi^{h}(U)=U+h\big{[}(1-\tfrac{1}{2p_{1}})k_{1}+(\tfrac{1}{2p_{1}})k_{2}\big{]}. \tag{62}\]

We have thus derived a family of LIRK methods of order two, which is parametrized by two parameters \(p_{1}\) and \(p_{2}\). We use this family in Section 3.2.2.
### The special case of the Crank-Nicolson explicit Euler method

The scheme in (60)-(62) includes as a special case the well-known Crank-Nicolson explicit midpoint scheme. Specifically, note that in the special case where \(p_{1}=p_{2}=\frac{1}{2}\) we have for all \(h\in[0,\infty)\), \(U\in\mathbb{R}^{d}\) that

\[\Phi^{h}(U)=(I_{d}-\tfrac{h}{2}A)^{-1}\big{(}(I_{d}+\tfrac{h}{2}A)U+hf\big{(}( I_{d}-\tfrac{h}{2}A)^{-1}(U+\tfrac{h}{2}f(u))\big{)}\big{)}. \tag{63}\]

## Acknowledgments

This work has been partially funded by the National Science Foundation of China (NSFC) under grant number 12250610192. This work has been partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC 2044-390685587, Mathematics Munster: Dynamics-Geometry-Structure.
