I am interested in using a numerical wave tank (NWT) to study the performance of various water wave dampers using MATLAB. I am looking at nonlinear water waves. My current NWT (2d, periodic) is based on the Zakharov/Craig-Sulem formulation, where the Dirichlet-Neumann map is represented via Cauchy integrals and spatial derivatives are computed spectrally via FFT. I'm currently working on setting up some other NWTs using other formulations, but that's a bit beside the point I suppose.

For the time integrator, I've mostly been using MATLAB's ode45. I have coded up a couple other very basic Runge-Kutta integrators, but they are currently very slow (e.g., they use fixed time steps). This leads to my question. As one component of evaluating the performance of various dampers, I am interested in looking at dissipation of energy. As such, I would like the time integrator to come as close as possible to conserving the energy. If we let $E = E(t)$ denote the energy, my current scheme has errors on the order of $E(t) - E(0) \sim 10^{-6}$ after approximately $1000$ time steps. I would really like to reduce this as much as possible. Theoretically, a symplectic integrator (e.g., the implicit midpoint method) would seem like a good choice, but I've read some not great things about them in the context of water wave simulations. However, perhaps they would indeed be a good choice for my rather narrow purpose?

Anyway, my question is: could anyone recommend a time integrator given the above objective(s)? Of course, I would like it to be as fast as possible, but I will give a greater weight to reducing the energy drift. Other comments or ideas are very much welcomed as well. I will give one final caveat: I have a lot more experience in the theoretical analysis of water waves and am a relative novice when it comes to the numerical side. Perhaps this is obvious from my posting, but I wanted to put that out there.

Edit: This edit is to address the (first) comment below.

(1) Yes, the energy is conserved in the exact solution.

(2) Here my relative lack of experience with numerics will show: I'm not entirely sure of the answer off the top of my head. My suspicion is that the energy would not necessarily be conserved up to rounding error by the spatial discretization, but the energy drift from the spatial discretization would be much smaller than that from the timestepping (which is why I've been focused on trying to improve my time integrator). However, I could be wrong here.

(3) The Zakharov formulation is a Hamiltonian system with canonical variables $\tilde{\phi}$ and $\eta$, where $\tilde{\phi}$ is the trace of the velocity potential $\phi$ along the free surface and the free surface is given by the graph of $\eta$. The Hamiltonian energy is then given, in my case, by $$E = \frac{1}{2}\int_0^{2\pi} \tilde{\phi}G(\eta)\tilde{\phi} \ dx + \frac{g}{2}\int_0^{2\pi} \eta^2 \ dx + \tau\int_0^{2\pi} \sqrt{1 + (\partial_x\eta)^2} - 1 \ dx.$$ In the above, $G(\eta)$ is the (normalized) Dirichlet-Neumann map associated to the fluid domain (i.e., $G(\eta)\tilde{\phi} = \sqrt{1 + (\partial_x \eta)^2}\partial_n\phi\rvert_{y = \eta}$), $g$ is acceleration due to gravity and $\tau$ is the coefficient of surface tension. So, the first term is the kinetic energy and is equal to $\frac12\iint_{\Omega_t} |\nabla \phi|^2 \ dydx$, where $\Omega_t$ is the fluid domain at time $t$. The second term is gravitational potential energy and the third term is potential energy due to stretching of the free surface.

  • $\begingroup$ It would help if you could edit your answer to include the following: 1. Is energy conserved in the exact solution of the PDE you are discretizing? 2. Is energy conserved (up to rounding errors, not just truncation errors) by the spatial discretization you are using? 3. What is the expression for the energy in terms of the variables used in your numerical simulation? $\endgroup$ May 18, 2022 at 3:21

1 Answer 1


Time discretizations that conserve energy are designed to be applied to ODE systems that have a conserved energy. In your case, it seems that the ODE system (i.e., the spatial discretization of your PDE) is not energy-conserving. I think your first step should be to understand how dissipative your spatial discretization is. You can do this by fixing the spatial grid and running with a sequence of smaller and smaller time steps. You should see that the loss of energy at the final time converges to some (non-zero) value. Now you know that you can't expect to do better than that value by changing the time discretization.

For linear problems, it's possible to design time discretizations that "compensate" for dissipation in the space discretization, so that the full discretization has little or no dissipation. But for nonlinear problems I don't think anyone knows a reliable way to do this.

If you decide that energy conservation is important enough, you could try to develop an energy-conserving spatial discretization, for instance by using summation-by-parts techniques. If you then wanted to conserve energy in time while using explicit methods, I recommend relaxation Runge-Kutta methods (disclosure: these methods were developed by myself and my collaborators):

The basic idea is simply to adjust the size of your Runge-Kutta step (not the time step, but the solution increment) in order to satisfy conservation. In your case, this will require the solution of a scalar nonlinear algebraic equation at each step.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.