Solving difference equations

[seadra]

Is it possible to solve a difference equation with this package?

Difference equations are mentioned in the documentation and README, so it seems this is possible, but it's not clear how I can do this and I can't find any examples anywhere.

[ChrisRackauckas]

Let me know if it has a discoverability issue.

https://diffeq.sciml.ai/stable/types/discrete_types/ https://diffeq.sciml.ai/stable/solvers/discrete_solve/

[seadra]

Thanks for the response! Ah, I see. Those pages don't refer to them as "difference equation", which is why I failed to find by using search. Adding "difference equation" would sure help with discoverability.

Is it limited to 1st order discrete difference equations? I'm actually looking to solve a continuous difference equation

u(x) = A(x)u(x+a) + B(x)u(-x+b)

which can be solved after discretizing the x-axis, maybe with a dx = a/100 or b/100, in order to get a a roughly OK solution. However, this would mean solving

U(n) = A(n dx)U(n+100) + B(n dx)U(n+200)

(taking 2a=b as an example). If I'm not missing anything, solving this with the discrete solver means creating a 200-dimensional coupled ODE problem.

[ChrisRackauckas]

Does scale_by_time do waht you're looking for?

https://diffeq.sciml.ai/stable/solvers/discrete_solve/#OrdinaryDiffEq.jl

[seadra]

Hi Chris,

So I've been pondering about this, but I just couldn't see how scale_by_time can help me find a solution to this continuous equation:

u(t) = A(t)u(t+a) + B(t)u(-t+b)

(Note the sign of the argument for the second RHS term is negative) It doesn't seem to be possible to access u at different times when specifying the problem.

[seadra]

Or is it possible to somehow use DiffEqFlux to solve an equation like u(t) = A(t)u(t+a) + B(t)u(-t+b)?

It's like a delay differential equation, but alas, the LHS is not the derive but the function itself :(

[ChrisRackauckas]

It uses points from the future? Sounds like you need to solve that as a fully implicit system via something like NLsolve.jl

[seadra]

I ended up discretizing the problem in t and using IterativeSolvers.jl. Getting a reasonably smooth end result turned out to be quite inefficient though, and required ~100000 points. I was hoping DifferentEquations.jl might provide an adaptive solver, or at least help me discretize the problem (like ModelingToolkit.jl --I suppose I can't use it directly?), but maybe they won't work for such continuous finite difference equations.

Thanks for mentioning NLsolve.jl, I checked it but if I'm not missing anything, it outputs a number/vector (the x satisfying F(x)=0), whereas in my problem the output is a function. Perhaps you mean I can use nlsolve after manually discretizing the problem?

[ChrisRackauckas]

Yes, you'd have to discretize (unless you solve it as a neural network or something like that).

[seadra]

That's an interesting direction, I'd like to give it a try with Flux!

With regards to discretization though, I'm wondering if ModelingToolkit.jl can help. Right now, I'm manually generating the matrix that encodes the equation u(t) = A(t)u(t+a) + B(t)u(-t+b) + C(t)

M u = v

where u contains the discrete points u(t_i), and v encodes the "source" term C(t_i).

Would it be possible to obtain the (sparse) M matrix from the equation u(t) = A(t)u(t+a) + B(t)u(-t+b) + C(t) using ModelingToolkit.jl? I know it's not exactly a PDE, but I believe it discretized by the same method.

The examples on ModelingToolkit.jl don't seem to allow inclusion of the time parameter in the equation, but maybe it's somehow possible?

[ChrisRackauckas]

Yes, and the Runge-Kutta discretizer in https://github.com/SciML/ModelingToolkit.jl/pull/562 is a somewhat more complicated proof of concept of that. You'd create the array of time points and then churn through it building the expressions, and then calculate the Jacobian to get the M.

The examples on ModelingToolkit.jl don't seem to allow inclusion of the time parameter in the equation, but maybe it's somehow possible?

What do you mean?

[seadra]

Sorry for the lack of clarity. What I mean is, the way we express an equation with ModelingToolkit is like this:

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

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

In the equations, it seems time is implicit and we can't use the unknown function (x,y,z in this case, or u in my case) at two different times. So something like this wouldn't work, right?

0 ~ x(t) + x(-t+10)

[ChrisRackauckas]

You can do @variables x(..) and then do 0 ~ x(t) + x(-t+10). Doing @variables x(t) is what defaults it to mean x(t) whenever you say x. The .. leaves it as an open function.

[seadra]

Oh, wow, thanks. That's so nice!

I'll check the documents to figure out how I can extract the matrix such that I can obtain the matrix M for the discretized problem M x = v corresponding to that equation

[ChrisRackauckas]

You can probably just use the Jacobian to extract the matrix form.

[seadra]

For example, if we limit t to be from 0 to 100, and use 10000 grid points Δt = 100/10000. Then sparse matrix corresponding to the equation u(t) + u(t+10) == 0 would be M = δ_{i, i+10/Δt} where δ_{i,j} is Kronecker delta, such that M*[u_1 ... u_10000] = 0 holds.

I tried using Jacobian like this

@variables t
@variables u(..)

eqs = [0 ~ u(t) + u(t+10)]

ns = NonlinearSystem(eqs, [u,t], [])
jac = ModelingToolkit.calculate_jacobian(ns)

but the output is just the derivative of the equation

1×2 Array{Expression,2}:
 Constant(0)  derivative(u(10 + t), t) + derivative(u(t), t)

Am I misusing the API?

[seadra]

(Also, if I can figure out how I can use ModelingToolkit.jl, can I contribute a continuous difference equation solver for such equations, that uses IterativeSolvers.jl as its backend?)

[ChrisRackauckas]

You would create the array of u: @variables u[1:n](t) and write down the generating formula for the equations on the discretization. The implementation of runge_kutta_discretize might be informative:

https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/control/controlsystem.jl#L117-L168

[seadra]

Thanks for the pointers!

However, if I'm not misunderstanding, runge_kutta_discretize does the generating formula/discretization manually which is exactly the thing I'm trying to avoid doing using ModelingToolkit.jl. Does it mean I can't use ModelingToolkit to handle the discretization automatically (which is the only thing I actually need from ModelingToolkit)? (There are no derivatives in my equation so there is no stencil involved)

I'm also a bit confused about the @variables u[1:n](t), the index of u already denotes discretized time, so why would each variable still have t?

[ChrisRackauckas]

What are you trying to do then? You either need to write down a discretization somewhere, or handle it continuously via a neural network.

[seadra]

I was hoping to leverage ModelingToolkit to do discretization, which is the hard part of the problem. Since it is capable of discretizing PDEs, I thought it could also maybe do it for a continuous, non-differential "delay". Apparently not, so thanks for clarifying that.

As I said, I already solved my problem by doing the discretization by hand (which was some error-prone labor, but I did it for my own problem), and then used IterativeSolvers.jl to obtain a solution. I thought it would be nice to have a general solver for this kind of problems somewhere in SciML someday though.

[ChrisRackauckas]

The ModelingToolkit PDE discretization stuff is trying to get there, but indeed it doesn't automate this yet. My hope is that in a year this will be automatic though, where you specify the PDE and boom it spits out a discretization that you solve with IterativeSolvers.jl, NLsolve.jl, or DifferentialEquations.jl depending on whether the discretization built a linear problem, a nonlinear problem, or an ODE. We're building that as a multi-repo effort right now, but it's not all ready just yet.

[azev77]

FYI, solving (stochastic) difference equations is at the heart of macroeconomics. @RJDennis's excellent package, SolveDSGE, has many options for this.