Mathematical modelling and

computer simulations

Projekt
Mathematical modelling and computer simulations winter semester: 2024/2025

Practical informations

The simulation (R code) and descriptive (pdf file) are evaluated. The project (both parts)
should be attached to the course of the subject with the appropriate activity in the Project
component.
Detailed information on the implementation of both parts of the project is provided below.

1. Simulation part.

– Each code attached to eLearning has to be signed.

– Please provide clear breakdowns into individual parts of the task in the code.
– You should check the compilation of each code fragment, pay attention to the variable
names (so that they do not overlap).
– The name of the simulation file: NameSurname Project1.R.

2. Theoretical part.

– The report should include: names and surnames, specialty, date of completion.
– We describe the project report in a specialized TeX or LaTeX mathematical environ-
ment, and during the course we attach a PDF file named NameSurname [Link].
– Please make sure that the required information is provided in each subsection of the
report. In addition, please provide conclusions and/or descriptions of the results or
graphs obtained. Please do not forget to sign each graph (especially for which data
it was made).
– Please note that the report’s execution method is also evaluated, so do not submit
“dry ”results or graphs.
– Please attach the corresponding program code to each point and all the code with
the description at the end of the report.

1
Mathematical modelling and computer simulations winter semester: 2024/2025

Pseudo-random numbers generation

A random number is a number X belonging to the set of Σ = {X1 , ..., Xn } values selected
with a certain probability. If X can be any number from the set with the same probability
p(X) = 1/n, then we are talking about the even probability distribution of random numbers
from the set Σ.
Simulating purely random number generation is an extremely difficult task. This is due to
the deterministic nature of the computer and the operations it performs. When a person rolls a
die, he does not know and cannot predict what will happen. The same operation on a computer
requires an action, the result of which is unpredictable - none of the operations performed by
the processor has this feature.
On the other hand, random numbers are commonly used in programming: games of chance,
simulations of various random or statistical processes, testing algorithms for random data sets,
encrypting images, etc. Since we cannot easily obtain real random numbers, we have to work
on their artificial counterpart - pseudorandom numbers. Pseudorandom numbers look like ran-
dom, but are not, because we use algorithms (a predetermined finite set of steps) to generate
them. This means that knowing the generational formula and a few consecutive pseudo-random
numbers, we can easily generate the rest of the necessary numbers - random numbers do not
have this feature, otherwise the lottery would lose its sense.
Therefore, computer pseudorandom number generators are in fact deterministic functions
that determine successive elements of the sequence of pseudorandom numbers as successive
iterations of the established algorithm. To start iteration of pseudorandom numbers, you need
a start value from which to start generating. We call this value pseudorandom seed.

PART I

I. Generating pseudo-random numbers with a given distribution. Many program-

ming languages, programming environments or statistical packages have built-in functions that
allow you to generate pseudo-random numbers. At the core of stochastic simulations is the
issue of generating pseudo-random numbers that imitate the behavior of uniformly distributed
random variables.

Therefore, we assume that we have at our disposal a potentially infinite sequence of inde-
pendent pseudo-random variables X1 , . . . , Xn , . . . with a uniform distribution of U (0, 1). There
is a function in the R package

runif(n) (random unif orm distribution).

Zadanie 1. Using the inverse distribution method, generate a sequence of pseudorandom num-
bers with the Extreme values distribution at a given density
(x−2)
f (x) = ex−2 e−e , x ∈ R.

1.1 Determine the cumulative distribution function and the inverse distribution function in
an analytical way.

2
Mathematical modelling and computer simulations winter semester: 2024/2025

1.2 Present both cumulative distribution functions in one diagram.

1.3 Generate at least three n-element samples, e.g. n = 100, 1000, 10000.

−⋆−

After determining the sequence of pseudo-random numbers with a given distribution, we

can proceed to its use / analysis. In our case, we will assume that these are experimental data.
We will make a statistical analysis of them.
Zadanie 2. For all generated samples:
2.1 Determine the selected position measures and compare with the parameters of the the-
oretical distribution.

2.2 Plot in one plot: the histogram (frequency plot) and the theoretical distribution density
plot (probability mass plot).

2.3 Plot in one graph: empirical and Extreme values distribution.

What can we infer from the above results?

−⋆−

II. Chaotic maps with a homogeneous distribution. Having a given sequence of pseudo
random numbers X = {x0 , x1 , . . . , xn } with any distribution (even unknown), we can generate
values from a distribution close to homogeneous using chaotic maps with homogeneous de-
composition.

One of the chaotic maps with a homogeneous distribution is skew tent map (skew tent
map), which is defined by a recursive relationship:
 zk

 , zk ∈ [0, p),

 p
zk+1 = f (zk ) = (1)
1 − zk
, zk ∈ [p, 1],



1−p


where p ∈ (0, 1), k = 0, 1, . . . , N , N is the number of iterations performed with the chaotic
map. Note that
z2 = f (z1 ) = f (f (z0 )) = f 2 (z0 ).
Thus, f N (x) stands for the N iteration of the chaotic map for the initial value of z0 = x. Let
f = {x
X en } normalize a given pseudo-random number sequence X.
e0 , . . . , x
Then the sequence U = {u0 , u1 , . . . , un }, the individual elements of which are determined
using the N -th iteration of the chaotic map with the appropriate initial value:

ui = f N (xei ), i = 0, 1, . . . , n, (2)

has a distribution close to homogeneous for the appropriate number of iterations N .

3
Mathematical modelling and computer simulations winter semester: 2024/2025

Zadanie 3. For the generated samples of sequences of pseudorandom numbers with the Extreme
values distribution, perform:

3.1 Normalize the samples.

3.2 For each normalized sample, generate a sequence of pseudorandom U using (2), where f N
denotes the N -th iteration of the chaotic map (1) for p = 0.45.

3.3 Generate U distribution histograms for each of the samples, for three N values (e.g. for
N = 1, 5, 12, it may be necessary to adjust the values experimentally).

3.4 In each case generate graphs of the dependence of pairs of points (uk , uk+1 ).

3.5 For what number of iterations of the chaotic map we get ”flattening”(approximation to a
homogeneous distribution)?

4
Mathematical modelling and computer simulations winter semester: 2024/2025

PART II

Chaotic maps are subtypes of nonlinear dynamic systems of a deterministic nature, the
solutions of which have behavior similar to a noise signal. The high sensitivity of the chaotic
map to the initial value is one of its most important features - butterfly effect. The butterfly
effect means that making a small change in the initial value results in a dramatic change in
the result. This property gives high effectiveness and unpredictability of the chaotic map in
applications, among others, from the point of view of security in the Internet network. Chaotic
maps are widely used to encode information sent over the network, e.g. encoding images, which
you will deal with in the second part of the project.
Image coding using chaotic maps can be done in many different ways. This issue is still
developing intensively and we find more and more proposed coding algorithms in the literature,
which achieve better and better results in tests.
In the project, you will use a chaotic map to encode the image parrot [Link] (available
on the course website in the project materials) by shuffling the color intensity positions. Each
color image is N · M · 3 pixels, i.e. height·width·channels, the image is stored as an array
of dimension N × M ×3. This means that each pixel stores its own RGB intensities. The task
description provides instructions for completing the project in R.
The procedure algorithm can be found in the article [1].

Image encoding and decoding:

1. We save and load the image (you can use the readImage() function from the OpenImageR
package, e.g. img = readImage("parrots [Link]")).

2. We extract three matrices {R, G, B} (e.g. using readImage() we get the red intensity
matrix for each pixel using img[,,1]). We write each matrix {R, G, B} as a vector of
length N · M , which gives {Rw, Gw, Bw}.

3. For each color in the original image, we can display the intensity (in bar form, example
for all colors) and present histograms (example for red).

4. Using a logistic map

xn+1 = rxn (1 − xn ), n = 0, . . . , K,

5
Mathematical modelling and computer simulations winter semester: 2024/2025
300

300

300
250

250

250
200

200

200
Gw * 255
Rw * 255

Bw * 255
150

150

150
100

100

100
50

50

50
0

0
0e+00 2e+04 4e+04 6e+04 8e+04 1e+05 0e+00 2e+04 4e+04 6e+04 8e+04 1e+05 0e+00 2e+04 4e+04 6e+04 8e+04 1e+05

Index Index Index

Histogram of Rw
2000
Frequency

1000
500
0

0.2 0.4 0.6 0.8 1.0

Rw

for the selected r ∈ [3.6, 4] and x0 ∈ (0, 1) we construct a vector of values of length
K = 3 · N · M . We divide the obtained vector into three subvectors {xR , xG , xB }. The
encoding key consists of two pieces of information

(x0 , r).

5. We sort (descending or ascending) the vectors {xR , xG , xB } separately, obtaining {xsR , xsG , xsB }.
We then create three position vectors {PR , PG , PB } that store the item position numbers
of {xsR , xsG , xsB } in the vectors {xR , xG , xB }.

7. (Coding) We shuffle the values of vectors {Rw, Gw, Bw} using position vectors {PR , PG , PB },
i.e. we create new vectors {TR , TG , TB }, into which we insert the values from {Rw, Gw, Bw}
as follows: we take the element from the i−th position of the vector Rw, i.e. Rw[i], we
look for what position i is in the PR vector (let’s assume it’s j), then we put Rw[i] in
position j of the TR vector, i.e. TR [j].

6
Mathematical modelling and computer simulations winter semester: 2024/2025

8. We combine the color vectors {TR , TG , TB } into an array (function array()) with dimen-
sion N × M × 3. The encoded image looks like this:

9. We plot the intensity of each color of the encoded image

300

300

300
250

250

250
200

200

200
TG * 255
TR * 255

TB * 255
150

150

150
100

100

100
50

50

50
0

0e+00 2e+04 4e+04 6e+04 8e+04 1e+05 0e+00 2e+04 4e+04 6e+04 8e+04 1e+05 0e+00 2e+04 4e+04 6e+04 8e+04 1e+05

Index Index Index

10. (Decode) The receiver, having the encoding key (x0 , r), can create position vectors {PR , PG , PB },
which will allow to arrange {DR , DG , DB }innewvectors(decodedcolors). values from {TR , TG , TB }.
Then we combine {DR , DG , DB } into an array to obtain the decoded image:

7
Mathematical modelling and computer simulations winter semester: 2024/2025

11. We determine the correlation coefficient between the individual colors of the original
image and the colors of the encoded image, e.g. for the color red we have
N
X ·M   
Rwi − Rw (TR )i − TR
i=1
corR = v
uN ·M 
v
uX 2 u ·M 
uNX 2
t Rwi − Rw t (TR )i − TR
i=1 i=1

12. We check the sensitivity of the encoding in relation to the key selection, for example, we
encode the image with the key (x0 = 0.1, r = 3.8), then we try to decode the image using
the key (x0 = 0.1000000001, r = 3.8) and (x0 = 0.1, r = 3.800000001).

Literatura
[1] M. Prasad, K.L. Sudha, Chaos Image Encryption using Pixel shuffling (artykuł umiesz-
czony na stronie kursu w materiałach do projektu).