Please suggest how to extend this to include (u^2){tt} and u{3t}. This problem occurs in non-linear acoustics when solving the Westervelt Equation. See e.g. Nonlinear acoustics - Wikipedia .

The Westervelt equation looks like a partial differential equation so you would need to discretize the spatial dimensions somehow (finite elements, finite volume, spectral methods etc).

I see more clearly now, what you mean. After you discretize the spatial dimensions you arrive at a third order nonlinear ODE problem. You can rewrite this as a nonlinear first order problem by including the first and second order time derivative in your state vector. (I admit this is not straight forward but I don’t see any technical limitations and the end result should be an explicit equation of the form du=f(u,p,t)dt which you can solve as usual.

Note, \partial_t^2p^2 can be exactly rewritten as 2(\partial_t p)^2 + 2p\partial_t^2p.

The challenge, however, I see is that after rearranging terms, we obtain

(1 + 2p) \partial_t^2 p = f(p_t, p, t).

The (1+2p) factor in the LHS bring the equation in non-standard/non-canonical form (in which we have 1 \partial_t^2 p = f() ). I wonder how to treat this non-canonical form. I imaging some coding to be required.

\partial_t
\begin{pmatrix}
p\\\dot p\\ \ddot p
\end{pmatrix}
=
\begin{pmatrix}
\dot p\\
\ddot p\\
\frac{c_0^4}{\delta} \left(-\nabla^2 p + c_0^{-2} \ddot p - \frac{2\beta }{\rho_0c_0^4} \left[ \dot p^2 + p\ddot p \right]\right)
\end{pmatrix}

So the state vector is (p, \dot p, \ddot p) and the right hand side is an explicit function of those quantities. (The implementation of \nabla^2 follows from your spatial discretization.)

I would also note that most terms on the r.h.s. are linear, including the Laplacian term which can be quite stiff. If you can somehow diagonalize the linear operator you can use a SplitODEProblem in DifferentialEquations.jl which should improve performance significantly.

\partial_t
\begin{pmatrix}
p\\\dot p\\ \ddot p
\end{pmatrix}
=
\begin{pmatrix}
0&1&0\\
0&0&1\\
-\frac{c_0^4}{\delta}\nabla^2&0&\frac{c_0^2}{\delta}\\
\end{pmatrix}
\begin{pmatrix}
p\\\dot p\\ \ddot p
\end{pmatrix}
+\begin{pmatrix}
0\\
0\\
-\frac{2\beta}{\rho_0\delta}\left[ \dot p^2 + p\ddot p \right]
\end{pmatrix}

Is my understanding correct that SplitODEProblem would allow implicit treatment of the term including the Laplacian (first term in your notation) and explicit treatment of the other term (second term in your notation)?

Yes, exactly. If you find a way to calculate the exponential of the linear operator efficiently (for example, by using its sparse structure), you can achieve better performance.

Hello, thanks for all the responses. I have a question about calculating the exponential of the linear operator. I assume its related to this: Matrix exponential. But i dont see how i would use this to improve the performance of a splitODE solver. Do you have some more information regarding this?

If you manage to diagonalize the matrix, then its matrix exponential is very fast to compute because you just exponentiate each entry on the diagonal. With a split ode solver you could, for example, use the Diagonal type from the LinearAlgebra module to achieve this.

The SplitODE solvers in DifferentialEquations.jl take advantage of that.

Furthermore, the SplitODE solvers can give good speed ups when the largest rate of change comes from the linear operator (because you can increase your timestep without sacrificing accuracy).

Matrix diagonalization and matrix exponential is likely to be too expensive. Would using algebraic multigrid on the diffusion part on the right-hand side make sense to apply here?