SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
556 stars 210 forks source link

Allow using LinearAlgebra.I(n) as is #2090

Open homocomputeris opened 10 months ago

homocomputeris commented 10 months ago

Is your feature request related to a problem? Please describe.

When passed to a problem constructor, LinearAlgebra.I(n) should be treated as is, without manual type conversion by the user. For now, it requires some boilerplate to avoid errors that appear in various usages, e.g.:

using DifferentialEquations
using LinearAlgebra

function f(u, p, t)
    return [0.0 1.0;
        -1.0 0.0] * u
end

tspan = (0.0, pi)
h = pi

# u0 = LinearAlgebra.I
# # MethodError: no method matching length(::UniformScaling{Bool})

# u0 = LinearAlgebra.I(2)
# # InexactError: Bool(1//1000000)

# u0 = float(LinearAlgebra.I(2))
# # ERROR: MethodError: Cannot `convert` an object of type Matrix{Float64} to an object of type Diagonal{Float64, Vector{Float64}}

u0 = float(Matrix(LinearAlgebra.I(2)))
# this works but not obvious why the (kinda) automatic and universal `I` requires this boilerplate, especially given the fact than `eye` is deprecated

prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Euler(), dt=h).u[end]

I've seen this issue only with IVs but it might appear in other cases.

Describe the solution you’d like

Probably, problem constructors or solvers should just typecast to an appropriate type. Probably, always casting to float(Matrix(LinearAlgebra.I(2))) would be safe.

Describe alternatives you’ve considered

Copy-pasting the boilerplate.

ChrisRackauckas commented 10 months ago

Is there a reason you cannot use I inside of the f? Functionally there's effectively no difference because I*x doesn't allocate, but it would then make u0 just a vector.

homocomputeris commented 10 months ago

Passing the identity matrix as the IVs is convenient to get the solution in form of the corresponding fundamental matrices. Or maybe there is another way to get the solution in this form?

ChrisRackauckas commented 10 months ago

2*I though is a very specific object though. It's not a matrix or a number, and it doesn't conform to any standard interfaces like that. If you were to keep this form then it also would only allow multiples of the identity matrix as the solution, which wouldn't hold in almost any case. In the cases which it did hold, you could just make it a scalar valued differential equation.

homocomputeris commented 10 months ago

it doesn't conform to any standard interfaces like that.

This is why I suggest it may be type-cast to a matrix at problem construction similar to this https://github.com/SciML/Integrals.jl/issues/208.

But maybe passing I(n) to a problem is a very rare use case that doesn't require changes.