# Solving a Diagonal ODE using DifferentialEquations.jl
**URL:** https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956
**Category:** General Usage
**Tags:** ode, differentialequation
**Created:** [October 5, 2024, 11:00am UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956 "2024-10-05T11:00:45Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)
#### Post date: [October 5, 2024, 11:00am UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/1 "2024-10-05T11:00:45Z")
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.
---
### Author:  [@dawbarton](https://discourse.julialang.org/u/dawbarton)
#### Post date: [October 5, 2024, 11:18am UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/2 "2024-10-05T11:18:48Z")
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?)
---
### Author:  [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)
#### Post date: [October 5, 2024, 11:25am UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/3 "2024-10-05T11:25:16Z")
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.
---
### Author:  [@dawbarton](https://discourse.julialang.org/u/dawbarton)
#### Post date: [October 5, 2024, 12:23pm UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/4 "2024-10-05T12:23:45Z")
Something like
```julia
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
```
---
### Author:  [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)
#### Post date: [October 5, 2024, 1:16pm UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/5 "2024-10-05T13:16:00Z")
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.
---
### Author:  [@dawbarton](https://discourse.julialang.org/u/dawbarton)
#### Post date: [October 5, 2024, 2:38pm UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/6 "2024-10-05T14:38:14Z")
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.
---
### Author:  [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)
#### Post date: [October 5, 2024, 4:55pm UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/7 "2024-10-05T16:55:46Z")
Thanks a lot!
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [October 5, 2024, 7:32pm UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/8 "2024-10-05T19:32:42Z")
> [@albertomercurio](#):
>
> I don’t need to construct the Krylov subspace
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-repl
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.
---
### Author:  [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)
#### Post date: [October 6, 2024, 11:43am UTC](https://discourse.julialang.org/t/solving-a-diagonal-ode-using-differentialequations-jl/120956/9 "2024-10-06T11:43:21Z")
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)`?