# Solving difference equation

**URL:** <https://discourse.julialang.org/t/solving-difference-equation/66977>\
**Category:** Performance\
**Tags:** question, difference-equations\
**Created:** [August 25, 2021, 12:44pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977 "2021-08-25T12:44:26Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![MathGuy](https://avatars.discourse-cdn.com/v4/letter/m/58f4c7/32.png) [@MathGuy](https://discourse.julialang.org/u/MathGuy)\
**Post date:** [August 25, 2021, 12:44pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/1 "2021-08-25T12:44:26Z")

</div>

This is my code to solve a system of difference equations:

```julia
function trajectoryDiscrete(sys_eq, init_cond, NIter, par1, par2)
    neq = length(init_cond)
    traj= zeros(Float64,NIter+1,neq)
    traj[1,:] = init_cond
    sys(x) = sys_eq(x, par1, par2);
     for i in 1:NIter
        init_cond = sys(init_cond)
        traj[i+1,:] = init_cond
    end
    return traj
end

```

Here is a sample system of difference equations:

```julia
function test_difference_eq(x,r,k)
    a = 5.0; mu = 0.5; d = 0.2; 
    dx = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
    dy = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
    return @SVector [dx, dy]
end

```

Here is the result:

```julia
tr = @time trajectoryDiscrete(test_difference_eq, [0.5, 0.3], 10000, 1, 2)

```

```julia
  0.000414 seconds (3 allocations: 156.484 KiB)
10001×2 Matrix{Float64}:
  0.5 0.3
 13.2216 0.380422
  0.000982828 2586.64
  ⋮
  0.0 2.5e-323
  0.0 2.5e-323

```

Query: If I remove `@SVector` from the return statement of `test_difference_eq` then it gives

```julia
0.000655 seconds (10.00 k allocations: 1.068 MiB)
10001×2 Matrix{Float64}:
 0.5 0.3
 0.99971 0.257598
 1.5792 0.229228
 ⋮
 2.0 2.5e-323
 2.0 2.5e-323

```

Ah! you guessed it right, I am new to Julia, trying to settle here ditching MATLAB. My queries are:

1. What magic does `@SVector` do so that allocation is reduced. Can you give me an idea when, where and how I can use @SVector and @SMatrix?

2. Can this code be optimized further?

My ultimate goal is to solve the system for a range of `r=range(1, 5, length = 501)` and `k=range(2, 4, length = 501)` in nested for loop, store last few point of `tr` to find its period using another user-defined function `seqper` ([see my previous post](https://discourse.julialang.org/t/julia-implementation-of-seqperiod-from-matlab/66791/2)) and finally store each `r-k` pair along with the associated period in a text file.

What are your suggestions?

Any form of help is appreciated.

Blockquote

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [August 25, 2021, 12:50pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/2 "2021-08-25T12:50:25Z")

</div>

> [@MathGuy](#):
>
> What magic does `@SVector` do so that allocation is reduced

`SVector` is allocated on stack when it’s small, think of it as a `NTuple` basically, it’s especially useful when you need to use a small matrix / vector over and over again.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 25, 2021, 1:17pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/3 "2021-08-25T13:17:36Z")

</div>

> [@MathGuy](#):
>
> Can this code be optimized further?

Your code is row major. You may want to take a look at the 18.337 course notes which include optimizing a difference equation as one of its early lectures.

> **[GitHub - mitmath/18337: 18.337 - Parallel Computing and Scientific Machine...](https://github.com/mitmath/18337)**
>
> 18.337 - Parallel Computing and Scientific Machine Learning - GitHub - mitmath/18337: 18.337 - Parallel Computing and Scientific Machine Learning

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 25, 2021, 1:32pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/4 "2021-08-25T13:32:43Z")

</div>

@MathGuy SolveDSGE.jl solves infinite horizon systems of stochastic & deterministic difference equations

@ChrisRackauckas will SciML ever be able to solve difference equations? Or is it only intended for continuous time?

---

<div class="post-metadata">

**Author:** ![MathGuy](https://avatars.discourse-cdn.com/v4/letter/m/58f4c7/32.png) [@MathGuy](https://discourse.julialang.org/u/MathGuy)\
**Post date:** [August 25, 2021, 1:36pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/5 "2021-08-25T13:36:24Z")

</div>

See `DynamicalSystems.jl`. It solves dicrete problmes like above. Just use `trajectory`. For further details see this [link](https://juliadynamics.github.io/DynamicalSystems.jl/latest/ds/general/).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 25, 2021, 7:17pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/6 "2021-08-25T19:17:24Z")

</div>

> [@Albert\_Zevelev](#):
>
> @ChrisRackauckas will SciML ever be able to solve difference equations? Or is it only intended for continuous time?

> [@MathGuy](#):
>
> See `DynamicalSystems.jl` . It solves dicrete problmes like above. Just use `trajectory` . For further details see this [link](https://juliadynamics.github.io/DynamicalSystems.jl/latest/ds/general/).

DynamicalSystems.jl uses the DiffEq `DiscreteProblem` interface for that.

[https://diffeq.sciml.ai/stable/types/discrete\_types/](https://diffeq.sciml.ai/stable/types/discrete_types/)

It’s not the most talked about one because, well, it’s not so hard to write the loops, but it exists and puts everything into the same interface as the DiffEq solutions. It’s mostly used as the basis for the discrete jump problems, but you can also just directly define and solve a discrete dynamical system with it. It’s also what ModelingToolkit.jl’s `DiscreteSystem` targets.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [August 25, 2021, 7:32pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/7 "2021-08-25T19:32:30Z")

</div>

@ChrisRackauckas the simplest system economists solve is an infinite horizon, 2nd order, boundary value (not initial value), non-linear system of difference equations:  
\left[\begin{array}{c} (c\_{t})^{-\sigma } = \beta (\alpha (k\_{t+1})^{\alpha-1} + (1-\delta)) (c\_{t+1})^{-\sigma } \\ k\_{t+1} = k\_t^\alpha + (1-\delta)k\_{t} - c\_t \end{array}\right]  
Boundary values:  
Initial condition: k(0)=k\_0 given  
Terminal condition: \lim\_{t\to \infty} \beta^t (c\_{t})^{-\sigma } k\_{t+1} =0

This is easy to solve w/ `SolveDSGE.jl`. (include the stochastic extensions…)  
Can SciML handle this yet?  
Does SciML want to eventually?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 25, 2021, 7:42pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/8 "2021-08-25T19:42:09Z")

</div>

Yeah, I owe @jlperla some work getting state space models onto the SciML interface.

Is that implicit discrete?

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [August 25, 2021, 7:58pm UTC](https://discourse.julialang.org/t/solving-difference-equation/66977/9 "2021-08-25T19:58:26Z")

</div>

> [@Albert\_Zevelev](#):
>
> Can SciML handle this yet?  
> Does SciML want to eventually?

No, but it would be nice at some point. The mixed initial value/boundary value problems here are a little non-standard for people outside of economics though, so I think having solvers that layer on top of the sciml stuff is probably better for now. The steady state transversality conditions are also very different from the standard sciml cases, but will hopefully get there.

My guess is that the first step is to get the simulation of state space models (and their differentiation) correct first, which is something that we will hopefully make progress on this fall.