FMint: Bridging Human Designed and Data Pretrained Models for Differential Equation Foundation Model

## Abstract

Abstract The fast simulation of dynamical systems is a key challenge in many scientific and engineering applications, such as weather forecasting, disease control, and drug discovery. With the recent success of deep learning, there is increasing interest in using neural networks to solve differential equations in a data-driven manner. However, existing methods are either limited to specific types of differential equations or require large amounts of data for training. This restricts their practicality in many real-world applications, where data is often scarce or expensive to obtain. To address this, we propose a novel multi-modal foundation model, named FMint ( F oundation M odel based on In i t ialization), to bridge the gap between human-designed and data-driven models for the fast simulation of dynamical systems. Built on a decoder-only transformer architecture with in-context learning, FMint utilizes both numerical and textual data to learn a universal error correction scheme for dynamical systems, using prompted sequences of coarse solutions from traditional solvers. The model is pre-trained on a corpus of 40K ODEs, and we perform extensive experiments on challenging ODEs that exhibit chaotic behavior and of high dimensionality. Our results demonstrate the effectiveness of the proposed model in terms of both accuracy and efficiency compared to classical numerical solvers, highlighting FMint’s potential as a general-purpose solver for dynamical systems. Our approach achieves an accuracy improvement of 1 to 2 orders of magnitude over state-of-the-art dynamical system simulators, and delivers a 5X speedup compared to traditional numerical algorithms. The code for FMint is available at https://github.com/margotyjx/FMint .

## 1 Introduction

Dynamical systems characterize the evolution of physical states over time. They are fundamental in describing the change of physical states across a wide range of disciplines, including physics [, , ], chemistry [, ], engineering [, , ], and finance [, ].
Typically, these systems are formulated as systems of ordinary differential equations (ODEs):

$$ $\frac{d{\bm{u}}(t)}{dt}=f[{\bm{u}}(t)],\quad{\bm{u}}(0)={\bm{c}}_{0},$ (1) $$

where ${\bm{c}}_{0}$ denotes the initial condition of the system.
To solve these systems numerically, one usually employs a human-designed numerical integration algorithm such as the Euler method or Runge-Kutta methods.
These methods can be adapted easily to solve different types of ODEs that share the same format with guaranteed accuracy.
The implementation is given as

$$ ${\bm{u}}_{n+1}={\bm{u}}_{n}+S(f,{\bm{u}}_{n},\Delta t_{n}),\quad{\bm{u}}_{0}=c _{0},\quad n=0,1,\cdots,$ (2) $$

where $S$ represents the numerical integration scheme, $\Delta t_{n}$ is the step size at the $n$-th time step, and ${\bm{u}}_{n}\in\mathbb{R}^{n}$ is the approximated solution at the cumulative time $\sum_{i=0}^{n}\Delta t_{i}$.

One obstacle of these human-designed algorithm is the trade-off between accuracy and efficiency. This makes the large-scale simulation using these numerical schemes impossible.
In fact, in many real-world scenarios, high-volume simulation that produces forecasts on a set of initial conditions simultaneously plays a significant role in various applications.
For example, simulations of virus propagation during an epidemic given different circumstances are necessary for formulating health regulations [].
In these scenarios, it is practical to standardize the time step $\Delta t:=\Delta t_{1}=\Delta t_{2}=\cdots$ across simulations, facilitating batch processing.
Yet, this standardization introduces a trade-off between accuracy and efficiency: a larger time step speeds up the simulation at the cost of increased simulation error, while a smaller time step reduces the error but slows down the simulation.
Therefore, the long runtime makes these traditional algorithms unsuitable for wide range simulations in many practical situations.

Recently, deep learning methods have demonstrated remarkable success across various scientific domains, including solving partial differential equations (PDEs) [, ], learning operators [], and addressing inverse problems [, , ].
However, they typically underperform in data-scarce environments and may lack essential domain knowledge. Therefore, we ask an important question:

Is it possible to combine the best of both worlds, leveraging the efficiency of human-designed algorithms and the accuracy of data-driven methods?

To address this, we adopt a multi-modal approach to develop a foundation model, FMint (Foundation Model based on Initialization), a pre-trained foundation model designed to speed up large-scale simulations of dynamical systems with high accuracy via error correction.
We integrate human expertise i.e., traditional ODE solvers into modern data-driven methods and adapt the idea of in-context learning to obtain refined solutions based on the initialization of coarse solutions that are computed using human-designed integration method for various differential equations.

In the scientific computing field, most models use solely numerical data as input. However, we ask: Can we integrate other data modalities in scientific computing to enrich the information available? Various characteristic behaviors of dynamical systems can be generalized by textual descriptions, providing more context to the systems. For example, for the Schrödinger equation $i\hbar\frac{\partial}{\partial t}\Psi(x,t)=-\frac{\hbar^{2}}{2m}\nabla^{2}\Psi
(x,t)+V(x)\Psi(x,t),$
where $\Psi(x,t)$ is the wave function, $\hbar$
is the Planck constant, and $V(x)$ is the potential energy, the possible energy levels of the electron, its spatial distribution, and its interaction with the environment has been analyzed.
We hence further incorporate textual modalities in our model that provide guidance to the numerical data.

Modal 1: Numerical Data
For the numerical modal data, we input the pairs of coarse solutions from traditional numerical solver and
the corresponding error correction term to the model. Based on in-context learning, given prompted sequences of such pairs, the model is trained to predict the error correction term for the query coarse solution to achieve high accuracy. Such a paradigm allows the model to generalize to various downstream tasks with a short prompt of examples, making it efficient and accurate in a data-scarce environment.

Modal 2: Textual Data (optional) FMint can also take textual information as supplemental input to the model. It includes the descriptive information about the equations such as the mathematical expression, or physical meaning of the parameters, etc. We would like to mention here that the textual data is only served as an optional input and is not necessary if not available.

Through extensive experiments involving ODEs with qualitatively different behaviors (periodic, chaotic, etc),
we demonstrate the effectiveness of FMint in terms of both accuracy and efficiency over classical numerical solvers and deep learning based methods.

We summarize our contributions as follows:

(1) To the best of our knowledge, FMint is the first multi-modal foundation model that synthesizes human-designed algorithms and deep learning framework. Back-boned on the decoder-only transformer with in-context learning scheme, FMint achieves competitive accuracy and efficiency in simulating high-dimensional ODEs with qualitatively different behaviors.

(2) We obtained 10 to 100 times higher accuracy than state-of-the-art dynamical system simulators, and 5X speedup compared to traditional numerical algorithms with remarkable generalization ability.

## 2 Methodology

### 2.1 Multi-Modal Data Preparation

FMint is a multi-modal foundation modal that bridges human-designed algorithms and data-driven methods. It takes two types of modalities: 1) Numerical Data and, 2) Textual data. The details of the data preparation is summarized below.

Modality 1: Numerical Data
In solving 1 for large-scale simulations, we consider selecting a numerical integration scheme that utilizes a large time step size. This can be written in stride $k\in\{1,2,...\}$ and step size $\Delta t$ that results in desired accuracy, denoted as $k\Delta t$. For illustrative purposes, we consider the Euler method, which yields the following numerical simulation scheme:

$$ $\hat{{\bm{u}}}(t+k\Delta t)=\hat{{\bm{u}}}(t)+f[\hat{{\bm{u}}}(t)]\cdot k \Delta t.$ (3) $$

However, solving the dynamical system 1 with numerical scheme 3 and large step size $k\Delta t$ unavoidably causes large simulation errors. From the Taylor expansion

$$ ${\bm{u}}(t+k\Delta t)=\underbrace{{\bm{u}}(t)+f[{\bm{u}}(t)]\cdot k\Delta t}_{ \text{For Euler method }}+\sum_{n=2}^{\infty}\underbrace{\frac{1}{n!}\frac{ \mathrm{d}^{n}}{\mathrm{~{}d}t^{n}}{\bm{u}}(t)\cdot[k\Delta t]^{n}}_{\text{ err}_{n}(k,\Delta t,{\bm{u}}(t))},$ (4) $$

we see that the error term $\sum_{n=2}^{\infty}\text{err}_{n}(k,\Delta t,{\bm{u}}(t))$ is non-negligible and this limits the fast simulation of real-world dynamical systems. We therefore consider building a corrector foundation model that approximates $\sum_{n=2}^{\infty}\text{err}_{n}$ for various dynamical systems.
We call solutions obtained by vanilla numerical integration schemes 3 with time step $k\Delta t$ as “coarse solutions”. With coarse solutions as an initialization, our goal is to produce highly accurate solution with fast inference time on a diverse set of dynamical systems, i.e.,

$$ $\hat{{\bm{u}}}_{k(n+1)}=\hat{{\bm{u}}}_{kn}+S\left(f,\hat{{\bm{u}}}_{kn},k \Delta t\right)+\operatorname{FMint}\left(\hat{{\bm{u}}}_{kn};\Theta\right), \quad\hat{{\bm{u}}}_{0}={\bm{c}}_{0},\quad n=0,1,\cdots,$ (5) $$

where $\Theta$ represents all the model parameters. We designed our model using a decoder-only transformer backbone [].
The model is trained to perform in-context learning such that it predicts the error correction term in examples based on previous demonstrations. The training is done in a similar manner to the next-token-prediction scheme.

Training Data Preparation. We construct FMint to learn the corrector from multiple demos from the same ODE system, each consists of coarse solutions and their corresponding correction term.
To prepare for the training data, for $i$-th ODE equation, we first simulate using fine step size $\Delta t$ and obtain ODE $\{{\bm{u}}^{i}_{j}\}_{j=1}^{kn}$ where ${\bm{u}}^{i}_{j}$ represents the fine-grained solution for $i$-th ODE system at time step $j\Delta t$.
Then using coarse step size $k\Delta t$, we generate ODE results $\{\hat{{\bm{u}}}^{i}_{kj}\}_{j=1}^{n}$ where we denote $\hat{{\bm{u}}}^{i}_{kj}$ the coarse solution for $i$-th ODE equation at time step $kj\Delta t$ with predefined stride $k$.
The corresponding error correction term for each coarse solutions are computed from the difference

$$ $\text{err}_{\hat{{\bm{u}}}_{kj}}={\bm{u}}_{k(j+1)}-\hat{{\bm{u}}}_{kj}-S\left( f,\hat{{\bm{u}}}_{kj},k\Delta t\right).$ (6) $$

One pair of coarse solutions $\hat{{\bm{u}}}^{i}=\{\hat{{\bm{u}}}^{i}_{kj}\}_{j=1}^{n}$ and error term $\text{err}^{i}=\{\text{err}_{\hat{{\bm{u}}}^{i}_{kj}}\}_{j=1}^{n}$ composes one demo. The model takes a collection of demos of size $d$, a query data sequence $\hat{{\bm{u}}}^{t}$ and outputs an error correction term $\text{err}^{t}$ for the query data

$$ $\left\{\{\hat{{\bm{u}}}^{1},\text{err}^{1}\},\{\hat{{\bm{u}}}^{2},\text{err}^{ 2}\},\ldots,\{\hat{{\bm{u}}}^{d},\text{err}^{d}\},\hat{{\bm{u}}}^{t}\right\} \rightarrow\text{err}^{t}.$ (7) $$

All the demo information will be tokenized before passed into the FMint model. We describe the tokenization in detail in Subsection 2.2.

Pretraining Data.
We pretrain the FMint model with four types of ODEs that exhibit qualitatively different characteristics to cover a wide range of diversity in data:

(1) Lorenz model. The Lorenz system is a 3D system, well-known for its chaotic behavior, originally developed to model atmospheric convection. It is governed by three coupled, nonlinear differential equations.
(2) Damped oscillator. A damped oscillator is a 2nd order system that describes a system where the motion of the oscillator is subject to a force that reduces its amplitude over time. This leads to a decay in oscillation.
(3) Van der Pol oscillator. The Van der Pol oscillator is a 2nd-order nonlinear oscillator. It can sustain oscillations indefinitely, known as limit cycles.
(4) Lotka-Volterra. The Lotka-Volterra system is a 2D system that describes the nonlinear interaction between two species: the prey and the predator.

Modality 2: Textual Data.
For the optional multi-modal training, we produced 30 descriptions per ODE as supplemental textual data for training and testing. These data contains the mathematical expression, exact parameters used for each ODEs, or behaviors under different parameter ranges.
These data are generated with the help of GPT-4. We provide the details of the textual data generation in Appendix A.5 and provide some examples listed in Appendix A.6.

### 2.2 Model Design and Modality Fusion

In this subsection, we describe the model architecture and the fusion of modalities.

**Table 1: Input tokens for a 2D ODE demo.**
|  | Coarse solution | Error term | Query |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| key | 0 | $\ldots$ | $t_{n}$ | 0 | $\ldots$ | $t_{n}$ | 0 | $\ldots$ | $t_{q}$ |
| value | $\hat{u}(0)$ | $\ldots$ | $\hat{u}(t_{n})$ | $\text{err}_{\hat{u}}(0)$ | $\ldots$ | $\text{err}_{\hat{u}}(t_{n})$ | $\hat{u}^{q}(0)$ | $\ldots$ | $\hat{u}^{q}(t_{q})$ |
|  | $\hat{v}(0)$ | $\ldots$ | $\hat{v}(t_{n})$ | $\text{err}_{\hat{v}}(0)$ | $\ldots$ | $\text{err}_{\hat{v}}(t_{n})$ | $\hat{v}^{q}(0)$ | $\ldots$ | $\hat{v}^{q}(t_{q})$ |

Tokenization FMint processes two modalities of data: numerical and textual. Like other language models, FMint requires tokenization of the input data before passing them into the transformer. The numerical data can be classified into three categories: coarse solutions, error corrections, and query solutions.
Query solutions consists of the coarse solution of the query trajectory.
As previously discussed, FMint employs in-context learning using pairs of coarse solutions and their associated error corrections to predict the error corrections for the query locations. As shown in Table 1, each column represents a token, with each token belonging to one of the three categories. To handle these three categories, we utilize a shared embedding layer that maps them into embedding vectors. Additionally, a learnable positional embedding is appended to the embedding vectors of each category. Notably, the positional embeddings are consistent within each data category but differ across the categories. For tokenizing textual data, we employ a separate embedding layer, and the resulting embedding vectors are appended with positional embeddings.

Model architecture.
Similar to language models, FMint is a decoder-only transformer model, where the key design aspect lies in the appropriate masking of the attention mechanism. The mask must satisfy the following requirements: (1) When predicting query locations, the model should be invariant to the order of queries, allowing predictions to be made independently and in parallel. (2) When predicting queries in the current example, the query tokens must have access to the tokens of coarse solutions and error corrections from both previous and the current examples. To effectively manage these constraints, we use a specialized masking technique introduced by  [].
In particular, the mask is designed such that when predicting the current queried error correction, the model "sees" the prompt, all previous coarse solutions with error terms and queries with the queried error terms, but not the current queried error term. After passing through multiple attention blocks, FMint is connected to a head layer (multi-layer perceptron), which outputs the error correction predictions for the query locations of the system.
The architecture of FMint is shown in Figure 1.

Figure: Figure 1: Work flow of FMint model.
Refer to caption: x1.png

Training and inference.
We pretrain FMint on 40K ODE simulation data using an in-context learning scheme. During training, ODE trajectories with the same parameters but different initial conditions are passed to the model as demo pairs. The query tokens share the same time variables as the error correction tokens, but their values are set to 0. The training loss function is the mean squared error (MSE) 8, which minimizes the difference between the predicted error corrections for the query data and the ground truth error. Thanks to the mask design, the model outputs predictions in a language generation-like manner, using information from both the current and previous examples.

$$ $\text{MSE}=\frac{1}{N}\sum_{i=1}^{N}\|\tilde{\text{err}}^{i}-\text{err}^{i}\|_ {2}$ (8) $$

During inference, we provide FMint with a few demo pairs of coarse solutions and corresponding error corrections. Once given the query locations, FMint accurately predicts the error corrections, achieving high accuracy even from coarse, low-accuracy simulation data. The generation of query locations is highly flexible and order-invariant, allowing users to perform data interpolation. Notably, for relatively simple dynamical systems, FMint can perform zero-shot or few-shot learning. For more complex systems, FMint still achieves high accuracy with minimal fine-tuning, as shown in the numerical results section 3.

## 3 Experimental Results

In this section, we conduct extensive experiments to evaluate the effectiveness of FMint and compare its results with SOTA baslines [, , ] on a wide range of ODEs from 1D to 3D. We aim to answer the following research questions:

RQ1. How does FMint perform compared to baseline models in terms of accuracy and efficiency.

RQ2. How does FMint adapt to ODEs’ different behavior under various data generation coefficients.

RQ3. Whether the optional textual data improves the performance of FMint.

### 3.1 Basic Set-Up

Pretraining data preparation. The pretraining data consists of 40K ODEs that are commonly observed in important applications in engineering and science: (1) Lorenz model, (2) Damped oscillator, (3) Van der Pol oscillator, and (4) Lotka-Volterra (LV) dynamics. For each ODE system, we created 1000 variations with different parameters and for each variation, we produce 100 trajectories with different initial conditions. Consequently, our dataset comprises trajectories of 100,000 ODEs for each dynamical system, differentiated by varying coefficients and initial conditions.

We have included more details on these ODE systems in Appendix A.1.
For all four ODE systems, the coefficients’ range are provided in Table 8;
the time step size $\Delta t$, the value of strides $k$, the range of initial conditions (IC), and the numerical integration scheme used for data generation are summarized in Table 7.

Implementation details.
As a decoder-only transformer model, FMint is configured with approximately 15.8M parameters.
The model features six heads for multi-head attention, with an input/output dimension of 256 for each layer.
The embedding vector dimensions for the coarse solution, error correction, and query are set to 256, while the hidden dimension of the feed-forward networks is set to 1024.
Pretraining is conducted on a NVIDIA A100 GPU with 80 GB of memory and finetuning is conducted on a NVIDIA A6000 GPU with 48 GB of memory.
We use AdamW optimizer with a warmup-cosine-decay schedule, with peak learning rate 1e-4 and 60 training epochs.
The Adam $\beta_{1}$ and Adam $\beta_{2}$ are 0.9 and 0.999, respectively, and the weight decay is set to 1e-4.
Demo number for training and testing is five.

Baselines and tasks. We compare FMint with three baseline models: Neural ODE [], NeurVec [], and In-Context Operator Networks (ICON-LM) []. Neural ODEs model continuous-time dynamics by parameterizing the derivative of the hidden state with a neural network. This approach turns the forward pass into solving an initial value problem, offering a memory-efficient way to capture temporal patterns in data.
NeurVec is a deep learning-based corrector aimed to compensate for integration errors and enable larger time step sizes in simulations.
ICON-LM is a foundation model for operator learning, achieving better accuracy and data efficiency than Fourier Neural Operator [] and DeepONet [].

Among them, ICON-LM is a multi-task model trained on a large collection of examples while both Neural ODE and NeurVec fit one neural network per example. The configuration and training details of Neural ODE and NeurVec are provided in Appendix A.3 and A.4, while the settings for ICON-LM follow those outlined in the original work [].

Evaluation metrics. We use the mean absolute errors (MAE) and root mean square errors (RMSE) compared to fine-grained ODE solutions as the evaluation metric:

$$ $\text{MAE}=\frac{1}{N}\sum_{i=1}^{N}\|\tilde{{\bm{u}}}^{i}-{\bm{u}}^{i}\|_{2}, \quad\text{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\|\tilde{{\bm{u}}}^{i}-{\bm{u} }^{i}\|_{2}},$ (9) $$

where $\tilde{{\bm{u}}}^{i}$ is the predicted ODE solution for the $i$-th equation. For FMint, it can be computed via error correction $\tilde{{\bm{u}}}^{i}_{k}=\hat{{\bm{u}}}^{i}_{k}+\hat{\text{err}}^{i}$ such that $\hat{{\bm{u}}}^{i}_{k}$ is the coarse solution of the $i$-th equation, and $\hat{\text{err}}^{i}$ is the model output by FMint.

### 3.2 RQ1. FMint v.s. baselins models

Here we show that FMint outperforms baseline methods on both in-distribution and out-of-distribution data in terms of accuracy and efficiency.

In-distribution ODEs

Accuracy. We evaluate FMint on the test split of the pretraining dataset. This contains ODEs from the same ODE families with the same parameter range, but with different random parameters within the range and different initial conditions. For each ODE system, we finetuned our model for 1000 or 2000 epochs with trajectories of size 100. Details on the testing coefficient range and finetuning epochs are included in Table 8.

**Table 2: Performance for in-distribution ODEs via MAE and RMSE (lower is better)**
|  | MAE | RMSE |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Methods | Lorenz | Damped Osci | Van der Pol | LV | Lorenz | Damped Osci | Van der Pol | LV |
| ICON-LM | 2.87 | 0.095 | 1.73 | 5.74 | 4.98 | 0.14 | 2.03 | 9.07 |
| Neural ODE | 17.9017 | 0.0865 | 1.8476 | 6.3497 | 22.6992 | 0.3907 | 2.2656 | 9.2982 |
| NeurVec | 6.6465 | 0.0919 | 1.7583 | 4.9032 | 10.5565 | 0.4499 | 2.0390 | 8.5570 |
| FMint | 0.159 | 0.0044 | 0.00878 | 0.0304 | 0.33 | 0.00695 | 0.0119 | 0.0433 |

MAE and RMSE results are shown in Table 2 for all in-distribution ODE families with the best result in bold and the second best with underline.
Both metrics are averaged over 25 ODEs with initial conditions from the same family and the number of demos is five during the inference stage.
Example visualization of FMint’s performance on Lorenz and Van der Pol is shown in Figure (2(a)) (2(b)). The fine-grained solution $u_{j}$ is labeled as Fine ode, the coarse solution $\hat{u}_{kj}$ is labeled as coarse ode and FMint result is labeled as FMint ode.
For visualization of FMint on all tested ODEs, see Appendix A.7.

Figure: (a) Lorenz
Refer to caption: x2.png

We observe that FMint achieves an accuracy that is at least an order of magnitude higher than all baseline models. For the relatively simple model of damped oscillators, we observe that all models exhibit relatively small errors, with FMint achieving the highest accuracy. However, for the Lorenz model, where the dynamics are highly chaotic, Neural ODE completely fails to capture the dynamics at large time steps. This underscores the importance of human-designed algorithms in stabilizing the system, especially when the initial conditions vary. NeurVec performs better but remains unsatisfactory, while ICON-LM achieves the best results among the baselines. These observations suggest that pretraining on a diverse dataset is crucial; otherwise, variations in the initial conditions will cause the model to fail, which explains NeurVec’s weaker performance in this case. FMint achieves the best performance due to both its pretraining and the combination of data-driven and human-designed algorithms.

Efficiency.
We further examined the runtime of FMint in comparison with fine solution generation using RK4. The test is conducted on Lotka-Volterra system with 500 equations.
To display the runtime better, we use the runtime for obtaining coarse solutions using RK4 as one unit, and we report the result in Figure 3.
FMint is able to attain results with comparable accuracy to the fine solutions (RK-FINE) using less than 20% of its time.

Out-of-distribution ODEs

Figure: Figure 3: Normalized runtime for Lotka-Volterra of reference solution (RK4-FINE) coarse solution (RK4-COARSE) and FMint.
Refer to caption: extracted/5891009/FMint_exp_fig/FMInt_results_ode/runtime.jpg

To further demonstrate the performance of FMint on unseen ODEs from different ODE families, we use data simulated from the following dynamical systems: driven-damped pendulum (2D), falling object with air resistence (2D), FitzHugh-Nagumo systems (2D), Pendulum under gravity (2D), and the Rössler dynamics (3D). These test ODEs are qualitatively different from the training data, thus convincing to validate FMint’s capability as a foundation model for dynamical systems.
For more details about the OOD ODEs, see Appendix A.2.
We evaluate the transfer performance of FMint using trajectories of size $N\in\{25,50,200,1000\}$.
We fine-tune the model for 1000 to 2000 iterations on the new data of size $N$ and report the RMSE and MAE on the test data. Details on the coefficient range and finetuning epochs are included in Table 8.
As a comparison, RMSE and MAE are computed for NeurVec and Neural ODE using training set of size 50K. The results are shown in Table 3.2

We observe that with only 25 training trajectories, FMint is able to outperform NeurVec and Neural ODE trained on data with sample size 50000.
Although increasing the training sample size may improve the accuracy of FMint, finetuning with small sample size still gives us robust results.
This shows that FMint has great potential for large-scale simulations in many real-world scenarios where data collection is expensive.

Figure: Table 3: Performance of FMint on unseen ODEs in MAE.
Refer to caption: x4.png

MAE v.s. training epochs.
To better examine what factors may affect FMint’s performance for unseen ODEs other than training sample size, we fine tuned our model on Lorenz equation with unseen coefficient range $\rho\in(100,150)$. The system exhibits more complex chaotic dynamics with larger attractor regions and more intricate patterns under this range. We plotted the MAE on the testing split of trajectories of size 100 with respect to finetuning epochs in Figure 4 with one standard deviation. As expected, we see that MAE decreases as training epochs increases but starts to saturate after a certain point.

### 3.3 RQ2. Uniformity of FMint’s performance

On dynamical systems with qualitatively different behaviors.

For the same dynamical systems, different parameter ranges can lead to vastly different behaviors. For example, the Lorenz system is highly sensitive to initial conditions and parameter values, leading to chaotic behavior under certain conditions.
This is challenging for traditional numerical solvers to simulate with large time steps, since even small errors introduced by large time steps can result in significant deviations over time due to the sensitivity to initial conditions.
In addition, chaotic systems exhibit fine-scale behaviors, such as oscillations, sharp changes, or bifurcations.
Large time steps may skip over these critical behaviors, missing important features like turning points, sharp gradients, or periodicity.

Figure: Table 4: Performance of FMint on systems exhibiting qualitatively different behaviors, measured in MAE and RMSE. Driven damped pendulum and FitzHugh-Nagumo are shorten to DDP and FHN.
Refer to caption: x5.png

For example, in Lorenz system, when $\rho\in[13,24.74)$, the system exhibits converging oscillatory behavior toward either of the two fixed points. But when $\rho\in[24.74,100)$, the system exhibits chaotic behavior and sensitive to initial conditions. As $\rho$ increases further e.g., $\rho>100$, the chaotic behavior becomes more complex.
Figure 5 shows trajectories generated with $\sigma=12.69,\beta=2.59,\rho=24.33$ under two initial conditions $[1,1,1]$ and $[20,20,20]$ (top) and trajectories generated with $\sigma=12.69,\beta=2.59,\rho=78.06$ under two initial conditions $[1,1,1]$ and $[20,20,20]$ (bottom). Both trajectories are simulated for 10000 steps with $\Delta t=0.004$.

Here we conduct extensive experiments of FMint in such cases to test its generalization ability. Conretely, we test on three of Lorenz, three of Driven damped pendulum (DDP), four of FitzHugh-Nagumo (FHN) and three of Rössler that exhibits various behaviors under different range of coefficients.
Results are reported in Table 3.3.
For more details on the coefficient range for each behavior category, see Appendix A.8.

On pretraining ODEs with different strides.

In addition, we test our pretrained model on test data simulated using different strides $k$ than the pretraining data. This examines how adaptable FMint is to handle realistic circumstances in which the coarse solution simulation varies during the inference stage.
For consistency over various families, we generate test examples with new stride values proportional to the training strides: $\alpha k,\alpha=\{0.5,2.0\}$.

**Table 5: MAE and RMSE of FMint under strides differs from pretraining data.**
| Name | $k$ | MAE | RMSE | $k$ | MAE | RMSE |
| --- | --- | --- | --- | --- | --- | --- |
| Lorenz model | 50 | 4.5e-2 | 9.4e-2 | 200 | 3.15e-2 | 4.65e-2 |
| Damped oscillator | 50 | 3.58e-3 | 5.51e-3 | 200 | 4.03e-3 | 6.54e-3 |
| Van der Pol | 5 | 1.6e-3 | 2.11e-3 | 20 | 2.76e-3 | 3.60e-3 |
| Lotka-Volterra | 50 | 6.51e-3 | 1.01e-2 | 200 | 1.86e-2 | 2.64e-2 |

For each ODE system, we finetuned our model for 1000 or 2000 epochs with trajectories of size 100.
Table 5 reports MAE and RMSE of FMint for smaller or larger strides $k$. We observe that the accuracy of all four ODE families is consistent with the accuracy of the training strides.

### 3.4 RQ3. Multi-modal FMint learning

Performance of FMint with the supplemental textual data.

We further investigated the performance of FMint with the supplemental textual data for both pretraining and out-of-distribution ODEs.
We first pretrained the model on the pretraining data but with prepared textual data for each ODEs. Then, the model is further finetuned on each ODEs for 1000 or 2000 epochs on trajectories of size 100. For consistency, we used the exact same numerical dataset for both pretraining FMint without prompts and finetuning as in Section 3.2. For details on the generation of the textual data, see Appendix A.5.

Table 6 shows the MAE and RMSE of FMint with and without the prompts. We observe that FMint, combined with prompts, enhances performance in challenging dynamical systems such as the Lorenz, Van der Pol, Lotka-Volterra, and Rössler systems. However, in relatively simpler cases like the Damped Oscillator and Pendulum under gravity, the inclusion of prompts does not improve performance, remaining comparable to baseline FMint results. This suggests that the multi-modal foundation model is particularly advantageous when simulating complex and challenging systems.

**Table 6: Comparison of FMint and FMint with textual data in MAE and RMSE**
|  | MAE |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Methods | Lorenz | Damped Osci | Van der Pol | LV | Driven damped | Falling | Fitzhugh Nagumo | Pendulum | Rössler |
| FMint | 0.159 | 4.4e-3 | 8.78e-3 | 3.04e-2 | 6.77e-4 | 2.99e-3 | 9.79e-3 | 1.26e-3 | 1.61e-3 |
| FMint + prompt | 0.109 | 4.52e-3 | 4.29e-3 | 1.76e-2 | 8.12e-4 | 2.66e-3 | 9.90e-3 | 2.84e-3 | 1.40e-3 |
|  | RMSE |  |  |  |  |  |  |  |  |
| Methods | Lorenz | Damped Osci | Van der Pol | LV | Driven damped | Falling | Fitzhugh Nagumo | Pendulum | Rössler |
| FMint | 0.33 | 6.95e-3 | 1.19e-2 | 4.33e-2 | 9.68e-4 | 5.26e-3 | 5.8e-2 | 1.7e-3 | 2.56e-3 |
| FMint + prompt | 0.224 | 6.97e-3 | 5.87e-3 | 2.54e-2 | 1.17e-3 | 4.35e-3 | 9.90e-3 | 2.78e-3 | 2.19e-3 |

## 4 Related Work

Neural network for dynamical systems.
In recent years, neural network-based solvers have been widely applied to scientific problems such as solving ODEs, PDEs, operator learning, and inverse problems. One common approach parameterizes PDE solutions using feed-forward neural networks [, , , , , , , , ], where physical laws are enforced via hard or soft constraints in the loss function. The Finite Expression Method (FEX) [, , ] offers an interpretable alternative by representing PDE solutions in computer algebra, capturing solution structure with high accuracy. Neural operator learning [, , , , , ] maps varying parameters or initial conditions to solutions, achieving discretization invariance but requiring large amounts of high-quality data and lacking generalization to unseen distributions.

Recent research has explored integrating traditional numerical algorithms with deep learning to improve the accuracy of dynamical system simulations [, ]. For instance, NeurVec [] enables rapid ODE simulations with large time steps, achieving decent accuracy on several classic systems. However, its limited generalization to out-of-distribution systems restricts its practicality for large-scale real-world simulations.

Foundation model in scientific machine learning.

Large language models like GPT-4 [], DALL-E [], and Llama [] have achieved remarkable success across various domains [, , , , , , , ], including text-to-visual generation [, ] and information retrieval []. These models leverage extensive pre-training and adapt to downstream tasks via zero-shot or few-shot learning [, ], demonstrating impressive transfer learning capabilities. Inspired by these advances, foundation models have gained traction in scientific machine learning. For example, Subramanian
et al. [] explored the Fourier Neural Operator (FNO) [] for solving classical PDEs, while the Unified PDE Solver (UPS) [] extended this approach to 1D and 2D PDEs using pre-trained models. Additionally, McCabe et al. [] embedded PDEs with varying properties into a shared space, and Rahman et al. [] enhanced attention mechanisms for handling PDEs with diverse dimensions.

A growing area in scientific machine learning is in-context learning [, , , , ], where models are prompted with example pairs and trained to predict new queries based on observed patterns. The In-context Operator Network (ICON) [, , ] applies this approach to operator learning by leveraging example pairs with varying PDE parameters and solutions to predict solutions for new query data.

## 5 Conclusion

In this paper, we presented FMint, a novel multi-modal foundation model that speeds up large-scale simulations of dynamical systems.
Based on the architecture of decoder-only transformer, FMint takes textual and numerical data to deliver high-accuracy simulation based on initial solution from traditional numerical solvers. FMint incorporates the in-context learning for a universal error corrector for ODEs from given prompted sequences of coarse initialized solutions.
It is pre-trained using a diverse set of ODE families with qualitatively different behaviors.

We show that FMint achieves a significant improvement in accuracy over state-of-the-art dynamical system simulators and accelerates traditional integration schemes.
In comparison to direct ODE solvers,
we recognize the importance of integrating the strengths of human-designed algorithms and data-driven methods for the simulation of dynamical systems. The in-context learning scheme enables it to effectively interpolate to arbitrary time points, enhancing its versatility in handling temporal dynamics. Furthermore, we propose a multi-modal foundation model perspective to address scientific computing problems, challenging the mainstream numerical-data-only approach in the field. Given FMint’s performance, it shows promise for scaling up to simulate even more complex dynamics in real-world applications.

## Acknowledgments

H. Y. and Z. S. were partially supported by the US National Science Foundation under awards DMS-2244988, DMS-2206333, the Office of Naval Research Award N00014-23-1-2007, and the DARPA D24AP00325-00. J. Y. was partially supported by AFOSR
MURI grant FA9550-20-1-0397.

## Appendix A Appendix

### A.1 Pretraining ODE details

- •
The Lorenz Model is described by
$\left\{\begin{aligned} \frac{dx}{dt}&=\sigma(y-x),\\
\frac{dy}{dt}&=x(\rho-z)-y,\\
\frac{dz}{dt}&=xy-\beta z,\end{aligned}\right.$
where $x$, $y$, and $z$ represent the convection rate, horizontal temperature variation, and vertical temperature variation, respectively. The Prandtl number, $\sigma$, is a dimensionless quantity indicating the ratio of momentum diffusivity to thermal diffusivity, while the Rayleigh number, $\rho$, quantifies the temperature gradient driving convection. The parameter $\beta$ represents the geometric aspect ratio of the convective cells in the model.
- •
Damped oscillator equation is given by:
$\frac{d^{2}x}{dt^{2}}+2\zeta\omega\frac{dx}{dt}+\omega^{2}x=0,$
where $\zeta$ is the damping ratio, and $\omega$ is the natural frequency.
- •
The Van der Pol Oscillator is given by the following equation:
$\frac{d^{2}x}{dt^{2}}-\mu(1-x^{2})\frac{dx}{dt}+x=0,$
where $\mu$ controls the degree of nonlinearity and damping.
- •
The Lotka-Volterra system is given by the following equations:
$\left\{\begin{aligned} \frac{dx}{dt}&=\alpha x-\beta xy,\\
\frac{dy}{dt}&=\delta xy-\gamma y,\end{aligned}\right.$
where $x$ is the number of prey, $y$ is the number of predator, $\alpha$ is the natural growing rate of prey in the absense of predators, $\beta$ is the natural dying rate of prey due to predation, $\gamma$ is the natural dying rate of predators in the absence of prey, and $\delta$ is the rate at which predators increase by consuming prey.

### A.2 Testing ODE details

- •
Rössler attractor is a system of three nonlinear ODEs, which exhibits chaotic dynamics. The equation is given by
$\left\{\begin{aligned} \frac{dx}{dt}&=-y-z,\\
\frac{dy}{dt}&=x+ay,\\
\frac{dz}{dt}&=b+z(x-c),\end{aligned}\right.$
where $x,y,z$ are state variables, $a,b,c$ are parameters that control the behavior of the system.
- •
FitzHugh-Nagumo model is used for neuron activity dynamics. It captures the essential features of neuronal excitability, including the generation and propagation of action potentials. The equation is given by
$\left\{\begin{aligned} \frac{dv}{dt}&=v-\frac{v^{3}}{3}-w+I,\\
\frac{dw}{dt}&=\epsilon(v+a-bw),\end{aligned}\right.$
where $v$ is membrane potential, $w$ is a recovery variable, $I$ is an external stimulus current. Parameter ranges are listed as follows: $\epsilon$, $a$, $b$ and $I$.
- •
This equation describes the dynamics of the falling object with air resistance:
$\frac{d^{2}x}{dt^{2}}=g-c\frac{dx}{dt},$
where $g$ represents the gravitational constant, and $c$ is the coefficient of drag.
- •
This equation describes how the pendulum swings under the influence of gravity and damping:
$\frac{d^{2}x}{dt^{2}}=-\frac{g}{l}\sin(x)-b\frac{dx}{dt},$
where $g$ represents the gravitational constant, and $l$ is the length of pendulum, and $b$ is the damping coefficient.
- •
This is a classic problem in the study of dynamical systems and chaotic behavior, particularly under the influence of non-linear restoring forces and external driving forces. It is written as:
$\frac{d^{2}\theta}{dt^{2}}+b\frac{d\theta}{dt}+c\sin(\theta)=A\cos(\omega t),$
where $\theta$ represents the angular displacement from the vertical, $b$ is the damping coefficient, $c$ is the gravitational constant times the length of the pendulum, $A$ is the amplitude, and $\omega$ is the angular frequency of the drive.

**Table 7: Data generation set-up**
| Name | $k$ | $\Delta t$ | IC (1st dim) | IC (2nd dim) | IC (3rd dim) | Integration |
| --- | --- | --- | --- | --- | --- | --- |
| Lotka-Volterra | 100 | 0.001 | $(10,100)$ | $(5,10)$ | NA | RK4 |
| Van der Pol | 10 | 0.01 | $(-1.0,1.0)$ | $(-0.5,0.5)$ | NA | RK4 |
| Damped Osci. | 100 | 0.0001 | $(-2.0,2.0)$ | $(-0.1,0.1)$ | NA | RK4 |
| Lorenz | 100 | 0.0001 | $(-5,5)$ | $(-5,5)$ | $(0,25)$ | RK4 |
| Fitzhugh Nagumo | 100 | 0.005 | $(-1.0,1.0)$ | $(-0.5,0.5)$ | NA | RK4 |
| Falling object | 20 | 0.01 | $(0,100)$ | $(0,2)$ | NA | RK4 |
| Pendulum gravity | 20 | 0.01 | $(0,\frac{\pi}{4})$ | $(-\frac{\pi}{4},\frac{\pi}{4})$ | NA | RK4 |
| Driven damped pendulum | 20 | 0.01 | $(-\frac{\pi}{4},\frac{\pi}{4})$ | $(-0.5,0.5)$ | NA | RK4 |
| Rössler | 100 | 0.001 | $(-1,1)$ | $(-1,1)$ | $(-1,1)$ | RK4 |

**Table 8: Coefficient range for pre-training and testing and epochs for fine-tuning**
| Name | pretraining | testing | finetuning epochs |
| --- | --- | --- | --- |
| Lotka-Volterra | $\alpha\in[0.1,1.0]$ | $\alpha\in[0.05,1.2]$ | 2000 |
|  | $\beta\in[0.01,0.1]$ | $\beta\in[0.005,0.15]$ |  |
|  | $\gamma\in[0.1,1.0]$ | $\gamma\in[0.05,1.2]$ |  |
|  | $\delta\in[0.01,0.1]$ | $\delta\in[0.005,0.15]$ |  |
| Van der Pol | $\mu\in[0.1,5]$ | $\mu\in[0.05,10]$ | 1000 |
| Damped Osci. | $\zeta\in[0.1,2.0]$ | $\zeta\in[0.1,7.0]$ | 1000 |
|  | $\omega\in[0.5,5.0]$ | $\omega\in[0.5,8.0]$ |  |
| Lorenz | $\sigma\in[8,12]$ | $\sigma\in[5,15]$ | 2000 |
|  | $\rho\in[20,30]$ | $\rho\in[15,35]$ |  |
|  | $\beta\in[2,3]$ | $\beta\in[1.5,3.5]$ |  |
| Falling object | NA | $c\in[0.01,2.0]$ | 1000 |
| Fitzhugh Nagumo | NA | $\epsilon\in[0.005,0.15]$ | 2000 |
|  |  | $a\in[0.4,1.2]$ |  |
|  |  | $b\in[0.4,1.2]$ |  |
|  |  | $I\in[0.3,2.0]$ |  |
| Pendulum gravity | NA | $l\in[0.5,2]$ | 1000 |
|  |  | $b\in[0.05,1]$ |  |
| Driven damped pendulum | NA | $A\in[0.1,2]$ | 1000 |
|  |  | $b\in[0.1,2]$ |  |
|  |  | $c\in[1,5]$ |  |
|  |  | $\omega\in[0.5,3]$ |  |
| Rössler | NA | $a\in[0.05,0.35]$ | 2000 |
|  |  | $b\in[0.05,0.35]$ |  |
|  |  | $c\in[3.0,10.0]$ |  |

### A.3 Neural ODE Implementation details

Neural ODE Architecture.
The neural ODE model employed in our experiments consists of a fully connected feedforward neural network with the following architecture:

- •
Input layer: The state dimension of the system (varies per system).
- •
Hidden layers: Two hidden layers, each with 1024 neurons and Tanh activation functions.
- •
Output layer: The derivative of the state with respect to time (same dimension as the input state).

The model is formally described as:

$$ $f_{\theta}(t,{\bm{x}})=\text{Linear}(\text{Tanh}(\text{Linear}({\bm{x}}))).$ $$

Neural ODE Training Details.
We trained the neural ODE models on various dynamical systems using fine time step data and evaluated them on coarse time step data. The specific training parameters were as follows:

- •
Learning Rate: 0.001
- •
Optimizer: Adam
- •
Learning Rate Decay: StepLR with a decay factor of 0.5 every 20 epochs
- •
Batch Size: 500
- •
Number of Epochs: 100
- •
Loss Function: Mean Squared Error (MSE)

### A.4 NeurVec Implementation details

NeurVec Model Architecture.
The NeurVec model incorporates a multi-layer perceptron (MLP) with a custom activation function, defined as follows:

- •
Input layer: The state dimension of the system .
- •
Hidden layer: 1024 neurons with a custom rational activation function.
- •
Output layer: The error correction term for the state update.

The rational activation function is defined by:

$$ $f(x)=\frac{a_{3}x^{3}+a_{2}x^{2}+a_{1}x+a_{0}}{b_{2}x^{2}+b_{1}x+b_{0}},$ $$

where the parameters are: $a_{0}=0.0218,a_{1}=0.5000,a_{2}=1.5957,a_{3}=1.1915,b_{0}=1.0000,b_{1}=0.0000,$ and $b_{2}=2.3830$.

NeurVec Training Details.
We trained the NeurVec model on coarse time step data derived from various dynamical systems. The training parameters were as follows:

- •
Learning Rate: 0.001
- •
Optimizer: Adam
- •
Learning Rate Decay: cosine annealing schedule
- •
Batch Size: 500
- •
Number of Epochs: 100
- •
Loss Function: Mean Squared Error (MSE)

**Table 9: Performance of FMint on unseen ODEs in RMSE**
|  |  | Unseen ODE |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Method | # Samples | Driven damped | Falling | Fitzhugh Nagumo | Pendulum | Rossler |
| FMint | 25 | 1.64e-3 | 4.18e-3 | 5.59e-2 | 1.55e-3 | 2.57e-3 |
|  | 50 | 9.68e-4 | 5.26e-3 | 5.8e-2 | 1.7e-3 | 2.56e-3 |
|  | 200 | 1.34e-3 | 4.74e-3 | 5.3e-2 | 1.6e-3 | 2.5e-3 |
|  | 1000 | 9.02e-4 | 3.7e-3 | 4.7e-2 | 1.69e-3 | 2.4e-3 |
| NeurVec | 50000 | 0.0741 | 0.8351 | 0.3859 | 0.0974 | 0.0265 |
| Neural ODE | 50000 | 0.6667 | 1.3754 | 1.2139 | 1.0763 | 0.1571 |

### A.5 Textual data generation

We generated the supplemental textual data for each ODE systems with the assistance of GPT-4. For any ODE systems that we are interested in, we use the following prompt for generation.

Please use [numbers required] different ways to explain what [ODE system], [mathematical formula] is, where it can be applied and what the meaning of [parameters]. Also include what types of behavior will there be for different ranges of its parameters.

Take Lorenz attractor as an example, we use the following prompt:

Please use 30 different ways to explain what Lorenz attractor, $$\backslash$frac{dx}{dt} = $\backslash$sigma (y - x), $\backslash$frac{dy}{dt} = x ($\backslash$rho - z) - y, $\backslash$frac{dz}{dt} = xy - $\backslash$beta$ is, where it can be applied and what the meaning of $$\backslash$sigma$, $$\backslash$rho$ and $$\backslash$beta$. Include what types of behavior will there be for different ranges of its parameters.

We have generated 30 textual data for each ODE systems, with slight manual adjustments applied after GPT-4 generation. For pretrained ODEs, we reserved placeholders for the actual parameters used for training and testing. For out-of-distribution ODEs, we did not specify the exact parameters used in our textual data.

During pretraining stage, the textual data is randomly selected from all the provided prompts that correspond to the ODE system. The placeholders for parameters are replaced with the actual parameter values associated with the numerical data. During inference stage, textual data is also randomly selected but only the pretraining ODE systems are filled with exact parameters.

### A.6 Textual data example

Here we show several examples of the textual data for pretraining ODEs and out-of-distribution ODEs: pretraining ODE Lorenz attractor and out-of-distribution ODE Rössler attractor:

- •
Prompts for Lorenz attractor:
It models the evolution of three variables—$x,y,z$ over time, governed by three parameters: $$\backslash$sigma = 12.69, $\backslash$rho = 34.297, $\backslash$beta = 2.592$.
The Lorenz system of equations $$\backslash$frac{dx}{dt} = $\backslash$sigma (y - x), $\backslash$frac{dy}{dt} = x ($\backslash$rho - z) - y, $\backslash$frac{dz}{dt} = xy - $\backslash$beta$ serves as a classic example in chaos theory, showing how a deterministic system can exhibit aperiodic, non-repeating, and sensitive trajectories despite being governed by simple rules.
- •
Prompts for Rössler attractor:
The Rössler system $$\backslash$dot{x} = -y - z, $\backslash$dot{y} = x + a y, $\backslash$dot{z} = b + z (x - c)$ is a set of three differential equations used to model chaotic behavior. It has parameters $a$, $b$ and $c$ which control the system’s dynamics.
When $c$ is low, oscillations are periodic; as $c$ increases, the system becomes chaotic, which can be used to design chaotic communication systems.

### A.7 Visualization of FMint on all ODEs

Figure: (a) Lotka Volterra
Refer to caption: x7.png

Figure: (a) Rössler attractor
Refer to caption: x9.png

Figure: (a) Driven damped pendulum
Refer to caption: x11.png

### A.8 Details on coefficient range for each behavior category

Lorenz oscillatory. $\sigma\in[5,15],\rho\in[13,24.74],\beta\in[1.5,3.5]$.The system exhibits converging oscillatory behavior toward either of the two fixed points.

Lorenz chaotic. $\sigma\in[5,15],\rho\in[24.74,100],\beta\in[1.5,3.5]$. This is the famous region where the Lorenz attractor forms, with characteristic of never repeating chaotic motion and sensitive dependence on initial conditions.

Lorenz chaotic (complex). $\sigma\in[5,15],\rho\in[24.74,150],\beta\in[1.5,3.5]$. It contains more complicated chaotic dynamics with larger attractor regions and more intricate patterns.

Driven damped pendulum underdamped. $b\in[0.01,0.1],c\in[1.0,2.0],A\in[0.1,0.5],\omega\in[0.5,3.0]$. In this parameter range, the underdamped pendulum may achieve steady oscillations.

Driven damped pendulum overdamped. $b\in[0.5,2.0],c\in[1.0,2.0],A\in[0.1,0.5],\omega\in[0.5,3.0]$. It occurs when the damping is very strong. The pendulum returns to its equilibrium position without oscillating.

Driven damped pendulum chaotic. $b\in[0.1,0.2],c\in[1.0,2.0],A\in[1,3.0],\omega\in[0.5,3.0]$. The chaotic behavior emerges and the system exhibits nonlinear and non-periodic motion. It is highly sensitive to initial conditions.

FitzHugh-Nagumo resting. $I\in[0.0,0.1],\epsilon\in[0.08,0.1],a\in[0.5,0.7],b\in[0.1,0.2]$. In this parameter range, the membrane potential stays close to its equilibrium value.

FitzHugh-Nagumo spikes. $I\in[0.1,0.5],\epsilon\in[0.01,0.5],a\in[0.7,1.0],b\in[0.2,0.25].$ The system can exhibit a single action potential or spike and followed by recovery.

FitzHugh-Nagumo bursting. $I\in[0.5,1.5],\epsilon\in[0.01,0.02],a\in[0.9,1.1],b\in[0.15,0.2].$ In this state, the system may have transition between excitable and oscillatory behavior, with possible periodic spiking.

FitzHugh-Nagumo oscillatory. $I\in[0,2],\epsilon\in[0.01,0.1],a\in[0.5,1.2],b\in[0.1,0.3]$. The system can exhibit bistability, or mixed mode oscillations.

Rössler periodic. $a\in[0.1,0.2],b\in[0.1,0.2],c\in[4,5].$ The trajectory in phase space forms a closed loop, and we may see that the system returns to the same state after a fixed period.

Rössler chaotic. $a\in[0.2,0.3],b\in[0.2,0.25],c\in[5,9]$. In this state, the system exhibits unpredictable behavior that is sensitive to the initial conditions.

Rössler hyperchaos. $a\in[0.25,0.4],b\in[0.25,0.3],c\in[9,13].$ The system is chaotic under this coefficient range in more than one direction, and may have even more complex and unpredictable behavior than standard chaos.

Figure: Figure 9: Lorenz attractors - oscillatory. Trajectories generated with $\sigma=12.69,\beta=2.59,\rho=24.33$ under two initial conditions $[1,1,1]$ and $[20,20,20]$ for 10000 steps with $\Delta t=0.004$.
Refer to caption: x13.png

Figure: Figure 10: Lorenz attractors - chaotic. Trajectories generated with $\sigma=12.69,\beta=2.59,\rho=78.06$ under two initial conditions $[1,1,1]$ and $[20,20,20]$ for 10000 steps with $\Delta t=0.004$.
Refer to caption: x14.png

Figure: Figure 11: Lorenz attractors - chaotic (complex). Trajectories generated with $\sigma=12.69,\beta=2.59,\rho=78.06$ under two initial conditions $[1,1,1]$ and $[20,20,20]$ for 10000 steps with $\Delta t=0.004$.
Refer to caption: x15.png

Figure: Figure 12: Driven damped pendulum - underdamped. Trajectories generated with $b=0.005,c=1.0,A=0.25,\omega=2.0$ under two initial conditions $[0.1,0]$ and $[1.0,0]$ for 10000 steps with $\Delta t=0.01$.
Refer to caption: x16.png

Figure: Figure 13: Driven damped pendulum - overdamped. Trajectories generated with $b=1.0,c=1.0,A=0.25,\omega=2.0$ under two initial conditions $[0.1,0]$ and $[1.0,0]$ for 10000 steps with $\Delta t=0.01$.
Refer to caption: x17.png

Figure: Figure 14: Driven damped pendulum - chaotic. Trajectories generated with $b=0.15,c=1.0,A=0.25,\omega=2.0$ under two initial conditions $[0.1,0]$ and $[1.0,0]$ for 10000 steps with $\Delta t=0.01$.
Refer to caption: x18.png

Figure: Figure 15: Rössler - periodic. Trajectories generated with $a=0.179,b=0.159,c=4.81$ under two initial conditions $[1,1,1]$ and $[5,5,5]$ for 10000 steps with $\Delta t=0.01$.
Refer to caption: x19.png

Figure: Figure 16: Rössler - chaotic. Trajectories generated with $a=0.368,b=0.279,c=12.24$ under two initial conditions $[1,1,1]$ and $[5,5,5]$ for 10000 steps with $\Delta t=0.01$.
Refer to caption: x20.png

## Peer Reviews

confidence: 2
contribution: 2
initial_rating: 3
presentation: 3
questions: 1.The details of the natural language prompt are lacking. What kind of embedding was used?

2.A plausible hypothesis is that the language prompt mainly helps the model differentiate between system types. Without it, systems like LV and the Van der Pol oscillator, which have similar shapes, might be confused with each other if only numeric prompts are used. A useful experiment might be to use simple category prompts (e.g., one for each system type) as a control, to demonstrate that complex natural language prompts indeed add value and that the model extracts meaningful information from them.

3.Can this model truly leverage information as specific as the symbolic form of the system’s control equations? Are there further experiments to support this?

4.I suggest including the MAE for a coarse prediction as a comparison.

5.Transformers are powerful but computationally intensive. Has there been any consideration of alternative architectures that could achieve similar results with lower computational overhead?

6.FMint’s performance is primarily compared to Neural ODE, NeurVec, and ICON-LM. Could additional benchmarks, like recent neural operators, provide further insights into its relative strengths and limitations?
soundness: 3
strengths: 1.The FMint model uniquely combines numerical and textual data, enhancing its robustness in diverse ODE environments

2.The use of in-context learning for error correction in dynamical systems showcases a novel integration of human-designed and data-driven algorithms, enhancing the model’s ability to handle data-scarce environments
summary: The paper presents FMint, a multi-modal foundation model for fast simulation of dynamical systems. It combines human-designed and data-driven approaches using a decoder-only transformer architecture with in-context learning, leveraging both numerical and textual data for error correction. Pre-trained on 40,000 ordinary differential equations (ODEs), FMint demonstrates significant improvements in accuracy and efficiency, achieving 1 to 2 orders of magnitude better accuracy and a 5X speedup over traditional numerical solvers, making it a promising general-purpose solver for dynamical systems.
weaknesses: 1.Some details in the paper are insufficiently described. For example, the length of the data trajectories and the number of pre-training contexts are not specified. These details are crucial to the experimental design as they relate to the transformer’s input length and memory usage.

2.While the model benefits from multi-modal input, its performance may decrease in the absence of textual data, especially for complex systems. How does the model perform on high-dimensional ODEs with limited numerical data and no textual data?

3.Since this is in-context learning, it would be beneficial to conduct few-shot experiments without fine-tuning and to examine how context length impacts results.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: - Table 2. Most of the baselines scores are very bad (>1). Is it because such models should not be used for such tasks? Or is it because the training didn't go well in your experiments? 
- Table 2. I am having a hard time convincing myself that the comparison with the baselines is fair. Does the baseline take as input the coarse-grained trajectory as well? It seems to me that this is not the case (take for example neuralODE). Then, I am not surprised that a network which takes as an input the coarse-grained trajectory is performing better.
soundness: 2
strengths: The idea of integrating textual information with numerical coarse-grained trajectories to simulate dynamical systems is definitely interesting.
summary: For the simulation of dynamical systems, the paper proposes a multi-modal model that processes both the coarse-grained trajectories and textual data describing the system generated by ChatGPT to obtain fine-grained trajectories. The model is tested on various dynamical systems including Lorenz attractors, driven damped pendulum, FitzHugh-Nagumo and Rössler systems.
weaknesses: - the authors make strong claims that are hardly substantiated in the paper, in particular the comparison in terms of accuracy is hardly a fair comparison 
	- l103: "the first multi-modal foundation model that synthesizes human-designed algorithms and deep learning framework" is too generic and vague, since "human-designed" covers both the standard numerical solvers and machine-learning models 
	- l106: the better accuracy, is not correctly substantiated by Table 2. as explained in the following section "Questions"
- the paper contains numerous overstatements and confusing sentences, making it difficult to understand the core arguments. The paper contains too many generic statements in the flavor of arguing that the method is in general "better and faster", see the minor points below 

Minor comments 
- l025: "corpus" is used for textual data, but I guess you also refer to numerical data in that sentence
- l027: "that exhibit ... high dimensionality" could be better written
- l047: what is the meaning of "format" in that context?
- l053: "human-designed algorithm" any model you mention in the paper is "human-designed", this expression does not add anything 
- l054-055: there is a repetition "real-world scenarios" "various applications"
- l068-069: the question is not well formulated. For many equations, the classical solvers appear to be actually both more accurate and more efficient than their neural network counterpart. However, neural network used to solving ODEs are certainly more flexible and allow for generalization to unseen ODEs to a certain extent.
- l070-071: repetition "foundation" "foundation"
- l135-136: this sentence is misleading since the Vaswani et al. paper uses an encoder-decoder and not a "decoder-only" architecture
- l247: "baslines"
- l248: "research questions" should exist independently of your paper. However, the three research questions are related to your model
- l322: "higher than", do you mean "lower than"?
- l378: table 3. The title is misleading. You mention "unseen ODEs" while doing fine-tuning, so actually "seeing" these ODEs.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: For clarification, in Sec. 2.1 the data preparation seems to be based upon Euler stepping.  Is this correct?  If so, I find this highly problematics.  One should never use the Euler stepper.... it is highly unstable and very inaccurate.  And adaptive stepping RK45 should be the default.  This detail alone makes me questions the results/performance of the algorithm overall of the paper.  Specifically, I would not trust the results if they are training using Euler.
soundness: 3
strengths: There are some good things to thing to think about here.  I think everyone is interested in good foundation models for dynamical systems.  So the authors have proposed some interesting designs to here build upon.  From that point, I think the paper might be acceptable to ICLR.
summary: The authors propose a foundation model based upon multi-modal information for time-stepping of dynamical systems.
weaknesses: Although the authors claim to "generalize" the method on to new ODES:

"These test ODEs are qualitatively different from the training data, thus convincing to validate FMint’s capability as a foundation model for dynamical systems."

The fact is, these new ODEs have very much the same qualitative behavior as those on which it was trained.  This statement is just false.  Specifically, the "driven-damped pendulum (2D), falling object with air resistence (2D), FitzHugh-Nagumo systems (2D), Pendulum under gravity (2D), and the Ro ̈ssler dynamics (3D)." are very similar to the training data.  


A big concern for me:  the testing on training data or against the speed of ODE solvers seems false.  ONE WOULD ALMOST NEVER USE A FIXED TIME STEP!  Instead, everyone uses adaptive stepping which is often included in adaptive stepping RK45.  To fix the time stepper is a completely unfair way to compare to traditional time-steppers.  This raises serious concerns about the comparisons and actual performance claims for this model.

Also see the question below:  I'm not convinced on the performance frankly as I think the generalization data is very similar.  Moreover, I'm not sure they trained in a fair manner.

Also, there is no comparison to ODEFormer:  https://arxiv.org/abs/2310.05573
Which is also a foundation model in this direction.  

Also, there is no comparison to broader time-series foundation models or things like S4, Mamba etc.  These also can be used for characterizing the dynamics of these systems which are efficient.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: How to prove that the error correction term can effectively improve the limitations of traditional numerical methods?
soundness: 3
strengths: 1. Based on transformer architecture and using context learning method
2. The accuracy is obviously improved, which is 1-2 orders of magnitude higher than the existing numerical method
3. High computational efficiency, 5 times faster than traditional numerical methods
summary: This paper proposes a new method for fast simulation of dynamical systems, FMint, which is based on the decoder-only transformer architecture. The framework adopts the context learning method, uses numerical and text data to learn the general error correction scheme of dynamical systems, and uses the prompt sequence of rough solutions from traditional solvers. The method is pre-trained on 400K ODEs. Their results show that FMint improves accuracy by 1-2 orders of magnitude and computation speed compared to traditional numerical solvers.
weaknesses: 1. Lack of error bound analysis
2. Lack of discussion of long-term predictions for different types of systems