Use Difference operator for discrete system

[sharanry]

#1088 (comment)

TODO:

• Add utilities which enable handling Difference operators.
• Add difference function generators
• Add tests

[sharanry]

Current state of discrete function generator:

User side:

@parameters t σ ρ β
D = Difference(t; dt=0.01)
@variables x(t) y(t) z(t)

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

de = DiscreteSystem(eqs,t,[x,y,z],[σ,ρ,β])

Function produced:

function (ˍ₋out, ˍ₋arg1, ˍ₋arg2, t)
    let var"x(t)" =  @inbounds(ˍ₋arg1[1]), 
        var"y(t)" =  @inbounds(ˍ₋arg1[2]), 
        var"z(t)" =  @inbounds(ˍ₋arg1[3]), 
        σ =  @inbounds(ˍ₋arg2[1]), 
        ρ =  @inbounds(ˍ₋arg2[2]), 
        β =  @inbounds(ˍ₋arg2[3])
        @inbounds begin
            ˍ₋out[1] = (+)(var"x(t)", (*)(σ, (+)(var"y(t)", (*)(-1, var"x(t)"))))
            ˍ₋out[2] = (*)(var"x(t)", (+)(ρ, (*)(-1, var"z(t)")))
            ˍ₋out[3] = (+)(var"z(t)", (*)(var"x(t)", var"y(t)"), (*)(-1, β, var"z(t)"))
            nothing
        end
    end
end

[sharanry]

@ChrisRackauckas @YingboMa The simple SIR test now works fine with Difference operator. Is there any other more complex scenarios I can look at and support or is this sufficient for this initial PR?

[sharanry]

Also, I think currently FunctionMap defaults to dt=1 here. I guess we will have to change that?

[sharanry]

One other thing we are assuming right now is that we have a common difference operator on all states with equal timesteps. Do we intend to support multi-timestep processes?

[ChrisRackauckas]

Also, I think currently FunctionMap defaults to dt=1 here. I guess we will have to change that?

Yes.

[YingboMa]

Shouldn't there also be an in-place version as well?

[sharanry]

Yes. It does get generated. But not sure if we have any use for it in discrete systems as even solvers like FunctionMap uses the not-in-place function (not sure what it is called 😅 ).

[YingboMa]

The code copy-pasting seems easily avoidable.

[sharanry]

Yup. I will try and refactor the code.

[sharanry]

@YingboMa @ChrisRackauckas Does the ODESystem assume some kind of ordering of the equations in the differential system? It looks like the generate_function for ODESystem does not take into consideration what the LHS of each equation is. I also noticed this TODO comment.

(This problem is also propagating to DiscreteSystems with/without Difference operator.)

[YingboMa]

LGTM

[YingboMa]

Oh, remember to also generalize the vars! function.

[sharanry]

Discrete system doesn't require it. Will add it in the difference operator in ODE system PR on which I have started a discussion here?

[YingboMa]

Yes, handling vars! will be required for that issue, but let's do that in a future PR.