Thank you for your extensive feedback. We understand that you have several questions about the paper, and we address each of your points in turn below.

Questions 1 and Other Comments Or Suggestions

First and foremost, we apologize for the multiple typos and notational inconsistencies you encountered. We take full responsibility for these and will correct all of them in the final paper. In Figure 5, $u$ and $v$ represent velocity components in $\boldsymbol{u}=(u,v)$ which is the solution of the Navier-Stokes equation. Also, since the notation $k = (n_1, n_2, \dots, n_N)$ was confusing, it is better to use $k= (k_1, k_2, \dots, k_N)$ and $\|\boldsymbol{k}\| _{\text{mix}}=\Pi _{j=1}^d \max \lbrace 1,k _j\rbrace$.

Question 2

Transforming Fourier basis functions to other basis functions is a very mature and classical research field. However, only exploring the Fourier basis is not unfair. There are two reasons: 1) it is not practical to explore all basis functions in one paper. 2) exploring one kind of basis is enough. For example, Fourier Neural Operator only explores Fourier basis in Neural Operator, but after FNO, researchers propose Spectral Neural Operator, Wavelet Neural Operator, and so on.

We understand that you are worrying about other basis functions. We conducted a simple example with Chebyshev basis on 1D heat equations:

$$u_{xx}=\epsilon u_t, \quad (x,t)\in[-1,1]\times[0,T],$$ $$u(0,x)= \frac{e^{-\epsilon x ^ 2 / 12}} {\sqrt{3}},\quad x \in [-1,1]. \quad u(t,1)=u(t,-1)=\frac{e^{-\epsilon / 12}} {\sqrt{3}}.$$

with the solution $u(x,t)=e^{-\frac{\epsilon x ^ 2 } {4 (t + 3)}} / \sqrt{t + 3}$. In our experiments, $T=1.0, \epsilon=10$.

The results are:

Question 3

The $N$ in this experiment is the number of $\sin$ terms in the initial condition of Eq.20. To avoid confusion, it is better to use $N_k$ instead of $N$.

This experiment is conducted to verify the robustness of SINNs. As we stated in Section 4.3, 'For most problems, the structure of coefficients in the spectral domain is simpler than the structure of solutions in the physical domain.', low frequencies can represent most cases in the real world. But how about some specific cases in which the structure of coefficients is also complex?

Furthermore, more experiments on testing high-frequency terms are asked by other reviewers and can answer your concern of 'decay in the higher frequency components of the initial condition'. You can find the experiments in the reviewers of zt1d and 34of.

Application-driven ML

When comparing algorithms, precise relative error measurement is essential. Real-world cases are noisy, and full of unknowns, making it hard to isolate and assess algorithm performance. By using synthetic examples, we can better evaluate our algorithm's core capabilities, laying a solid groundwork for real-world applications. For example, because of its inherent mathematical complexity, the 2-D Taylor-Green vortex flow has been established as a representative benchmark in computational fluid dynamics for validating numerical methods, fincluding finite element penalty–projection method (doi: 10.1016/j.jcp.2006.01.019), PINN (doi: 10.1016/j.jcp.2021.110676), and so on.

Furthermore, the periodic boundary condition is not real, but it can be used to understand the real world's features, for example, coherent vortex evolution (doi: 10.1017/S0022112084001750), intrinsic physical dynamics(doi: 10.1017/S0022112095000012), and so on.

Theoretical Claims

'Is this statement ..': Yes, the conclusion is held by the universal approximation theorem, thanks for the question.

'Perhaps an NTK-style ..': The error analysis is based on a very ideal assumption. NTK analysis is based on the loss function, which is constructed by the residual equation rather than the error equation in both PINNs and SINNs. To our best knowledge, there is no research that reveals that the error and the residual have direct relationships. Furthermore, empirically speaking, a smaller loss function may not guarantee a smaller error (see Error Certification (arXiv:2305.10157)). However, Ref(doi:10.1016/j.jcp.2023.112527) proves that: adding some terms to the error yields an upper bound including the residual. Thus, if we can find or construct a connection between the residual (or other loss functions, for example, Astral (arXiv:2406.02645)) and the error, perhaps we can use NTK analysis.

Eq.26

Thanks for being concerned about the aliasing error. How to dealiase the error lies in how to deal with $\hat{N}$. In a trivial way, one can do inverse FFT back to the physical domain and do those multiplications. Additionally, one can use the 2/3 rule (doi: 10.1175/1520-0469(1971)028<1074:OTEOAI>2.0.CO;2)

Essential References

Thanks, we will give a further discussion on the relevance of those papers.

Thanks for your positive feedback and useful suggestions.

Firstly, we sincerely appreciate your having caught this typo of $\nabla$ which should be $\nabla$ $\cdot$.

Burgers’ equation with low viscosity

We already conducted an experiment on Burgers’ equation with a fairly low viscosity ($\mu = \pi/150 \approx 0.0209$) which already contrains high-frequency content. But we agree that exploring different viscosities is insightful.

Here we consider a more explicit Burgers' equation:

$$u_t = \nu u_{xx} - u u_x, \quad x\in \left[0,2\pi\right], t \in \left[0,T\right]$$ $$u(0,x) =\sin (x)+\sin(5x)$$

which contains a high frequency component and a low frequency component.

We set $\nu=\pi/r$. The results with different $r$ are as follows:

Helmholtz

Additionally, Reviewer zt1d suggests conducting experiments on Helmholtz equations which only contain a single frequency.

Here we consider a 1D Helmholtz equation:

$$u_{xx} + \lambda^2u=0, x\in [0,2\pi],$$

with the boundary condition

$$u(0)=u(2\pi)=1, \quad u_x(0)=u_x(2\pi)=\lambda.$$

The solution is $u=\cos(\lambda x)+\sin( \lambda x)$. The results are as follows:

Thank you for your detailed review. We address each of your points below:

Is there any particular reason why you did not compare with STDE?

STDE is indeed a recent work to speed up high-order derivatives in PINNs. However, STDE works by amortizing the cost of computing mixed partials in multivariate problems through random estimations. Thus, for univariate problems, its “stochastic” part isn’t needed and it essentially performs the same operations as Taylor-mode. In our manuscript, to intuitively and concisely demonstrate the results, we choose a univariate function with only one spatial derivative term. Herein, Taylor-mode performs the same as STDE. As for multivariate functions, because STDE is an estimator and still uses univariate Taylor-mode, STDE still scales with the order of derivatives and the dimensionality. However, SINNs remain in a constant number with both the order of derivatives and the dimensionality.

To provide a more convinced conclusion, we conducted a test on the efficiency of $d$-variate problems. If we show SINN is constant with the dimensionality, combining the result that SINNs is constant with the order of derivative, we can claim that SINNs is more efficient than STDE. We consider a Poisson's equation:

$$\Delta u =f, \boldsymbol{x} \in [0,2\pi]^d$$

with $u=\sum_{i=1}^d \cos(x_i)$ and $f=-\sum_{i=1}^d \cos(x_i)$. The results are as follows:

The memory increases because the input scales with $d$. But, because our operator is only multiplication, the iteration rate is almost constant as expected.

Theoretical Claims

We agree that we are considering an ideal setting that the network can express the truncated Fourier series accurately. And we admit that, in practice, the error we obtained is worse than the theorem results. We will clarify in the paper that the theoretical exponential convergence rate is an asymptotic result assuming sufficient capacity.

Additionally, we don't expect treating the approximation error as exactly zero. $\| u^{\theta^*}-u^*\| \ll \|u^{\theta^*_N}- u^{\theta^*}\|$ is enough in error analysis.

Irregular geometries

Thanks for pointing out the limitation that our experiments were on regular domains. However, we argue that our SINNs can be extended to irregular geometries with additional operations. There are several approaches to deal with this limitation, for example:

1. The moving mesh method can construct a coordinate transformation to map an irregular domain to a regular domain (or map subdomains respectively, for example, doi: 10.1016/j.jcp.2020.109835), and it has already been used in PINNs by PhyGeoNet (Gao et al., 2021). (impressive results in Fig. 5 and 7).

2. One can extend the irregular domain to the regular domain, for example, the fictitious domain method(doi: 10.1137/20M1345153).

3. Even, there is a simpler and more straightforward approach if the irregular region can be expected to be expanded into a regular region: one can just add the irregular geometries as a constraint by inverse Fourier transformation. Mathmatically speaking, for PINNs: we minimize the boundary loss by $\mathcal{L} [u_{NN}] (x_b)$, where $u_{NN}$ is the output; For SINNs: we minimize the boundary loss by $\mathcal{L} [\mathcal{F}^{-1} [\hat{u}_{NN}]] (x_b)$, where $\hat{u} _{NN}$ is the output.

Here we consider a 2D Helmholtz equation on a plate with a hole using the third approach:

$$\Delta u+ \lambda^2u=0, \boldsymbol{x}\in \Omega = [0,2\pi]^2\textbackslash B_\pi(\boldsymbol{x}),$$

with Dirichlet boundary condition:

$$u(\boldsymbol{x})=g(\boldsymbol{x}), \boldsymbol{x} \in \partial \Omega,$$

where $B_\pi(\boldsymbol{x})$ is the ball with the center at the origin (0,0) and a radius of $\pi$ and $g$ is derived by the solution $u=(\cos(\lambda x)+\sin(\lambda x))(\cos(\lambda y)+\sin(\lambda y))$.

Here are the results:

Thanks for the positive assessment and thoughtful questions.

Challenging Problems

Firstly, we consider a well-defined 1D Helmholtz equation as you suggested:

$$u_{xx} + \lambda^2u=0, x\in [0,2\pi],$$

with the boundary condition

$$u(0)=u(2\pi)=1, \quad u_x(0)=u_x(2\pi)=\lambda.$$

The solution is $u=\cos(\lambda x)+\sin( \lambda x)$. The results are as follows:

For a 2D Helmholtz equation:

$$\Delta u+ 8u=0, \boldsymbol{x}\in \Omega = [0,2\pi]^2,$$

with the condition:

$$u(\boldsymbol{x})=g(\boldsymbol{x}), \boldsymbol{x} \in \partial \Omega,$$ $$u_x(\boldsymbol{x})=p(\boldsymbol{x}), \boldsymbol{x} \in \lbrace 0 \rbrace \times[0,2\pi],$$ $$u_y(\boldsymbol{x})=q(\boldsymbol{x}), \boldsymbol{x} \in [0,2\pi] \times \lbrace 0 \rbrace.$$

$g,p,q$ are derived by $u=(\cos(2x)+\sin(2x))*(\cos(2y)+\sin(2y))$.

The results are:

Furthermore, we don't think this problem is a challenge for SINNs. Because solving the boundary-value problems are like optimizing:

$$\min_{x} |Ax-b|^2.$$

Thus, changing $\lambda$ is merely changing the location of the non-zero element in $x$ in SINNs.

Irregular geometries

Due to the limitations, please read it in reviewer PNMU.

High-dimensional problems

Thanks for pointing out this limitation. There is a difference between using SINNs to solve high-dimensional PDEs and low-dimensional PDEs: PINNs can handle high-dimensional problems because their error is dependent on $N$ which is the number of sampled points and independent of the dimensionality $d$. As SINNs sample $N$ frequencies, if the error is also independent of the dimensionality $d$, SINNs can address the curse of dimensionality (CoD). Fortunately, by utilizing the optimized hyperbolic cross, (Shen et al., 2011) proved that the error is also independent of the dimensionality $d$ (Key result is Corollary 8.3). Some examples can be found in (Shen & Yu, 2010). Herein, we will apply high-dim Fourier operators by a sparse matrix to avoid CoD.

And we show the scaling by an ideal example:

$$\Delta u =f, \boldsymbol{x} \in [0,2\pi]^d$$

with $u=\sum_{i=1}^d \cos(x_i)$ and $f=-\sum_{i=1}^d \cos(x_i)$. Here are the results (memory and rate are for SINN):

Theoretical Claims

The sampling error exists in PINNs because the Monte Carlo method is used to estimate the integration. For example, the residual term:

$$\min_{\theta} \int_{\Omega}|\mathcal{L}_r(x;\theta)|^2 d x \approx \min _{\theta}\mathcal{L}= \min _{\theta} \sum _{i=0}^M \left|\mathcal{L}_r\left(x_i;\theta\right)\right|^2$$

As the integration domain contains infinite points, using finite points to estimate it always has an error, i.e. the sampling error in Monte Carlo methods. However,SINNs don't need Monte Carlo methods because we don't have integrations and our input domain is finite: the frequencies. Herein, our loss function is directly:

$$\min_{\theta} \tilde{\mathcal{L}}= \min _{\theta} \sum _{i=0}^M|\tilde{\mathcal{L}}_r(k_i;\theta)|^2$$

With the assumption that the network is powerful enough (i.e. regardless of the sampling method selected, finally, the network learns the total dataset). Then $\tilde{\mathcal{L}}=0$ means the output of SINNs exactly equals to the coefficients of the frequencies. But $\mathcal{L}=0$ only means the output of PINNs is accurate at the points $\lbrace x_i \rbrace_{i=1}^M$ but still has the sampling error from Monte Carlo methods.

In conclusion, because we truncate the Fourier series to a finite number, SINNs have the spectral error but don't have the sampling error.

Essential References

Thanks for the given relevant papers. However, the frameworks of PINNs and Neural Operators are different. PINNs aim to learn the map from a definition domain to a functional space while NOs learn the map from a functional space to a functional space. Herein, we didn't write the review of Neural Operators in our manuscript. But, indeed, we should write a brief review of neural networks with spectral methods in other fields.