How to extend SecondOrderODEProblem in DifferentialEquations.jl?

Specific Domains Numerics
,
1

Greetings.

SecondOrderODEProblem in DifferentialEquations solves u_{tt} = f(u_t,u,p,t), see Dynamical, Hamiltonian and 2nd Order ODE Problems · DifferentialEquations.jl .

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 .

Thanks, Domenico.

2

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).

Also Please read: make it easier to help you

3

Yes, agreed, true.

We will first discretize in space and subsequently solve the semi-discrete system. The approach is referred to as method of lines by some.

The approach is documented for linear wave for instance here Wave equation · SummationByPartsOperators.jl and for the Korteweg-deVries equation for instance here Solving KdV Solitons with Upwinding Operators · DiffEqOperators.jl

Our wish is to adapt these examples to our needs. Does this clarify your doubt?

Thx, Domenico.

4

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.

5

Sincere thanks again for your valuable input.

I do agree in part.

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.

How do you view this matter?

6

I would write the problem like so:

\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.)

1 Like
7

Aha! Interesting. Allow me to think about this and get back to later.

Thank you. Domenico.

8

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}
2 Likes
9

Sincere thanks again.

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)?

10

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.

11

I arrived at a prototype implementation and asked my collaborator to give it a look. More later.

12

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?

13

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).

14

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?

15

It could, yes. See Solving Large Stiff Equations · DifferentialEquations.jl