# Solving difference equation: Part 2

**URL:** <https://discourse.julialang.org/t/solving-difference-equation-part-2/67057>\
**Category:** Performance\
**Tags:** difference-equations\
**Created:** [August 26, 2021, 8:36pm UTC](https://discourse.julialang.org/t/solving-difference-equation-part-2/67057 "2021-08-26T20:36:12Z")\
**Posts on this page:** 5\
**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 26, 2021, 8:36pm UTC](https://discourse.julialang.org/t/solving-difference-equation-part-2/67057/1 "2021-08-26T20:36:12Z")

</div>

In my [previous post](https://discourse.julialang.org/t/solving-difference-equation/66977) I was trying to get to reach a conclusion whether I can get more speed in solving a system of difference equation (basically a discrete dynamical system) through a hand-written code. I got useful suggestions, but today I found that the `trajectory()` function in `DynamicalSystems` runs faster. My ultimate goal was to write a code that solves the system in a biparametric space (varying two parameters in nested for loop) and store the corresponding period in a file. Here is the system and its jacobian

```julia
using DynamicalSystems, DelimitedFiles, .Threads, BenchmarkTools
## Components of a test DiscreteDynamicalSystem

function dds_constructor(u0 = [0.5, 0.7]; r=1.0, k=2.0)
    return DiscreteDynamicalSystem(dds_rule, u0, [r, k], dds_jac)
end
## equations of motion:
function dds_rule(x, par, n)
   r, k = par;
   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
## Jacobian:
function dds_jac( x, par, n)
    r, k = par;
    a = 5.0; mu = 0.5; d = 0.2;
    J11 = exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)) - x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1))*(r/k - (2*x[1]*x[2])/(x[1]^2 + a)^2)
    J12 = -(x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)))/(x[1]^2 + a)
    J21 = x[2]*exp((mu*x[1])/(x[1]^2 + a) - d)*(mu/(x[1]^2 + a) - (2*mu*x[1]^2)/(x[1]^2 + a)^2)
    J22 = exp((mu*x[1])/(x[1]^2 + a) - d)
    return @SMatrix [J11 J12; J21 J22]
end

```

Some functions I used:

```julia
function seqper_new(x; tol=1e-3) # function to calculate periodicity 
    @inbounds for k in 2:length(x)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:lastindex(x)) && return k - 1
        end
    end 
    return length(x)
end

```

```julia
function meshh(x,y) # create a meshgrid of x and y
    len_x = length(x);
    len_y = length(y);
    xmesh = Float64[]
    ymesh = Float64[]
    for i in 1:len_x
        for j in 1:len_y
            push!(xmesh, x[i])
            push!(ymesh, y[j])
        end
    end
    return xmesh, ymesh 
end

```

Here is my function that I want to optimize. It calculates the period of each pair of xpts\*ypts parameter values and store the corresponding x-values in a column of some preallocated variable.

```julia
function isoperiodic_test(ds, par_area::Matrix{Float64}, NIter::Int64, filename::String,
                         nxblock::Int64, nyblock::Int64, NTr::Int64, xpts::Int64, ypts::Int64)    

    io = open(iso_filename, "a") # opeing file to write output
    p1_st = par_area[1]; # beginning of first parameter 
    p1_nd = par_area[2]; # end of first parameter
    p2_st = par_area[3]; # beginning of second parameter
    p2_nd = par_area[4]; # end of second parameter

    p1_block = range(p1_st, p1_nd, length=nxblock + 1); # divide par_area in blocks 
    p2_block = range(p2_st, p2_nd, length=nyblock + 1); # divide par_area in blocks 
    l_bipar = xpts * ypts; # total number of pts in each block
    total_blocks = nxblock * nyblock # total number of blocks
    sol_last = Array{Float64,2}(undef, NIter + 1, l_bipar); # stores x values of solution for l_bipar (r,k) paris

    for ii in 1:nxblock
         step1size = (p1_block[ii + 1] - p1_block[ii]) / xpts;
         par1range = range(p1_block[ii], p1_block[ii + 1] - step1size, length=xpts);
        for jj in 1:nyblock
            display("Progress: " * string(nxblock * (ii - 1) + jj) * " out of " * string(total_blocks) * " steps.")
            step2size = (p2_block[jj + 1] - p2_block[jj]) / ypts;
            par2range = range(p2_block[jj], p2_block[jj + 1] - step2size, length=ypts);   
            Threads.@threads for i = 1:xpts         
                for j = 1:ypts
                     set_parameter!(ds, [par1range[i], par2range[j]]) # change the parameter values
                     tr = trajectory(ds, NIter; Ttr=NTr); # find the solution
                     sol_last[:,(i - 1) * xpts + j] = tr[:,1]; # store x-values
                end
            end
            periods = seqper_new.(eachcol(sol_last), tol=0.001) # seqper_new calculates periodicity 
            par1mesh, par2mesh = meshh(par1range, par2range) # meshgrid of parameter values
            writedlm(io, [par1mesh par2mesh periods]) # write data in txt file
        end 
    end
    close(io)
    return iso_filename
end

```

```julia
## Test
 
parameter_area = [1.0 5.0 2.0 5.0]
nxblock = 1;
nyblock = 1;
xpts = 100;
ypts = 100;
NIter = 2000;
NTr = 50000;
init = [0.4, 0.5];
ds = dds_constructor(init);
##
iso_filename = @time isoperiodic_test(ds, parameter_area, NIter, "isoperiodic_test_data.txt",
nxblock, nyblock, NTr, xpts, ypts)

```

Observation:

```julia
"Progress: 1 out of 1 steps."
  4.281922 seconds (100.33 k allocations: 628.321 MiB, 17.13% gc time)
"isoperiodic_test_data.txt_1.txt"

```

for 2 × 2 blocks

```julia
"Progress: 1 out of 4 steps."
"Progress: 2 out of 4 steps."
"Progress: 3 out of 4 steps."
"Progress: 4 out of 4 steps."
 16.914337 seconds (401.27 k allocations: 2.007 GiB, 10.98% gc time)
"isoperiodic_test_data.txt_1.txt"

```

Am I on the right track? I am still learning. Trying my best to find answers by myself but often find it difficult optimizing my own codes. Any form of help is appreciated.

If you find this difficult to read/understand let me know.

---

<div class="post-metadata">

**Author:** ![genkuroki](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/genkuroki/32/18030_2.png) [@genkuroki](https://discourse.julialang.org/u/genkuroki)\
**Post date:** [August 28, 2021, 7:42am UTC](https://discourse.julialang.org/t/solving-difference-equation-part-2/67057/2 "2021-08-28T07:42:01Z")

</div>

There are a few bugs in your code.

- The code `for k in 2:length(x)` in the function `seqper_new(x; tol=1e-3)` must be corrected as `for k in 2:(length(x) ÷ 2 + 1)`.
- `nxblock * (ii - 1) + jj` must be `nyblock * (ii - 1) + jj`.
- `(i - 1) * xpts + j` must be `(i - 1) * ypts + j`
- The `Threads.@threads for` loop in the function `isoperiodic_test` is not thread-safe, so that the result changes with each execution. Threads are conflicting at the line `set_parameter!(ds, [par1range[i], par2range[j]])`, where the same `ds` will be modified by multiple threads.

Memory allocations are also unnecessarily too much.

- The size of the array `sol_last` in memory is (NIter + 1) \* xpts \* ypts \* 8 bytes = 2001 \* 100 \* 100 \* 8 bytes = 153 MB, which is quite large, but it is discarded at the end of the function `isoperiodic_test`.
- Each execution of the `DynamicalSystems.trajectory` function results in a wasted memory allocation, which is repeated tens of thousands of times.

Other points that need improvement:

- The calculation results are saved to a file, and the results cannot be reused without going through the file.
- The code is unnecessarily complicated.
- You have created the complicated function `isoperiodic_test` to do most of all the work you want to do.

Furthermore, it is not necessary to use `DynamicalSystems.jl` in this case.

**My test implementation:**

```julia
using Base.Threads
using Plots
pyplot() # https://github.com/JuliaPlots/Plots.jl/issues/3560
using CSV
using DataFrames

```

```julia
function f(x, p)
   r, k = p
   a, mu, d = 5.0, 0.5, 0.2
   xnew1 = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
   xnew2 = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
   xnew1, xnew2
end

function traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
    x = x0
    for _ in 1:nwarmup
        x = f(x, p)
    end
    @inbounds tr1[1], tr2[1] = x
    for i in 2:niters+1
        x = f(x, p)
        @inbounds tr1[i], tr2[i] = x
    end
end

```

```julia
function findperiod(x; tol=1e-3)
    n = length(x)
    @inbounds for k in 2:(n ÷ 2 + 1)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:n) && return k - 1
        end
    end 
    return n
end

```

The function that calculates the period corresponding to the parameters (It does not save the result to a file):

```julia
function calcperiod2d(f = f, niters = 2000, nwarmup = 50000, x0 = (0.4, 0.5),
        r = (1:0.05:5)[1:end-1], k = (2:0.05:5)[1:end-1], tol=1e-3
    )
    tmp1 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    tmp2 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    period = Matrix{Int}(undef, length(r), length(k))
    @threads for j in eachindex(k)
        tid = threadid()
        tr1, tr2 = tmp1[tid], tmp2[tid]
        for i in eachindex(r)
            p = (r[i], k[j])
            traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
            @inbounds period[i, j] = findperiod(tr1; tol)
        end
    end
    period, r, k
end

```

Plotting function (I want to see the results immediately by plotting.):

```julia
xm1ox(x) = (x - 1)/x # x minus 1 over x : {x ≥ 1} -> {0 ≤ y < 1}
inv1my(y) = 1/(1 - y) # inverse of 1 minus 1 : {0 ≤ y < 1} -> {x ≥ 1}

function plot_period(r, k, period; kwargs...)
    cbtick = [1, 2, 3, 4, 5, 6, 10, 30]
    colorbar_ticks = (xm1ox.(cbtick), string.(cbtick))
    heatmap(r, k, xm1ox.(period)'; xlabel="r", ylabel="k", title="period",
        colorbar_ticks, kwargs...)
end

```

Save and load functions (I rarely use them, but many people may want to keep the results in a file.):

```julia
function save_period(r, k, period; fn = "result.csv")
    ii, jj = reim(complex.(axes(r, 1), axes(k, 1)'))
    rr, kk = reim(complex.(r, k'))
    df = DataFrame(i = vec(ii), j = vec(jj), r = vec(rr), k = vec(kk), period = vec(period))
    CSV.write(fn, df)
end

function load_period(fn = "result.csv")
    df = CSV.read(fn, DataFrame)
    imax = maximum(df.i)
    jmax = maximum(df.j)
    r = reshape(df.r, imax, jmax)[1:imax]
    k = reshape(df.k, imax, jmax)[1:imax:end]
    period = reshape(df.period, imax, jmax)
    period, r, k
end

```

**Example:**

```julia
period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
plot_period(r, k, period)

```

```julia
  3.109461 seconds (517.41 k allocations: 33.035 MiB, 7.54% compilation time)
  2.870107 seconds (228 allocations: 433.594 KiB)
  2.909849 seconds (221 allocations: 433.531 KiB)

```

![result](https://global.discourse-cdn.com/julialang/original/3X/9/f/9f1755b9d9ed54035a3c31eafa88f5da6025b770.png)

Calculate, save, load, and plot the result:

```julia
period, r, k = calcperiod2d()
fn = save_period(r, k, period; fn = "result.csv")
period, r, k = load_period(fn)
plot_period(r, k, period)

```

![result](https://global.discourse-cdn.com/julialang/original/3X/9/f/9f1755b9d9ed54035a3c31eafa88f5da6025b770.png)

Jupyter notebook: [https://github.com/genkuroki/public/blob/main/0018/DynamicalSystem/On%20Solving%20difference%20equation%20Part%202.ipynb](https://github.com/genkuroki/public/blob/main/0018/DynamicalSystem/On%20Solving%20difference%20equation%20Part%202.ipynb)

I love mathematics and undoubtedly, so do the members of this community.

If you find something mathematically interesting, please let us know. 😊

---

<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 28, 2021, 8:22am UTC](https://discourse.julialang.org/t/solving-difference-equation-part-2/67057/3 "2021-08-28T08:22:36Z")

</div>

Thank you very much @genkuroki , you are great.

- The reason I was saving the data is to plot them in MATLAB (by calling a user-defined function from Julia) as I am not very good at plotting them in Julia.
- I used `DynamicalSystems.jl` so that I can modify the code a little bit to also calculate the Lyapunov spectrum for all those parameter values and compare the iso-periodic diagram with the Lyapunov one. One way to do it is using `lyapunovspectrum` from `DynamicalSystems.jl` library, otherwise I have to write the code for finding lypunov spectrum myself.

What should I do? Should I open an issue in `DynamicalSystems.jl` to add another function `trajectory!` that takes a preallocated trajectory data and outputs the modified one. Will that be useful?

What if I want to use `DynamicalSystems.jl` along with `@threads`, can I make it threadsafe by defining the dynamical system inside the for loop where I put the new parameter values?

---

<div class="post-metadata">

**Author:** ![HerAdri](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heradri/32/5816_2.png) [@HerAdri](https://discourse.julialang.org/u/HerAdri)\
**Post date:** [August 29, 2021, 8:07am UTC](https://discourse.julialang.org/t/solving-difference-equation-part-2/67057/4 "2021-08-29T08:07:31Z")

</div>

> [@genkuroki](#):
>
> Calculate, save, load, and plot the result:  
> How to resolve this alert?
> 
> ```julia
> period, r, k = calcperiod2d()
> fn = save_period(r, k, period; fn = "result.csv")
> period, r, k = load_period(fn)
> plot_period(r, k, period)
> 
> ```

```julia
julia> plot_period(r, k, period)
Error showing value of type Plots.Plot{Plots.PyPlotBackend}:
ERROR: MethodError: no method matching getproperty(::Char, ::String)
Closest candidates are:
  getproperty(::DataFrames.GroupKey, ::AbstractString) at C:\Users\Hermesr\.julia\packages\DataFrames\vuMM8\src\groupeddataframe\groupeddataframe.jl:551
  getproperty(::DataFrames.DataFrameColumns, ::AbstractString) at C:\Users\Hermesr\.julia\packages\DataFrames\vuMM8\src\abstractdataframe\iteration.jl:212
  getproperty(::DataFrameRow, ::AbstractString) at C:\Users\Hermesr\.julia\packages\DataFrames\vuMM8\src\dataframerow\dataframerow.jl:300
  ...
Stacktrace:
  [1] py_set_axis_colors(sp::Plots.Subplot{Plots.PyPlotBackend}, ax::PyCall.PyObject, a::Plots.Axis)
    @ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\backends\pyplot.jl:815
  [2] _before_layout_calcs(plt::Plots.Plot{Plots.PyPlotBackend})
    @ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\backends\pyplot.jl:1100
  [3] prepare_output(plt::Plots.Plot{Plots.PyPlotBackend})
    @ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\plot.jl:177
  [4] display(d::VSCodeServer.InlineDisplay, m::MIME{Symbol("image/svg+xml")}, x::Plots.Plot{Plots.PyPlotBackend})
    @ VSCodeServer c:\Users\Hermesr\.vscode\extensions\julialang.language-julia-1.3.32\scripts\packages\VSCodeServer\src\display.jl:0
  [5] display(d::VSCodeServer.InlineDisplay, mime::String, x::Any)
    @ Base.Multimedia .\multimedia.jl:216
  [6] display(d::VSCodeServer.InlineDisplay, x::Plots.Plot{Plots.PyPlotBackend})
    @ VSCodeServer c:\Users\Hermesr\.vscode\extensions\julialang.language-julia-1.3.32\scripts\packages\VSCodeServer\src\display.jl:109
  [7] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:328
  [8] #invokelatest#2
    @ .\essentials.jl:708 [inlined]
  [9] invokelatest
    @ .\essentials.jl:706 [inlined]
 [10] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:247
 [11] (::REPL.var"#40#41"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:231
 [12] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:462
 [13] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:229
 [14] (::REPL.var"#do_respond#61"{Bool, Bool, REPL.var"#72#82"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:798
 [15] #invokelatest#2
    @ .\essentials.jl:708 [inlined]
 [16] invokelatest
    @ .\essentials.jl:706 [inlined]
 [17] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\LineEdit.jl:2441
 [18] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\REPL\src\REPL.jl:1126
 [19] (::REPL.var"#44#49"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL .\task.jl:411

```

installed version of matplotlib

```julia
Z:\>pip show matplotlib
Name: matplotlib
Version: 3.4.1
Summary: Python plotting package
Home-page: https://matplotlib.org
Author: John D. Hunter, Michael Droettboom
Author-email: matplotlib-users@python.org
License: PSF

```

---

<div class="post-metadata">

**Author:** ![genkuroki](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/genkuroki/32/18030_2.png) [@genkuroki](https://discourse.julialang.org/u/genkuroki)\
**Post date:** [August 29, 2021, 8:43am UTC](https://discourse.julialang.org/t/solving-difference-equation-part-2/67057/5 "2021-08-29T08:43:02Z")

</div>

Same as [Use Plots() library to graph a colored surface - #8 by genkuroki](https://discourse.julialang.org/t/use-plots-library-to-graph-a-colored-surface/66561/8)