# How to extend SecondOrderODEProblem in DifferentialEquations.jl?

**URL:** <https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831>\
**Category:** Numerics\
**Tags:** package, differentialequation\
**Created:** [December 30, 2021, 6:57pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831 "2021-12-30T18:57:44Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [December 30, 2021, 6:57pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/1 "2021-12-30T18:57:44Z")

</div>

Greetings.

SecondOrderODEProblem in DifferentialEquations solves u\_{tt} = f(u\_t,u,p,t), see [Dynamical, Hamiltonian and 2nd Order ODE Problems · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/types/dynamical_types/#Mathematical-Specification-of-a-Dynamical-ODE-Problem) .

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](https://en.wikipedia.org/wiki/Nonlinear_acoustics) .

Thanks, Domenico.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [December 30, 2021, 7:00pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/2 "2021-12-30T19:00:37Z")

</div>

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](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757)

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [December 30, 2021, 7:09pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/3 "2021-12-30T19:09:42Z")

</div>

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](https://ranocha.de/SummationByPartsOperators.jl/stable/tutorials/wave_equation/) and for the Korteweg-deVries equation for instance here [Solving KdV Solitons with Upwinding Operators · DiffEqOperators.jl](https://diffeqoperators.sciml.ai/stable/operator_tutorials/kdv/)

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

Thx, Domenico.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [December 30, 2021, 7:42pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/4 "2021-12-30T19:42:41Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [December 30, 2021, 8:03pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/5 "2021-12-30T20:03:59Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [December 30, 2021, 9:10pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/6 "2021-12-30T21:10:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [December 30, 2021, 9:17pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/7 "2021-12-30T21:17:33Z")

</div>

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

Thank you. Domenico.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [December 31, 2021, 11:08am UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/8 "2021-12-31T11:08:52Z")

</div>

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}

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [December 31, 2021, 12:45pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/9 "2021-12-31T12:45:04Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [December 31, 2021, 12:53pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/10 "2021-12-31T12:53:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [January 2, 2022, 6:18pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/11 "2022-01-02T18:18:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Siempie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/siempie/32/32790_2.png) [@Siempie](https://discourse.julialang.org/u/Siempie)\
**Post date:** [February 8, 2022, 12:15pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/12 "2022-02-08T12:15:59Z")

</div>

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](https://en.wikipedia.org/wiki/Matrix_exponential#Linear_differential_equation_systems). 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?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [February 8, 2022, 12:35pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/13 "2022-02-08T12:35:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [February 27, 2022, 9:43pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/14 "2022-02-27T21:43:05Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 1, 2022, 12:55pm UTC](https://discourse.julialang.org/t/how-to-extend-secondorderodeproblem-in-differentialequations-jl/73831/15 "2022-03-01T12:55:52Z")

</div>

It could, yes. See [Solving Large Stiff Equations · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/)