Solving a Diagonal ODE using DifferentialEquations.jl

General Usage
,
1

Hello,

I have the following ODE system

\dot{\mathbf{u}} = \mathbf{A} \ \mathbf{u}

where \mathbf{u} is a vector and \mathbf{A} is a Diagonal matrix. The following ODE system admits an analytical solution of the form

\mathbf{u} (t) = \exp [\mathrm{diag}(\mathbf{A}) \ t ] \cdot \mathbf{u}_0

so I wouldn’t need DifferentialEquations.jl to solve it. However, I would take advantage of its OrdinaryDiffEqCallbacks.jl framework to implement for example a ContinuousCallback with interpolation.

I was wondering if there exists a method to implement this. I tried the LinearExponential algorithm, but this is an overkill for my problem, because I don’t need to construct the Krylov subspace, since I already have the analytical solution. I was thinking of creating a custom solver like DiagonalExponential and defining only the methods needed from the integrator interface, but I don’t know how to proceed exactly.

2

Since you have the functional form of the solution, is there any reason you wouldn’t use NonlinearSolve.jl directly with your callback? (Or SimpleNonlinearSolve.jl?)

3

How exactly? I want to trigger a callback when the norm of the vector \mathbf{u} (t) is smaller than a random number r, then I basically change the value of \mathbf{u} and I keep evolving the state analytically until the next jump.

4

Something like

using SimpleNonlinearSolve
using LinearAlgebra: norm

u0 = [1.0, 2.0, 3.0]  # ICs
λ = [-1.0, -1.5, -2.0]  # Eigenvalues / diagonals
t = 0.0  # Start time
r = 0.5  # Required norm

f(t, p) = norm(exp.(p.λ .* t) .* p.u0) - p.r

prob = NonlinearProblem{false}(f, t, (; λ, u0, r))
sol = solve(prob, SimpleNewtonRaphson())
if SimpleNonlinearSolve.SciMLBase.successful_retcode(sol)
    t = sol.u
    println("Converged to time: $t")
    # Do callback
else
    println("Failed to converge")
end
1 Like
5

Thank you very much for this example. It seems in the direction I want.

Then I only need to repeat this procedure up to a predefined time t_f, right?

Do you think that auto differentiation using a loss function similar to your f^2 (t, p) would be more efficient? In principle it should take advantage of auto differentiation, finding the minimum of the function.

6

Yes, just iterate until you hit the required time.

I would be very surprised if formulating it as a loss function would be helpful. Generally speaking if the problem is a root finding problem (this is), using an optimiser on it will be very inefficient. The solver here will already be using auto diff to calculate the Jacobian. (I’m guessing that AutoForwardDiff is the default.)

The only way to get more performance out of something like this is to use StaticArrays (if the problem is reasonably small) or change the function to be in place rather than allocating.

1 Like
7

Thanks a lot!

8

You can pass krylov = :off, but you don’t have to—it’s the default. Adapting the example from the docs to a diagonal operator:

julia> using LinearAlgebra, OrdinaryDiffEq, SciMLOperators

julia> A = MatrixOperator(Diagonal([-0.2, -0.35]))
DiagonalOperator(2 × 2)

julia> prob = ODEProblem(A, [1.0, -1.0], (1.0, 6.0));

julia> sol = solve(prob, LinearExponential())
retcode: Success
Interpolation: 3rd order Hermite
t: 2-element Vector{Float64}:
 1.0
 6.0
u: 2-element Vector{Vector{Float64}}:
 [1.0, -1.0]
 [0.36787944117144217, -0.17377394345044517]

Add your CotinuousCallback and you’re good to go.

But of course, @dawbarton’s solution works just as well and is arguably more transparent while having fewer and smaller dependencies.

1 Like
9

Oh, thats true. In this way we directly apply the exp to a Diagonal matrix. Is this action non-allocating? I mean does it allocates anytime exp(A) u, or does it compute mul!(du, exp(A), u)?