Usage of OrdinaryDiffEqLinear.jl for non-autonomous problems

General Usage
, ,
1

I’m trying to use OrdinaryDiffEqLinear.jl 's the state-independent Magnus solvers to solve non-autonomous Linear ODEs.

To confirm the usage, I’m trying to solve

\dot{\boldsymbol{u}}(t)=\begin{pmatrix}0&-\alpha t\\\alpha t&0\end{pmatrix}\boldsymbol{u}(t),

where \alpha is a constant (fixed as 0.1).

Its exact solution may be given by

\boldsymbol{u}(t)=\begin{pmatrix}\cos\left(\frac{\alpha t^2}{2}\right)& -\sin\left(\frac{\alpha t^2}{2}\right)\\ \sin\left(\frac{\alpha t^2}{2}\right)& \cos\left(\frac{\alpha t^2}{2}\right)\end{pmatrix}\boldsymbol{u}(0).

I cannot reproduce the exact solution.

Here is my code:

using OrdinaryDiffEq
using Plots

alpha = 0.1; t0 = 0.0; t1 = 1.0; dt = 0.01
tspan = (t0, t1) 
u0 = ones(2) # initial condition
A0 = [0.0 1.0; 0.0 0.0] # time-independent part

function update_func(A, u, p, t)
    A[1, 2] = -alpha * t 
    A[2, 1] = alpha * t
end

# Magnus solver
A = MatrixOperator(A0; update_func)
prob = ODEProblem(A, u0, tspan)
sol = solve(prob, MagnusGL8(); dt)
plot(sol, label=["u[1] (magnus)" "u[2] (magnus)"])

# exact solutions
F(t0,t1) = alpha * (t1^2 - t0^2) / 2
u_exact1(t) = cos(F(t0, t)) * u0[1] - sin(F(t0, t)) * u0[2]
u_exact2(t) = sin(F(t0, t)) * u0[1] + cos(F(t0, t)) * u0[2]
plot!(sol.t, u_exact1.(sol.t), label="u[1] (exact)")
plot!(sol.t, u_exact2.(sol.t), label="u[2] (exact)")

# autonomous solver
prob = ODEProblem(MatrixOperator(A0), u0, tspan)
sol = solve(prob, LinearExponential())
plot!(sol, label=["u[1] (exponential)" "u[2] (exponential)"], color=[:lightblue :orange], alpha=0.8, linestyle=:dash, lw=3)

Here are the obtained results.

  • magnus is obtained by the Magnus solver (here is MagnusGL8()).
  • exact is obtained by hand. The formula is shown above.
  • exponential is obtained by exponentiation. This is a solver for autonomous cases.

result
result1192×784 40.9 KB

Although magnus and exact are expected to behave similarly, it’s actually magnus and exponential that do.
This suggests that update_func may not be passed correctly to the solver.

If I’m misusing it, please let me know how to pass it correctly.

Julia v1.11.4
OrdinaryDiffEqLinear v1.1.0

2

That looks right on your end. Open an issue.

3

Thank you for your quick reply. I’ll open it.

4

Opened Time dependence in a non-autonomous linear ODE problem does not appear to propagate correctly to the solver · Issue #2661 · SciML/OrdinaryDiffEq.jl · GitHub

5

Are you using MatrixOperator correctly? There is update_func and update_func! and the latter is the one you should be using I think.

using OrdinaryDiffEq
using Plots

alpha = 0.1; t0 = 0.0; t1 = 10.0; dt = 0.01
tspan = (t0, t1) 
u0 = ones(2) # initial condition
A0 = [0.0 1.0; 0.0 0.0] # time-independent part

function update_func!(A, u, p, t)
    A[1, 2] = -alpha * t 
    A[2, 1] = alpha * t
end

# Magnus solver
A = MatrixOperator(A0; update_func!)
prob = ODEProblem(A, u0, tspan)
sol = solve(prob, MagnusGL8(); dt)
plot(sol, label=["u[1] (magnus)" "u[2] (magnus)"], seriestype=:line, lw=7, opacity=0.3)

# exact solutions
F(t0,t1) = alpha * (t1^2 - t0^2) / 2
u_exact1(t) = cos(F(t0, t)) * u0[1] - sin(F(t0, t)) * u0[2]
u_exact2(t) = sin(F(t0, t)) * u0[1] + cos(F(t0, t)) * u0[2]
plot!(sol.t, u_exact1.(sol.t), label="u[1] (exact)")
plot!(sol.t, u_exact2.(sol.t), label="u[2] (exact)")

image
image600×400 33.5 KB

EDIT: Ah I just saw that you also wrote this in the linked issue. I’ll just leave this here for posterity.

1 Like
6

I appreciate your time and effort in resolving this issue.

1 Like