# Solving difference equation: Part 2

In my previous post 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

``````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*exp(r*(1-x/k)-x/(a+x^2))
dy = x*exp(mu*x/(a+x^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/(x^2 + a) - r*(x/k - 1)) - x*exp(- x/(x^2 + a) - r*(x/k - 1))*(r/k - (2*x*x)/(x^2 + a)^2)
J12 = -(x*exp(- x/(x^2 + a) - r*(x/k - 1)))/(x^2 + a)
J21 = x*exp((mu*x)/(x^2 + a) - d)*(mu/(x^2 + a) - (2*mu*x^2)/(x^2 + a)^2)
J22 = exp((mu*x)/(x^2 + a) - d)
return @SMatrix [J11  J12; J21 J22]
end
``````

Some functions I used:

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

``````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;          # beginning of first parameter
p1_nd =  par_area;          # end of first parameter
p2_st =  par_area;          # beginning of second parameter
p2_nd =  par_area;          # 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);
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
``````
``````##  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:

``````"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

``````"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.

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:

``````using Base.Threads
using Plots
pyplot() # https://github.com/JuliaPlots/Plots.jl/issues/3560
using CSV
using DataFrames
``````
``````function f(x, p)
r, k = p
a, mu, d = 5.0, 0.5, 0.2
xnew1 = x*exp(r*(1-x/k)-x/(a+x^2))
xnew2 = x*exp(mu*x/(a+x^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, tr2 = x
for i in 2:niters+1
x = f(x, p)
@inbounds tr1[i], tr2[i] = x
end
end
``````
``````function findperiod(x; tol=1e-3)
n = length(x)
@inbounds for k in 2:(n ÷ 2 + 1)
if abs(x[k] - x) ≤ 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):

``````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))
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.):

``````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.):

``````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

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:

``````period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
plot_period(r, k, period)
``````
``````  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)
`````` Calculate, save, load, and plot the result:

``````period, r, k = calcperiod2d()
fn = save_period(r, k, period; fn = "result.csv")
plot_period(r, k, period)
`````` I love mathematics and undoubtedly, so do the members of this community.

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

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?

``````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:
 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
 _before_layout_calcs(plt::Plots.Plot{Plots.PyPlotBackend})
@ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\backends\pyplot.jl:1100
 prepare_output(plt::Plots.Plot{Plots.PyPlotBackend})
@ Plots C:\Users\Hermesr\.julia\packages\Plots\hyS17\src\plot.jl:177
 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
 display(d::VSCodeServer.InlineDisplay, mime::String, x::Any)
@ Base.Multimedia .\multimedia.jl:216
 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
 display(x::Any)
@ Base.Multimedia .\multimedia.jl:328
 #invokelatest#2
@ .\essentials.jl:708 [inlined]
 invokelatest
@ .\essentials.jl:706 [inlined]
 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
 (::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
 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
 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
 (::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
 #invokelatest#2
@ .\essentials.jl:708 [inlined]
 invokelatest
@ .\essentials.jl:706 [inlined]
 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
 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
 (::REPL.var"#44#49"{REPL.LineEditREPL, REPL.REPLBackendRef})()
``````Z:\>pip show matplotlib