Delaunay Triangulation covering of a Rossler Attractor (using MATLAB)
Published:
This post discusses performing Delaunay triangulation on a Rossler Attractor system. Details are given in ee520-hw8.pdf.
```{r, echo=FALSE} <!DOCTYPE html>
Homework Assignment № 8
Kevin Mack, Clarkson University 12/14/2017
Question 1:
In the first task of the assignment, there is Matlab code provided from section A3 of the appendix
from Applied and Computational Measurable Dynamics by Erik M. Bollt and Naratip Santitis-
sadeekorn. The code is used to perform several tasks on data sets from chaotic systems.
The tasks include Delaunay Triangulation, and use of the Ulam-Galerkin matrix to estimate the
Frobenius-Perron operator. The methods will be used on a Lorenz attractor given by (1), the
standard map in (2), and finally the Rossler attractor in (5).
The Lorenz attractor is a system of ordinary differential equations which produce chaotic
solutions for certain parameter values of parameters (σ, ρ, and β) and initial conditions,
dx
dt
= σ(y − x),
dy
dt
= x(ρ − z) − y,
dz
dt
= xy − βz.
(1)
The code to produce the Delaunay triangulation of the Lorenz attractor is given in appendix
A.3 (runnerLorenzCover.m) and is shown in the first code listing.
Listing 1: Matlab code – runnerLorenzCover.m
1 c l e a r ; c l o s e a l l ;
2 load ’ LorenzDat . ma t ’ %This shou l d be n x 3 array o f r e a l numbers in t h i s ←-
ca s e being p oi nt s on ( near ) the l o r e n z a t t r a c t o r
3 f i g u r e
4 hold on
5 g r i d on
6 g r i d minor
7 p lo t 3 (X( : , 1 ) ,X( : , 2 ) ,X( : , 3 ) , ’ . ’ )
8 z=X;
9
10 xlow=f l o o r ( min (X) ) ; xhigh=c e i l (max(X) ) ; h = 2 . 5 ;
11 [ X1, X2, X3] = ndgrid ( xlow (1 ) : h : xhigh ( 1) , xlow ( 2 ) : h : xhigh ( 2 ) , xlow ( 3) : h : ←-
xhigh ( 3) ) ;
12 m=s i z e (X1) ;
13
14 x1=reshape (X1 ,m( 1 ) *m(2 ) *m( 3) , 1) ; x2=r e s h a p e (X2 ,m( 1) *m( 2) *m( 3) ,1 ) ; x3=←-
reshape (X3,m( 1 ) *m( 2) *m( 3 ) ,1 ) ;
Assignment № 8 Page 1
15
16
17 %Formulate Delaunay T r ia n gu la t io n o f r e g i on
18 dt= DelaunayTri ( [ x1 x2 , x3 ] ) ; %See Matlab Subroutine DelauneyTri f o r ←-
input / output i n fo r ma t ion
19 %dt i s the t r i a n g u l a t i o n c l a s s
20 %
21 % where
22 %
23 %dt . dt . T ri a ng u la t io n i s an m1 by 4 array o f
24 %i n t e g e r s l a b e l l i n g the v e rte x c orn e r numbers
25 %o f the t r i a n g l e s
26 %
27 % and
28 %
29 %dt .X i s an m2 by 3 a rray o f r e a l numbers
30 %d e f i n i n g p o s i t i o n s
31
32
33
34 %Count number of o r b i t p oi n t s i n z which cause a t r i a n g l e to be counted as
35 %occupied ( and o th e rw i se a t r i a n g l e i s not counted as i t i s empty u n t i l
36 %obser v e d o c c u p i e d
37
38 n o ttru e =0;
39 w h i le ( n ottrue <1)
40 no t true =1;
41 SI = p o i ntL o cat i on ( dt , z ) ; %Matlab command , Locate the simplex element ←-
in dt
42 %c on ta i ni ng the s p e c i f i e d l o c a t i o n s o f each
43 %o f the e l e m e n t s o f the a r r a y z o f o r b i t
44 %samples
45 %
46 l=unique ( SI ) ; %Matlab s ubro u t ing : Count the number o f unique ←-
i n s t a n c e s i n SI
47 k=1;
48 wh i l e ( i snan ( l (k ) )<1&&k<l en g t h ( l ) ) %Using Matlab sub r ou t ine True f o r ←-
Not - a - Numb er
49 [ i i , j ]= f i n d ( SI==l ( k ) ) ; %C o l l e c t t h o se l o c a t i o n s c o rre s pon d ing to ←-
each unique l .
50 cnt ( k )=sum ( j ) ;
51 k=k+1;
52 end
53 k=k - 1 ; l l =l ( 1 : k ) ; cnt=reshape ( cnt , s i z e ( l l ) ) ;
54 end
55 %%
56
Assignment № 8 Page 2
57 % P l o t t hose sim p l e x e l e m e n t s o f the dt which a r e occupied by an o r b i t
58 % i t e r a t e o f z - dt ( l l , : ) a r e t h o se occupied simplex elements
59 f i g u r e ;
60 h old on
61 g r i d on
62 g r i d minor
63 p lo t 3 (0 ,0 , 0 ) ; patch ( ’ f a c e s ’ , dt ( l l , : ) , ’ v e r t i c e s ’ , dt .X, ’ FaceColor ’ , ’ r ’ ) ;
64 N=l eng t h ( l l ) ;
65
66 %Creating and Edit i n g Delaunay T ri a n g u l a t io ns
The code called runnerLorenzCover.m produces Figure 1, which shows several angles of
the Lorenz attractor as it is covered by a Delaunay triangulation algorithm. The algorithm can
be adjusted based on how much resolution is needed for the cover (by changing the tessellation
size, h).
Further, the Ulam-Galerkin matrix is used to estimate the Frobenius-Perron operator and will
be used to evaluate the standard map. The standard map is an area-preserving chaotic map,
where p
n
and θ
n
are taken modulo 2π,
p
n+1
= p
n
+ Ksin(θ
n
),
θ
n+1
= θ
n
+ p
n+1
.
(2)
The code to produce the standard map is given in the listing below, and can be found in the
appendix A.1.2.
Listing 2: Matlab code – standard.m
1 f u n c t i o n xvecNew = s tandar d ( xvec , k )
2 % Matlab coded i n v e c t o r i z e d form o f stand ard
3 %Input :
4 % xvec - v e ct o r o f i n i t i a l c o n d i t i o n s
5 % k - stand ard map parameter
6 %
7 %Output :
8 % xvecNew - ve ct o r o f t h e i r images
9 %
10 xvecNew = z er os ( s i z e ( xvec ) ) ;
11 xvecNew ( 1 , : ) = xvec ( 1 , : )+xvec ( 2 , : ) -k* s i n (2* p i * xvec ( 1 , : ) ) /(2* p i ) ;
12 xvecNew ( 2 , : ) = xvec ( 2 , : ) -k* s i n ( 2* p i * xvec ( 1 , : ) ) /(2* p i ) ;
The code to create the stochastic matrix, A, of the standard map is called TransitionMa-
trix.m and is tested by runnerSimpleTriCover.m (appendix A.1.2), which are both given in the
listings below. The runnerSimpleTriCover.m code makes a sample orbit of the standard map,
then calls TransitionMatrix.m to use the Ulam-Galerkin method to approximate the Frobenius-
Perron operator in phase space after being thresholded with second eigenvector, v
2
. The output
Assignment № 8 Page 3
(a) Angle 1, Lorenz Attractor (b) Angle 1, Lorenz Attractor outer cover h = 2.5
(c) Angle 2, Lorenz Attractor (d) Angle 2, Lorenz Attractor outer cover h = 2.5
Figure 1: The images in the left column show plotted points of the Lorenz attractor for a specific
set of parameters and initial conditions. The images in the right column represent the outer
cover of the attractor which has been achieved through Delaunay triangulation.
from the code can be seen for different values of k in Figures 2 and 3 in the next section.
Listing 3: Matlab code – runnerSimpleTriCover.m
1 c l o s e a l l ; c l e a r ;
2 nPlot = 1 0000; % nPlot - i t e r a t e t e s t o r b i t t h i s many po i n t s a f t e r the ←-
t r a n s i e n t .
3 M=250; diam=10; a =0; b=1;
4 h=0.025
5
6 %%
7 %%%%%%%Make a Sample Orbit o f the Standard Map
8 k =1 . 2 ; %k i s chosen > 0 . 9 7 . . . f o r the magic breakup o f the " golde n mean" ←-
t o ru s
9 z = [ ] ;
Assignment № 8 Page 4
10 % Set some v a r i a b l e s
11 t r a n s i e n t S t e p s = 10 0 ; %l e ngt h o f i n i t i a l t r a n s i e n t o f t e s t o r b i t to ←-
i g n o r e
12 i n i t i a l X = rand ( 2 , 1) ;
13 x = i n i t i a l X ; hold on ;
14
15 %%
16 % Throw away t r a n s i e n t s :
17 f o r i =1: t r a n s i e n t S t e p s
18 x = stan dard ( x , k ) ; %Use as a t e s t example , the stand ard map .
19 end
20
21 %%
22 % P l o t the next nPlot p o i n t s v i s i t e d :
23 f o r i =1: nPlot
24 x = stan dard ( x , k ) ;
25 z=[ z ; x ‘ ] ;
26 end
27 z=mod( z , 1 ) ;
28 p l o t ( z ( : , 1 ) , z ( : , 2 ) , ‘ b . ‘ , ‘ markersize ‘ , 1 0 ) ; hold on ; %See Fig . 1 1 . 1 .
29 %
30 %%
31
32
33 %%%%%%%%%%%%%%%
34 %%%Build an Ula m - Galer kin ‘ s Matrix Based on a Test Orb it V i s i t i n g ←-
T r i a n g l e s
35 %%%of a T e s s e l at i o n
36 [ dt , l l ,A, zz ]= T r ansi t i onMat r i x ( z , h , a , b , a , b) ;
37 t r i p l o t ( dt ) %Matlab r o u ti ne t o draw the t r i a n g u l a t i o n simplex
38 %See Fig . 1 1 . 1 .
39 drawnow ;
40
41 %%
42 [ v , d]= e i g s (A‘ , 2 ) ; w=abs (v ( : , 1 ) ) ; %Compute 1 s t and send e i ge n v a l u e s / ve c t o rs ←-
o f Galerkin - Ulam Matrix A
43 %
44 f i g u r e ; stem3 ( zz ( : , 1 ) , zz ( : , 2 ) ,w ( : ) , ‘ f i l l ‘ ) ; %Show the dominant e i g e n v e ct o r ←-
( d ( 1) =1) meant to r oughly e sti m at e i n v a r i a n t d en si t y
45 %
46 w2=v ( : , 1 ) ; [ i ,ww]= fi n d (w2>0) ;
47 f i g u r e ; stem3 ( zz ( i , 1 ) , zz ( i , 2 ) ,w2( i ) , ‘ r ‘ , ‘ f i l l ‘ ) %See Fig . 1 1 . 2 .
48
49 %% Produce a r e v e r s i b l e Markov Chain R
50 P=A;
51 [ v , d]= e i g s (A‘ , 1 ) ;
52 N=s i z e (A, 1 ) ;
Assignment № 8 Page 5
53 f o r i =1:N
54 f o r j =1:N
55 Phat ( i , j )=v ( j , 1 ) *P( j , i ) /v ( i , 1 ) ;
56 end
57 end
58 R=(P+Phat ) / 2;
59
60 %%
61 %P a r t i t i o n
62 [w, l ]= e i g s (R, 4 ) ;
63 f i g u r e ; p l o t (w( : , 2 ) ) %See Fig . 1 1.3
64 c =0; eps =0.00 5 ;
65 [ i ,ww]= fi n d (w ( : , 2 )>c ) ; [ i i ,ww]= f i n d (w ( : , 2 )<=c ) ; [ i i i ,ww]= f i n d ( abs (w ( : , 2 ) )<←-
eps ) ;
66
67 %% These commands pl ot t he r e s u l t i n g p a r t i t i o n in phase space .
68 % See Fig . 1 1 . 4 .
69 f i g u r e ; patch ( ‘ f a c e s ‘ , dt ( l l ( i ) , : ) , ‘ v e r t i c e s ‘ , dt .X, ‘ FaceColor ‘ , ‘ r ‘ ) ; ←-
hold on ;
70 patch ( ‘ f ac e s ‘ , dt ( l l ( i i ) , : ) , ‘ v e r t i c e s ‘ , dt .X, ‘ FaceColor ‘ , ‘ b ‘ ) ;
71 patch ( ‘ f ac e s ‘ , dt ( l l ( i i i ) , : ) , ‘ v e r t i c e s ‘ , dt .X, ‘ FaceColor ‘ , ‘ k ‘ ) ;
Listing 4: Matlab code – TransitionMatrix.m
1 %%%%%%%%%%%%%%%
2 %% by Erik B o l l t
3 %%%Build an Ula m - Galer kin ‘ s Matrix Based on a Test Orb it V i s i t i n g ←-
T r i a n g l e s
4 %%%of a T e s s e l at i o n
5 %%%%%%%%%%%%%%%
6 % Input :
7 % z - Test o r b i t i s n x 2 f o r an n - i t e r a t e o r b i t sample
8 % h - i s the s i d e l eng t h o f the t r i a n g l e e d ges which s h are a r i g h t
9 % an g l e
10 % ax , bx , ay , by - the low and high ends o f a box bounding data
11 % Output :
12 % dt - a DelaunayTri
13 %%
14
15 f u n c t i o n [ dt , l l ,A, zz ]= Tran s i t ionM a t rix ( z , h , ax , bx , ay , by )
16
17 %%
18
19 %low = -2; high =2; low =0; high =1;
20 [ X1, X2] = ndgrid ( ax : h : bx , ay : h : by ) ;
21 [m, n]= s i z e (X1) ;
Assignment № 8 Page 6
22 x1=reshape (X1,m*n , 1 ) ; x2=reshape (X2 ,m*n , 1 ) ;
23
24 %Formulate Delaunay T r ia n gu la t io n o f r e g i on
25 dt= DelaunayTri ( [ x1 x2 ] ) ; %See Matlab Su b r o u tine DelauneyTri f o r input /←-
output i nfo r ma t io n
26 %dt i s the t r i a n g u l a t i o n c l a s s
27 %
28 % where
29 %
30 %dt . dt . T ri a ng u la t io n i s an m1 by 3 array o f
31 %i n t e g e r s l a b e l l i n g the v e rte x c orn e r numbers
32 %o f the t r i a n g l e s
33 %
34 % and
35 %
36 %dt .X i s an m2 by 2 ar r a y o f r e a l numbers
37 %d e f i n i n g p o s i t i o n s
38
39 %t r i p l o t ( dt ) %Opti onal p l o t command o f t h i s t r i a n g u l a t i o n
40
41 %%
42
43
44 %Count number of o r b i t p oi n t s i n z which cause a t r i a n g l e to be counted as
45 %occupied ( and o th e rw i se a t r i a n g l e i s not counted as i t i s empty u n t i l
46 %obser v e d o c c u p i e d
47
48 n o ttru e =0;
49 w h i le ( n ottrue <1)
50 no t true =1;
51 SI = p o i ntL o cat i on ( dt , z ) ; %Matlab command , Locate the simplex element ←-
in dt
52 %c on ta i ni ng the s p e c i f i e d l o c a t i o n s o f each
53 %o f the e l e m e n t s o f the a r r a y z o f o r b i t
54 %samples
55 %
56 l=unique ( SI ) ; %Matlab s ubro u t ing : Count the number o f unique ←-
i n s t a n c e s i n SI
57 k=1;
58 wh i l e ( i snan ( l ( k ) )<1&&k<l eng t h ( l ) ) %Using Matlab su b rou t in e True f o r ←-
Not - a - Numb er
59 [ i i , j ]= f i n d ( SI==l ( k ) ) ; %C o l l e c t t h o se l o c a t i o n s c o rre s pon d ing to ←-
each unique l .
60 cnt ( k )=sum ( j ) ;
61 k=k+1;
62 end
63 k=k - 1 ; l l =l ( 1 : k ) ; cnt=reshape ( cnt , s i z e ( l l ) ) ;
Assignment № 8 Page 7
64 end
65 %%
66
67 % P l o t t hose sim p l e x e l e m e n t s o f the dt which a r e occupied by an o r b i t
68 % i t e r a t e o f z - dt ( l l , : ) a r e t h o se occupied simplex elements
69 patch ( ‘ f ac e s ‘ , dt ( l l , : ) , ‘ v e r t i c e s ‘ , dt . X, ‘ FaceColor ‘ , ‘ r ‘ ) ;
70 N=l eng t h ( l l ) ;
71
72 %%
73
74 %Translateba c k from l l back to phase space p o s i t i o n s z by usi n g the c en te r
75 %p o s i t i o n s with Matlab s u b ro u tin e "mean"
76 zz=ze r o s (N, 2 ) ;
77 f o r i =1:N
78 zz ( i , : ) =mean ( [ dt .X( dt . T r ia n gu l at io n ( l l ( i ) , 1 ) , : )
79 dt .X( dt . T r ia n gu la t io n ( l l ( i ) , 2 ) , : )
80 dt .X( dt . T r ia n gu la t io n ( l l ( i ) , 3 ) , : ) ] ) ;
81 end
82
83 %%Now b ui l d the t r a n s i t i o n matrix A
84 %
85 %So that A( i , j )>0 i f f t h er e i s an element in s i mplex element i such t h a t ←-
t h er e i s an i t e r a t e z ( k , : ) and t h a t
86 % the next i t e r a t e , z ( k +1 , : ) , t r a n s i t i o n s to s i m p l e x element j
87 A=ze r o s (N,N) ;
88 f o r i =1:( le n gt h ( SI ) - 1 )
89 i i = f i n d ( l l==SI ( i ) ) ;
90 j j =f i n d ( l l==SI ( i +1) ) ;
91 %[ i i j j s i z e (A) ]
92 A( i i , j j )=A( i i , j j ) +1;
93 end
94 %Now make A i n t o a s t o c h a s t i c matrix by row n orm a liz i ng
95 f o r i =1:N
96 q=l eng t h ( f i n d ( SI==l l ( i ) ) ) ;
97 A( i , : ) =A( i , : ) . / q ;
98 A( i , : ) =A( i , : ) . / sum(A( i , : ) ) ;
99 end
Question 2:
In Question 2, the task is to investigate the almost invariant sets of the standard map as the
critical parameter, k, is stepped from 0.9 to 1.5. The standard map is given by (2), and is known
to be non-linear (for k 6= 0) and chaotic. The graphs in Figure 2 show similar properties, as com-
pared with the graphs in Figure 3. It should be noted that for a value of k = 0 the standard map
Assignment № 8 Page 8
is linear, and increases in non-linearity as k is increased. The graphs were created with code
from Question 1 (standard.m, runnerSimpleTriCover.m, and TransitionMatrix.m). Accord-
ing to the reference material A Glance at the Standard Map by Ryan Tobin, as the non-linearity
factor, k, is increased, the nonintegrable components of the standard map become more ap-
parent. Also, the final invariant torus is terminated when the irrationality of the winding number
corresponding to the last quasiperiodic orbit is farthest from rational, which occurs at the golden
mean (approximately 1.61). This occurs at k > 0.97 and can be seen by comparing the graphs
in the left column of Figures 2 and 3, which include the standard map from k = 0.9 to k = 1.5.
In summary, for each different value of k, the codes work together to produce an Ulam-Galerkin
estimate of the Frobenius-Perron operator. The major steps of the code are as follows:
Step 1: Create a test orbit of the standard map (standard.m and runnerSimpleTriCover.m)
Step 2: Perform a Delaunay triangulation in the specified domain for the test orbit for a set k-value
(runnerSimpleTriCover.m and built-in function DelaunayTri.m)
Step 3: Prune the initial triangulation to a smaller collection of triangles which include only those
triangles which have been visited by the test orbit. (TransitionMatrix.m)
Step 4: Form stochastic matrix, A (which is an M × M matrix where M is the smaller number of
triangles), by using the Ulam-Galerkin estimate of the Frobenius-Perron operator. Transi-
tionMatrix.m)
The Frobenius-Perron operator can be defined as an infinitely large stochastic matrix acting
on an infinite-dimensional linear space, however this representation is not particularly useful in
practice. Instead, a finite-rank approximation of the operator is obtained through the process
described above. The matrix approximation we are after is given by,
P
i,j
=
m(B
i
∩ F
−1
(B
j
))
m(B
i
)
, (3)
where m denotes the Lebesgue measure on M. However, since we are working with a test
orbit (4) then becomes,
P
i,j
≈
#({x
k
|x
k
∈ B
i
and F (x
k
) ∈ B
j
})
#({x
k
∈ B
i
})
, (4)
Assignment № 8 Page 9
(a) A Delaunay triangulation for the test
orbit, k = 0.9
(b) The eigenvector v
2
for k = 0.9 (c) The thresholded weakly invariant sets
of the test orbit for k = 0.9
(d) A Delaunay triangulation for the test
orbit, k = 1.0
(e) The eigenvector v
2
for k = 1.0 (f) The thresholded weakly invariant sets
of the test orbit for k = 1.0
(g) A Delaunay triangulation for the test
orbit, k = 1.1
(h) The eigenvector v
2
for k = 1.1 (i) The thresholded weakly invariant sets
of the test orbit for k = 1.1
Figure 2: The images in the left column represent the Delaunay triangulation of the standard
map. The blue dots represent the test orbit, while the red triangles are the grid elements visited
by the test orbit. The images in the middle column are the second eigenvector, v
2
, of the re-
versible Markov chain. Lastly, the images in the right column are the almost invariant sets of the
standard map test orbit on the pruned tessellation (shown in the left column). The eigenvector
v
2
can be thresholded to produce the tessellation elements corresponding to weakly transitive
comnponents (red and blue), with the triangles on the boundary (black).
Assignment № 8 Page 10
(a) A Delaunay triangulation for the test
orbit, k = 1.2
(b) The eigenvector v
2
for k = 1.2 (c) The thresholded weakly invariant sets
of the test orbit for k = 1.2
(d) A Delaunay triangulation for the test
orbit, k = 1.4
(e) The eigenvector v
2
for k = 1.4 (f) The thresholded weakly invariant sets
of the test orbit for k = 1.4
(g) A Delaunay triangulation for the test
orbit, k = 1.5
(h) The eigenvector v
2
for k = 1.5 (i) The thresholded weakly invariant sets
of the test orbit for k = 1.5
Figure 3: The images here represent the same concepts as in Figure 2, but for values of k =
{1.2, 1.4, 1.5}. It can be seen here in the middle column of images of eigenvector v
2
that the
higher the value of k, the more the oscillation is present between sets (or sets). It can also
be seen that the higher k-values produce more chaotic behaviour, as is expected. (Note that
graphs for k = 1.3 have been omitted simply to save space, since their behaviour is similar to
the other graphs present)
Assignment № 8 Page 11
Question 3:
In this section the task is to apply the same methods in Part 1 to the Rossler Attractor (instead
of the Lorenz attractor). The Matlab code will attempt to create an outer cover of the chaotic
attractor. The Rossler Equations are given in (5), with parameters a = 0.2, b = 0.2, and c = 8.0.
dx
dt
= −y − z
dy
dt
= x + ay
dz
dt
= b + z(x − c)
(5)
The Rossler attractor is modeled by the code rossler.m, which makes use of Matlab‘s
ode45 to solve a set of ordinary differential equations. The data from rossler.m is then used
by runnerRosslerCover.m to create a cover of the attractor (similar to runnerLorenzCover.m
in Question 1). The results, using the same Delaunay Triangulation parameters as in Question
1, are given by Figure 4.
(a) Angle 1, Rossler Attractor (b) Angle 2, Rossler Attractor (c) Angle 3, Rossler Attractor
(d) Angle 1, Rossler Attractor outer cover
h = 2.5
(e) Angle 2, Rossler Attractor outer cover
h = 2.5
(f) Angle 3, Rossler Attractor outer cover
h = 2.5
Figure 4: The figure depicts the accuracy of computing the outer cover of the Rossler attractor
using Delaunay triangulation and basic grid element size h = 2.5. It is clear from the graphs
that the grid elements are too large to give a proper representation of the outer cover, causing
the figures on the bottom row to appear to be a very low representation of the attractor. This
can be fixed by adjusting the h parameter, and the results of doing so are represented in Figure
5.
Assignment № 8 Page 12
Listing 5: Matlab code – rossler.m
1 %% Create Data f o r R o s s l e r Equations
2 c l e a r ;
3 a = 0 . 2 ; b = 0 . 2 ; c = 8 . 0 ; % s e t parame ters f o r R os sl e r a t t r a c t o r
4 t = 0 : . 0 0 0 1 : 1 0 0 ; % s e t time s t ep and i n t e r v a l
5
6 x_init = 1 ; y_init = 1 ; z _ ini t = 1 ; % s e t i n i t i a l c o nd i t i o n s
7
8 % R os s l e r ODE
9 f = @( t , x ) [ - x (2 ) - x (3 ) ; x (1 )+a*x ( 2 ) ; b+x ( 1) *x ( 3) - c*x ( 3) ] ;
10
11 [ t1 , x1 ] = ode45 ( f , t , [ x _ i n i t y _ i n i t z _in i t ] ) ;
12
13 % p lo t r e s u l t s
14 f i g u r e
15 h old on
16 g r i d on
17 g r i d minor
18 t i t l e ( ‘ R o s sl er Equations ‘ )
19 p lo t 3 ( x1 ( : , 1 ) , x1 ( : , 2 ) , x1 ( : , 3 ) , ‘ k - ‘ )
20
21 % Save R os s l e r e q uat i on s data
22 X = x1 ;
23 X_data = X‘ ;
24 f il e na m e = ‘ r o s s l e r . mat ‘ ;
25 sav e ( fi l e n a me , ‘X‘ , ‘ X_data ‘ ) ;
It can be seen that in the Figure 4, the tessellation size is not low enough to give a very
descriptive cover for this particular attractor. This can be fixed by adjusting the size of the
elements (decreasing the tessellation size). The results from this adjustment can be seen in
Figure 5. The desired accuracy can be achieved through this adjustment ( tessellation size h),
with the caveat that the size of the sparse matrix, A, scales in size according to
1
h
d
for domain
dimension d (meaning more computation time and resources).
Assignment № 8 Page 13
(a) Angle 1, Rossler Attractor outer cover h = 1.5 (b) Angle 1, Rossler Attractor outer cover h = 0.5
(c) Angle 2, Rossler Attractor outer cover h = 1.5 (d) Angle 2, Rossler Attractor outer cover h = 0.5
(e) Angle 3, Rossler Attractor outer cover h = 1.5 (f) Angle 3, Rossler Attractor outer cover h = 0.5
Figure 5: Results of varying the parameter h in the cover computation. The outer cover of the
Rossler attractor becomes more clear as h is decreased, however this requires more compu-
tation resources. The h parameter corresponds to the tessellation size. The solution with the
best accuracy (of the tested h values) is when h = 0.5. This attractor required more resolution
to achieve results that are visibly satisfying than the Lorenz attractor from Question 1.
Assignment № 8 Page 14
```