Characterization of attractors in a model for boundary-driven nonlinear optical waveguide arrays with disorder, gain and damping

(1)

Characterization of attractors in a model for

boundary-driven nonlinear optical waveguide arrays with

disorder, gain and damping

Felix Faber

2013-05-17

Abstract

The purpose of this thesis is to study the effects of gain and damping on a nonlin-ear waveguide array with a strong disorder that is driven in the first site, and try to find regimes which have stable stationary solutions. This has been done with a modi-fied DNLS (Discrete nonlinear Schr¨odinger equation). Stable stationary solutions were mainly found when the damping was stronger than the gain, but some stable stationary regimes were also found when the gain was stronger than the damping.

(2)

1 Introduction

Studies have shown that linear waveguide arrays with randomly distributed potentials have a stronger localization in its solutions than linear waveguide arrays without random potentials. This is due to the so called Anderson localization [1].

It has also been shown that the Anderson localization effect can still survive in nonlinear waveguide arrays if the amplitude of the incident light is small enough [2]. This study is dedicated to see how the stability of the system is affected if the waveguide array also has gain and damping effects.

The DNLS (Discrete nonlinear Schr¨odinger equation) is a system of nonlinear cou-pled ODEs that in its simplest form looks like

id

dtθj+ γ|θj|

2_θ

j+ ε(θj+1+ θj−1) = 0 (1)

where j is a positive integer. This equation describes the motion of coupled an-harmonic oscillators [3]. The DNLS equation can be used to give a simplified model of a more complicated system, as a model for a Bose-Einstein condensate in certain settings [4] or optical waveguide arrays [5]. This thesis will be focused on the optical waveguide array model.

2 Theory

2.1 A short motivation why we use DNLS to describe the

system

The waveguide array is a system where the light is bounded in parallel 1D cavities along some z-axis in which we denote the complex light in the j : th cavity as θj. The cavities

are so close to each other so that the light can tunnel between the cavities, and the medium in which the wave is traveling has a refractive index that is dependent on the light intensity in the medium. This gives the waveguide a self-focusing property (Kerr effect)[5].

These two properties are effectively modeled by i_dzdθj+ ε(θj+1+ θj−1) + γ|θj|2θj = 0

[5]. Note that this looks like eq(1) but θj is a function of distance z in the waveguide

and not time t. The waveguide arrays studied in this work have a random potential in each cavity, which is why we introduce the term Vjθj where Vj is a real random

constant in the interval [−V, V ].

Since we want to see how light travels in a waveguide array affected by damping (e.g. heat loss) and gain (e.g resonance in lasers), we will add two terms which simulate the effects of gain and damping in the system. The damping effect will be simulated by ihθj

where h is a real constant and our gain to be simulated by −ig[θj]θj, g[θj] = _1+|θG j|2,

where G is a real constant. We have chosen g[θj] to decrease when the intensity |θj|2

in the waveguide increases [6, 7]. This has been done so that the total intensity of the system P

j|θj|2 does not tend to infinity (saturable). The waveguide is also driven in

one of the boundary waveguides by a wave with amplitude A and wave-number ω. Combining the above equations will yield the DNLS equation

i d

dzθj+ |θj|

2_θ

j− Vjθj+ i(h − g[θj])θj+ (θj+1+ θj−1) = Aj, (2)

where Aj = δ1,jAeiωz is the driving force with an amplitude A and wave-number ω.

This equation simulates the system in a time-independent regime and θ0 = θN = 0

(rigid boundary condition). Also note that γ and ε from eq(1) have been set equal to 1 since they can be removed by a rescaling of z and θj.

(4)

The term eiωz in eq(2) can be removed by transforming the system to a rotating frame θj → θjeiωz which yields:

id

dzθj + |θj|

2_θ

j− (ω + Vj)θj+ i(h − g[θj])θj+ (θj+1+ θj−1) = Aj. (3)

2.2 Uncoupled DNLS

We can neglect the coupling of the system in the first waveguide if g(θj), A, h, ω and

V >> 1, yielding in the first waveguide the equation id

dzθ + |θ|

2_{θ − (ω + V )θ + i(h − g[θ])θ = A.} ₍₄₎

An expression for the first site when _dzdθj = 0 can now be obtained for θj by solving

the following equations [8]:

|θ|2= A 2 (g(θ) − h)2_{+ (ω + V − |θ|}2₎2 (5) Re[θ] |θ|2 = |θ|2_{− ω − V} A (6) Im[θ] |θ|2 = g(θ) − h A . (7)

Eq(5) will have three solutions if the following inequalities are fulfilled

V + ω > 0 (8) | ω + V g[θ] − h| > √ 3, (9) " 1 + 9 g[θ] − h ω + V 2 − 27A 2 2(ω + V )3 #2 < " 1 − 3 g[θ] − h ω + V 2#3 (10) else it will have only one solution [8]. This is illustrated in Figure(1).

It is not possible to express eq(8,9,10) in terms of A only since we can not express |θ| as an explicit function of A. However, we can express A as an explicit function of |θ| from eq(5), so A2 _{in the inequality from eq(10) can be replaced by |θ|}2_{((g[θ] − h)}2₊

(ω + V − |θ|2)2) to give " 1 + 9 g[θ] − h ω + V 2 −27|θ| 2_{((g[θ] − h)}2_{+ (ω + V − |θ|}2₎2₎ 2(ω + V )3 #2 < " 1 − 3 g[θ] − h ω + V 2#3 . (11) When eq(5) has three solutions, it can be shown that one of the solutions will always be unstable [8]. If A = 0 eq(5) can be rewritten into

|θ|2((g(θ) − h)2+ (ω + V − |θ|2)2) = 0 (12) which will have the solutions

(5)

y x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.1 0.2 0.3 0.4 0.5 0.6

Figure 1: G has been set equal to zero which makes it possible to make the variable substitution, x = _ω+Vh and y = _Ah32.The shadowed area is the domain where eq(8,9,10)

has three solutions and outside the shadowed area is the domain where it only has one solution. Note that when h = 0, there will exist regimes where there are three solutions but this is not shown in the figure since x and y will be independent of ω

and V if h = 0.

The second solution in eq(13) will only exist if the following conditions are fulfilled:

G

h > 1, G/h − 1 = ω + V. ω will be a free variable since the system has no driving force

(A = 0) which means that the solution will oscillate with an angular wave-number ω. If we start with the initial condition |θ| = 0 and A = 0 when z = 0 and then increase A slowly, |θ| will follow that stationary solution (the linear solution) if it is stable towards a larger |θ| until a critical driving strength where the solution becomes unstable and bifurcates towards some other solution. We will use

A = A0(1 − e− z

τ) (14)

where τ is a real number between 100 to 1000.

2.3 Stationary solutions and turning point

If _dzdθj = 0 (stationary solution), we can rewrite eq(3) into

(ω + Vj− |θj|2)θj+ i(g[θj] − h)θj− θj+1= θj−1 (15)

for j 6= 1 and

−(ω + V1− |θ1|2)θ1− i(g[θ1] − h)θ1+ θ2 = A. (16)

A manifold with all stationary solutions for A can be found by assuming an initial condition (θN −1, θN) = (ϕ, 0) with N ' 100 and around 1000 to 10000 values of ϕ with

positive values and an order of magnitude around 10−13, then iterating backwards from j = N . An example of such manifold can be seen in Figure(2).

Note that if ϕ is replaced with ϕeiα where α is a real constant, iterating eq(15) backwards yields θj = θ0jeiα where0 denotes the case where α = 0. Then we can insert

θ1 and θ2 in eq(16) to yield A = A0eiα. This means that we can find a real solution for

(6)

A turning point (TP) of A is a local minimum or maximum of |A| and we will denote it as A[x] if it is plotted against a variable x. Each TP of A should correspond to a threshold-like behavior in the system [9].

0.0 0.2 0.4 0.6 0.8 1.0 1.2ÈAÈ 0.0 0.2 0.4 0.6 0.8 1.0 Èq1È 0.0 0.1 0.2 0.3 0.4ÈAÈ 0.00 0.05 0.10 0.15 0.20 Èq1È

Figure 2: The figures show a manifold with stationary solutions for a random DNLS with |θ1| plotted against |A|, the constants have been set equal to G = h = 0, ω = 0.2

and V = 1.9 and the length of the chain is 100. The left figure shows a bigger part of the manifold whereas the figure to the right is focused on the first turning point for

|A|.

3 Simulations etc

3.1 Total norm, participation number and lyapunov

ex-ponent

Total norm and participation number are tools that are useful for seeing how the system behaves as a whole.

The total norm

η =X

j

|θj|2 (17)

measures the amount of excitation of the system, which is the total power in optics [5]. The participation number

P = η 2 P j|θj|4 = ( P j|θj|2)2 P j|θj|4 (18) gives a good indication of the spreading in the system. The intensity is localized around a few sites if P is small and spread out on many sites if P is large.

The Lyapunov exponent is a real quantity that describes the separation of two trajectories with an infinitesimally small initial separation. A system with dimension N will have N Lyapunov exponents [10], this set of Lyapunov exponents is called a Lyapunov spectrum. The largest Lyapunov exponent of the Lyapunov spectrum determines stability of the system, which in most cases can be found using eq(19).

λ = lim z→∞δθ[0]→0lim 1 zln| δθ[z] δθ[0]| (19)

where δθ[z] is the difference between two trajectories of θ1 in eq(3) that have an initial

separation δθ[0]. This means that that δθ[z] = θ1[z] − θ01[z] where θ01[z] is the trajectory

(7)

by the dynamics in the first site, and it is hard to find a way to calculate the Lyapunov exponent for the whole system. This is why I have chosen to calculate the Lyapunov exponent in the first site only.

Since we can not do infinitely long simulations, or let δθ[0] go to zero, we have calculated the Lyapunov exponent as

λ = 1 z X n ln|δθ[zn+1] δθ[zn] | (20)

where zn = nδz, δz is a real constant and δθ[zn] is set equal to 10−5eia where a is a

random real constant in the interval [0, 2π). ln|δθ[zn+1]

δθ[zn] | is then calculated after a fixed

zn+1and δθ[zn+1] is then rescaled to 10−5eiaand the same procedure is done again. The

average value of all ln|δθ[zn+1]

δθ[zn] | exponents (eq(20)) should yield a good approximation of

the system Lyapunov exponent. The length between the rescaling of δθ[zn] has been set

equal to one in the simulations. The Lyapunov exponent is always calculated between z = 6000 and z = 10000 since we want to know the Lyapunov exponent of the system when A ' A0.

There have been some problems with calculating the Lyapunov exponent. The simulations take up a lot of RAM (around 14 Gb) and it appears to sometimes be negative in regimes where the dynamics appears to be chaotic but it still appears to reflect some of the dynamics of the system. There is not enough time in this study to see if there is something wrong with the code or if this is actually the case. The exact codes for the simulations and the values of the random chain Vj that was used will be

attached in the appendix.

All code for simulations were written in Mathematica. The disorder V has been set equal to 1.9, τ = 1000 and one set of values for the random chain Vj have been used

for all simulations since there was not enough time to vary all the parameters in this study. The time taken for calculating the Lyapunov exponent, the total norm and the participation number varies between 2 and 10 hours (With an intel i7 k2600 processor) depending on the input values and the rest of the simulations took around 5-20 minutes depending on the input values.

3.2 Without gain or damping

The TP of an infinitely long random DNLS A[ω] without any gain or damping (h = G = 0) is thought to yield an upper semi-continuous function [2], as seen in Figure(3) where the discontinuities correspond to resonances in the random sites and the strong resonances come from the sites close to the driving site and the weak come from sites located far away from the driving site.

If the driving strength is far below the first TP of A and if A is increased sufficiently slowly (eq(14)), the solution will remain localized around the driving site and will essentially be stationary with some small periodic oscillations [2].

If the driving strength is just below the first turning point, the small oscillations from when A were far below the first TP will be weakly chaotic and the system will at some z change its dynamics, it will still be bounded but with larger oscillations. The solution will then again change its dynamics at an even greater z, jumping into a continuously spreading and strongly chaotic regime [2].

The solution will jump into the weakly chaotic state as soon as A passes the first TP [2].

More detailed information about the system without gain or damping can be ob-tained by reading [2].

(8)

-1 1 2 3 4 5 Ω 1 2 3 4 5 A

Figure 3: The figure shows the first TP of the manifold of the stationary solutions for a random DNLS plotted against ω with h = G = 0, V = 1.9, n = 100. The figure is

plotted with a resolution 10−3 (1000 points per unit of ω).

3.3 With only damping

Adding only damping to the system appears to negate the discontinuities caused by the resonances in the sites in A[ω] which we can see in Figure(4); this effect increases as the damping is increased.

-1 1 2 3 4 5 Ω 1 2 3 4 5 A -1 1 2 3 4 5 Ω 1 2 3 4 5 A

Figure 4: The figures show the first TP of the manifold of the stationary solutions for a random DNLS plotted against ω with G = 0, V = 1.9, n = 100 in both figures and h = 0.01 in the left figure and h = 0.1 in the right. Both figures are plotted with the

same resolution as the plot in Figure(3).

If we look at how the first TP A[h] Figure(5) behaves in one of the ”resonance valleys” from the first TP A[ω], then we can see that h has a relatively small effect on the TP until a specific value of h when it jumps to a significantly higher value.

This behavior can be explained if we look at the manifold itself Figure(6) where we can see that the first TP disappears so that the system can follow the linear solution to the next TP.

In Figure(7) we can see to the left a plot of the Lyapunov exponent of the system with respect to h with an A0 that is above the first TP before the jump at h ' 0.031

and below the first TP after the jump, and a plot of the Lyapunov exponent of the system with respect to h with an A0 that is above the first TP before and after the

(9)

0.02 0.04 0.06 0.08 0.10 h 0.2 0.4 0.6 0.8 1.0 A

Figure 5: The figure shows the first TP of the manifold of the stationary solutions for a random DNLS plotted against h with G = 0, V = 1.9, n = 100 and ω = 0.2. We can clearly see a discontinuity at h ' 0.031 where A[h] jumps to a significantly higher

value. 0.430 0.435 0.440 0.445 0.450ÈAÈ 0.162 0.164 0.166 0.168 0.170 0.172 0.174 ÈΘ1È 0.430 0.435 0.440 0.445 0.450ÈAÈ 0.162 0.164 0.166 0.168 0.170 0.172 0.174 ÈΘ1È

Figure 6: Two plots of the manifold of stationary solutions to eq(3) with G = 0, V = 1.9, ω = 0.2, n = 100 with h = 0.031 to the left and h = 0.032 to the right. One

can see that the TP exists in the left manifold but not in the right, this explains why there is a discontinuity if you plot A[h].

Note that the Lyapunov exponent is negative in both cases even for very small values of h, though it is smaller when A0 is above the first TP before and after the jump in

Figure(5). λ does not appear to be affected when the system passes the jump in the first TP at h ' 0.031 and the solutions which can be seen in Figure(8) are practically the same before and after the jump when A0 is below the first TP after the jump at

h ' 0.031.

Figure(9) shows the first site plotted against h, we can see that the system appears to smoothly converge when A0 is above the first TP before the jump and below after

the jump, whereas the system appears to make a more or less discrete jump to its final value when A0 is above the first TP before and after the jump. We can also see that

the value that the first site converges to is significantly smaller when A0 is above the

first TP before the jump and below after the jump than the value that the first site jumps to when A0 is above the first TP before and after the jump.

The participation number at z = 10000 is plotted against different values of h in Figure(10), with A0 above the first TP before and below after the jump plotted to the

(10)

0.01 0.02 0.03 0.04 0.05 0.06 0.07 h -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.05 l 0.01 0.02 0.03 0.04 0.05 0.06 0.07 h -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 l

Figure 7: Lyapunov exponent of eq(3) plotted against h with G = 0, V = 1.9, n = 102, ω = 0.2 and A0 = 1 in the left figure and A0 = 1.2 in the right figure. Note

that both plots only consist of 80 interpolated points which means that the resolution is very low, but λ is not positive when h goes to zero and A0 = 1.2.

0 50 100 n 0 5000 10 000 z 0.0 0.1 0.2 0.3 ÈqÈ^2 0 50 100 n 0 5000 10 000 z 0.0 0.1 0.2 0.3 ÈqÈ^2

Figure 8: Simulation of eq(3) before the jump of A[h] in Figure(5) to the left (h = 0.03), and after the jump to the right (h = 0.04). The other values are A0 = 1,

G = 0, V = 1.9, ω = 0.2 and n = 102.

can see that P is fairly small in both cases, except when h . 0.0008 where P starts to increase drastically.

Also note that the system is more localized when A0 is above the first TP before

and after the jump than when A0 is above the first TP before the jump and below

after. This can be seen even clearer if we plot the participation number at z = 10000 against different values of A0 with a value of h which lies after the jump in the first TP

of A. In Figure(11), a strong decrease of P can be seen when A0 passes the TP of A at

A ' 1.15.

There is a steep increase in the Lyapunov exponent when h . 0.0008 in both cases, and will become positive if A0is above the first TP before the jump and below after, but

not if A0is above the first TP before and after the jump. The reason why the Lyapunov

exponent is not positive when A0is above the first TP before and after the jump may be

due to the fact that I did not calculate the Lyapunov exponent over a sufficiently long length. One good indication that the Lyapunov exponent is positive when h is small is that the first site starts to oscillate at small h in both cases in Figure(9), however there was not enough time to try to calculate the Lyapunov exponent over longer lengths.

A comparison of the dynamics of the system before and after the steep increase in Figure(7) when A0 is above the first TP before and below after the jump at h = 0.031

(11)

0.02 0.04 0.06 h 0 5000 10 000 z 0.00 0.05 0.10 0.15 0.20 ÈqÈ^2 0.00 0.02 0.04 0.06 h 0 5000 10 000 z 0 2 4 ÈqÈ^2

Figure 9: A plot of the first site for different values of h with A0 = 1 to the left and

A0 = 1.2 to the right, other values are ω = 0.2, G = 0, V = 1.9 and n = 102.

0.00 0.01 0.02 0.03 0.04 0.05 0.06h 0 10 20 30 40 50 60 70 P 0.000 0.05 0.10 0.15 0.20 0.25 0.30h 10 20 30 40 50 60 70 P

Figure 10: Participation number of eq(3) plotted against h with G = 0, V = 1.9, ω = 0.2, n = 102, A0 = 1 in the left figure and A0 = 1.2 in the right figure.

similar to the continuously spreading and strongly chaotic regime when A had passed the first TP in [2], if h . 0.0008 but not if h & 0.0008. The system seems to jump into this regime as well when h . 0.0008 if A0 is above the first TP before and after the

jump, this is not shown though.

3.4 With gain and damping

Adding only gain negates the resonances in the first TP A[ω] system similarly to adding only damping.

If G = h and G, h is not too large (G = h ∼ 0 − 0.5) then the first TP A[ω] will remain very similar to the case without gain or damping, as seen in Figures(3, 13). This suggests that g = _1+|θG

j|2 ' h which means that |θj| must be small which in turn

means that η should be relatively small. If we plot the first TP of A[h] with a nonzero gain or A[g] with a nonzero damping Figure(14) (A[g] is not shown), we can see that there is a zone around G = h (0.06 . h . 0.13) in Figure(14) where A[h] or A[g] is significantly lower.

(12)

0.0 0.5 1.0 1.5 2.0 A0 0 2 4 6 8 10 P

Figure 11: Participation number of eq(3) plotted against A0 with h = 0.05, G = 0,

V = 1.9, ω = 0.2 and n = 102. 0 50 100 n 0 5000 10 000 z 0 2 4 ÈqÈ^2 0 50 100 n 0 5000 10 000 z 0.0 0.2 0.4 ÈqÈ^2

Figure 12: Simulation of eq(3) with (h = 0.0007) to the left and (h = 0.0008) to the right. The other values are A0 = 1, G = 0, V = 1.9, ω = 0.2 and n = 102.

In Figure(15) we can see to the left a plot of the Lyapunov exponent of the system with respect to h with an A0 that is above the first TP between the jumps (0.06 .

h . 0.13 in Figure(14)) and below the first TP when not, and to the right a plot of the Lyapunov exponent of the system with respect to h with an A0 above the first TP

for all h. We can see that there is practically no difference between λ in the two cases when h . 0.08, but somewhere between h = 0.08 and h = 0.09 λ jumps to a higher value when A0 is above the first TP between the jumps and below the first TP when

not, which it does not when A0 is above the first TP for all h.

We can further see that the dynamics is very similar in the two cases when h . 0.09 and then deviates if we look at the first site plotted against h, which is seen in Figure(16).

The dynamics when 0.09 . h is very similar to the case without any gain which we can see in Figure(9). The first site (the right figure of Figure(16)) does a jump that is very similar to the case with no gain and an A0 that is above the first TP for all h.

Though it can not be seen in Figure(16), the dynamics is also basically the same as in the left graph in Figure(9) when A0 is above the first TP before and below after the

jump and 0.09 . h.

(13)

-1 1 2 3 4 5 Ω 1 2 3 4 5 A -1 1 2 3 4 5 Ω 1 2 3 4 5 A

Figure 13: The figure shows the first TP of the manifold of the stationary solutions for a random DNLS plotted against ω with h = G = 0.1 to the left and h = G = 0.5 to the right. Other values are V = 1.9, n = 100. The resolution is the same as in

Figure(3). 0.00 0.05 0.10 0.15 0.20 h 0.2 0.4 0.6 0.8 1.0 A

Figure 14: The figure shows the first TP of the manifold of the stationary solutions for a random DNLS plotted against h with ω = 0.2, G = 0.1, V = 1.9, n = 100. The

resolution is 10−3.

some z = z0 jump to a significantly higher value when h . 0.09, which can be seen in

Figure(17).

P will go towards the number of sites N as h is decreased further, however there exist stable stationary solutions when h < G. There is one in the small dip in the Lyapunov exponent in Figure(15) at h ' 0.03 and the system appears to become stationary when h goes to zero. This can be seen in Figure(18) where we plot the difference between the maximum and minimum value of |θ| against h. We can see these two solutions in Figure(19).

In Figure(20) we see a plot of the real part against the imaginary part of the first site in the stationary/semi stationary solutions (all sites appear to have the same dynamics with different amplitude but it is not shown). We can see that the stationary solution at h ' 0.03 is practically constant whereas the other appears to have a motion that is bounded in a circle with a center that is not in the origin, but shifted slightly to the positive side on the Re[θ1]-axis. This may explain why λ in Figure(15) increases and

becomes positive but the intensity converges towards a constant value as h goes to zero, which means that the chaotic motion must only lie in the phase dynamics.

(14)

0.02 0.04 0.06 0.08 0.10 0.12 h -0.5 0.5 1.0 1.5 2.0 2.5 l 0.02 0.04 0.06 0.08 0.10 0.12 h -0.5 0.5 1.0 1.5 2.0 2.5 l

Figure 15: Lyapunov exponent of eq(3) plotted against h with A0 = 1 in the left figure

and A0 = 1.2 in the right. Other values are ω = 0.2, G = 0.1, V = 1.9 and n = 102.

Note that both plots only consist of 60 interpolated points which means that the resolution is very low.

0.05 0.10 h 0 5000 10 000 z 0 5 10 15 ÈqÈ^2 0.05 0.10 h 0 5000 10 000 z 0 5 10 15 20 ÈqÈ^2

Figure 16: A plot of the first site for different values of h with ω = 0.2, G = 0.1, V = 1.9 and n = 102 and A0 = 1 to the left and A0 = 1.2 to the right.

Similar stationary large amplitude solutions where found in [6] when the system had no disorder. These solutions amplitudes were found to be constant over all sites, which they are not when disorder is added, which we have seen in Figure(19).

A series of plots of the imaginary part of θ1 against the real part of θ1 for different

values of h will be attached in appendix so one can see how this manifold behaves in different regimes, and how it changes when h is varied.

4 Conclusion

Stationary solutions were found in the system when both gain and damping was added even when the driving force (A0) is above the turning point of the linear solution. The

damping seems to have a stabilizing effect on the system, making the stationary solution stable almost every time the damping is greater than the gain except in some cases when the damping is very small or when the gain is almost as strong as the damping.

The gain on the other hand seems to have a destabilizing effect on the stationary solutions, making the linear solution unstable when the gain in the system is greater than the damping. However two stable stationary solutions when the gain is greater than the damping were found, and these two solutions appear when the intensity have

(15)

0 50 100 n 0 5000 10 000 z 0.0 0.1 0.2 0.3 ÈqÈ^2 0 50 100 n 0 5000 10 000 z 0 1 2 ÈqÈ^2

Figure 17: Simulation of eq(3) with A0 = 1, G = 0.1, V = 1.9, ω = 0.2 and n = 102

and h = 0.095 to the left and h = 0.085 to the right.

spread out over all sites. One of them seems to come from some sort of resonance effect between the sites and is constant (_dzdθj = 0) and seems to only appear in a

small gap when the gain is in a certain relation with the damping. The other solution only has a constant intensity (_dzd|θj| = 0) and appears when the damping goes to

zero. The Lyapunov exponent is positive in this regime which can either be caused by something wrong in the code that calculates the Lyapunov exponent, or the phase of this solution has a chaotic dynamic. There was not enough time in this study to examine which one of the two cases it is but it needs to be done. Another thing that needs to be examined further is why the Lyapunov exponent is negative in regimes where the dynamics suggests that it is positive.

(16)

0.01 0.02 0.03 0.04 0.05 0.06 h 0.5 1.0 1.5 2.0 2.5 Èdq1È

Figure 18: The graph shows the difference between the maximum and minimum value of |θ1| plotted against h, with A0 = 1, ω = 0.2, G = 0.1, V = 1.9 and n = 102.This

graph has been plotted with 120 points and the smallest value of h is 10−3 . We can see that the difference appears to go toward zero when h goes to zero.

0 50 100 n 0 5000 10 000 z 0 2 4 ÈqÈ^2 0 50 100 n 0 5000 10 000 z 0 10 20 30 ÈqÈ^2

Figure 19: A plot of the system with A0 = 1, ω = 0.2, G = 0.1, V = 1.9 and n = 102

and h = 0.03 to the left and h = 0.003 to the right. If we look closely at the first site in the right figure, we see some small oscillations but these become smaller as h

(17)

1.447664 1.447669 1.447674Re@q1D 0.547582 0.547589 0.547596 Im@q1D -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D

Figure 20: A plot of the real part versus the imaginary part of the first site when 9500 < z < 10000 and with a sample rate equal to 1. A0 = 1, ω = 0.2, G = 0.1,

(18)

References

[1] K. Staliunas, Spatial solitons and Anderson localization, Phys. Rev, A 68 (2003), 013801 Page 1-7

[2] M. Johansson, G. Kopidakis, S. Lepri and S. Aubry, Transmission thresholds in time-periodically driven nonlinear disordered systems, EPL, Vol 86 (2009), 10009 [3] S Flach, A. V. Gorbach, Discrete breathers — Advances in theory and

applica-tions, Phys. Rep, Vol 467 (2008) Page 1-116

[4] R. Franzosi, R. Livi, G.-L. Oppo and A. Politi, Discrete breathers in Bose–Einstein condensates, Nonlinearity, Vol 24, (2011) Page R89-R122

[5] F. Lederer, G. I. Stegeman, D. N. Christodoulides, G. Assanto , M. Segev and Y. Silberberg, Discrete solitons in optics, Phys. Rep, Vol 463 (2008) Page 1-126 [6] J. E. Prilepsky, A. V. Yulin, M. Johansson and S. A. Derevyanko, Discrete solitons

in coupled active lasing cavities, Opt. Lett., Vol 37 (2012) Page 4600-4602 [7] S.K. Turitsyn et al, Random Distributed feedback fiber laser, Nature Photonics,

Vol 4 (2010) Page 231

[8] M. Rigo, G. Alber, F. Furtado and P. F. O’Mahony, Quantum-state diffusion model and the driven damped nonlinear oscillator, Phys. Rev A, Vol 55 (1997) Page 1665-1673

[9] P. Maniadis, G. Kopidakis and S. Aubry, Energy dissipation threshold and self-induced transparency in systems with discrete breathers, Physica D, Vol 216 (2006) Page 121-135

[10] G. Ohl´en, S. ˚Aberg, P. ¨Ostborn, Chaos, Division of Mathematical Physics LTH, Lund 2007

5 Appendix

A

Imaginary part of θ

1

against the real part of

θ

1

A series of plots of the imaginary part of θ1 against the real part of θ1 in the interval

9500 < z < 10000. There are a total of 120 graphs and they are chronologically ordered, starting with the lowest damping h = 0.003 and then increasing by 0.0003 between each graph. The other values are: G = 0.1,A0 = 1,ω = 0.2, V = 1.9 and n = 102.

(19)

Out[65]=: -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -4 -2 2 4 Re@q1D -4 -2 2 4 Im@q1D , -3-2-1 1 2 3 4Re@q1D -3 -2 -1 1 2 3 4 Im@q1D , -3-2-1 1 2 3 4Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3-2-1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3-2-1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , , , , Out[65]= -3-2-1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3-2-1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2-1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2-1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2-1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -3 -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 Im@q1D , -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 3 Im@q1D , -2 -1 1 2 3 Re@q1D -3 -2 -1 1 2 Im@q1D , , , , 2 Untitled-3 Out[65]= -2 -1 1 2 3Re@q1D -3 -2 -1 1 2 Im@q1D , -2 -1 1 2 3Re@q1D -3 -2 -1 1 2 Im@q1D , -2 -1 1 2 3Re@q1D -3 -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , Untitled-3 3 Out[65]= -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -2 -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , 4 Untitled-3

(20)

Out[65]= -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -1 1 2 Re@q1D -2 -1 1 2 Im@q1D , -1.0 -0.5 0.51.01.52.0 Re@q1D -2 -1 1 2 Im@q1D , -1.0 -0.5 0.5 1.0 1.52.0 Re@q1D -2.0 -1.5 -1.0 -0.5 0.5 1.0 1.5 Im@q1D , -1.0 -0.5 0.5 1.0 1.5 2.0 Re@q1D -2.0 -1.5 -1.0 -0.5 0.5 1.0 1.5 Im@q1D , -1.0 -0.5 0.5 1.01.5 2.0Re@q1D -2.0 -1.5 -1.0 -0.5 0.5 1.0 1.5 Im@q1D , -1.0 -0.5 0.5 1.0 1.5 2.0Re@q1D -2.0 -1.5 -1.0 -0.5 0.5 1.0 1.5 Im@q1D , -1.0 -0.5 0.5 1.0 1.5 2.0 Re@q1D -2.0 -1.5 -1.0 -0.5 0.5 1.0 1.5 Im@q1D , -0.5 0.5 1.0 1.5 2.0Re@q1D -2.0 -1.5 -1.0 -0.5 0.5 1.0 1.5 Im@q1D , -0.5 0.5 1.0 1.5 2.0Re@q1D -1.5 -1.0 -0.5 0.5 1.0 Im@q1D , 0.5 1.0 1.5 2.0Re@q1D -1.5 -1.0 -0.5 0.5 1.0 Im@q1D , 0.5 1.0 1.5 2.0Re@q1D -1.5 -1.0 -0.5 0.5 1.0 Im@q1D , 1.0 1.5 2.0Re@q1D -0.5 0.5 1.0 Im@q1D , , , , Untitled-3 5 Out[65]= 0.8 1.0 1.2 1.4 1.6 1.8 Re@q1D -0.2 0.2 0.4 0.6 0.8 1.0 1.2 Im@q1D , 0.8 1.0 1.21.4 1.6 1.8Re@q1D -0.2 0.2 0.4 0.6 0.8 1.0 Im@q1D , 1.0 1.2 1.4 1.6 1.8 Re@q1D 0.2 0.4 0.6 0.8 1.0 Im@q1D , 1.345091.345091.34510Re@q1D 0.707245 0.707250 0.707255 0.707260 Im@q1D , 1.402811.402811.402811.402821.402821.40282Re@q1D 0.626500 0.626505 0.626510 Im@q1D , 1.447671.447671.447671.447671.44767Re@q1D 0.547585 0.547590 0.547595 Im@q1D , 1.483541.483541.483551.483551.48355Re@q1D 0.470198 0.470200 0.470202 0.470204 0.470206 0.470208 0.470210 Im@q1D , 1.512221.512221.512221.51223Re@q1D 0.394220 0.394222 0.394224 0.394226 0.394228 Im@q1D , 1.534781.534781.534781.534781.534781.53478Re@q1D 0.319590 0.319592 0.319594 0.319596 0.319598 Im@q1D , 1.551941.551941.551951.551951.55195Re@q1D 0.246284 0.246285 0.246286 0.246287 0.246288 0.246289 0.246290 0.246291 Im@q1D , 1.564271.564271.564271.564271.56427Re@q1D 0.174301 0.174302 0.174303 0.174304 0.174305 Im@q1D , 1.572171.572171.572171.57217 Re@q1D 0.103649 0.103649 0.103650 0.103650 0.103651 0.103651 0.103652 Im@q1D , 1.575981.575981.575981.57598Re@q1D 0.0343526 0.0343528 0.0343530 0.0343532 0.0343534 0.0343536 Im@q1D , 1.575971.575971.575971.57597 Re@q1D -0.0335569 -0.0335568 -0.0335567 -0.0335566 -0.0335565 -0.0335564 Im@q1D , 1.572371.572371.572371.57237 Re@q1D -0.100037 -0.100037 -0.100036 -0.100036 -0.100035 Im@q1D , , , , 6 Untitled-3 Out[65]= 1.54 1.56 1.58 1.60 1.62Re@q1D -0.2 -0.1 0.1 Im@q1D , 1.52 1.54 1.56 1.58 1.60 1.62Re@q1D -0.3 -0.2 -0.1 Im@q1D , 1.50 1.52 1.54 1.56 1.58 1.60 1.62Re@q1D -0.3 -0.2 -0.1 0.0 Im@q1D , 1.52 1.54 1.56 1.58 1.60 Re@q1D -0.3 -0.2 -0.1 0.0 Im@q1D , 1.48 1.50 1.52 1.54 1.56 1.58 1.60Re@q1D -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 Im@q1D , 1.50 1.55 1.60Re@q1D -0.5 -0.4 -0.3 -0.2 -0.1 Im@q1D , 1.45 1.50 1.55 1.60Re@q1D -0.6 -0.5 -0.4 -0.3 -0.2 Im@q1D , 1.35 1.40 1.45 1.50 1.55 1.60Re@q1D -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 Im@q1D , 1.40 1.45 1.50 1.55 Re@q1D -0.6 -0.5 -0.4 -0.3 -0.2 Im@q1D , 1.45 1.50 1.55 Re@q1D -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 Im@q1D , 1.44 1.46 1.48 1.50 1.52 1.54 Re@q1D -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 Im@q1D , 1.40 1.45 1.50 Re@q1D -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 Im@q1D , 1.35 1.40 1.45 1.50 Re@q1D -0.7 -0.6 -0.5 -0.4 Im@q1D , 1.3 1.4 1.5 1.6 1.7 Re@q1D -0.8 -0.6 -0.4 -0.2 0.0 Im@q1D , 1.31.4 1.5 1.61.7 Re@q1D -0.8 -0.6 -0.4 -0.2 0.2 Im@q1D , Untitled-3 7 Out[65]= 1.2 1.31.41.51.61.7 Re@q1D -0.8 -0.6 -0.4 -0.2 0.2 Im@q1D , 1.21.3 1.41.51.61.7 Re@q1D -1.0 -0.8 -0.6 -0.4 -0.2 0.2 Im@q1D , 1.2 1.3 1.4 1.5 1.6 1.7 Re@q1D -1.0 -0.8 -0.6 -0.4 -0.2 0.2 0.4 Im@q1D , 1.1 1.2 1.3 1.4 1.5 1.6 1.7 Re@q1D -1.0 -0.8 -0.6 -0.4 -0.2 0.2 0.4 Im@q1D , 1.2 1.4 1.6 Re@q1D -1.0 -0.8 -0.6 -0.4 -0.2 0.2 0.4 Im@q1D , 1.2 1.4 1.6 1.8 Re@q1D -1.0 -0.8 -0.6 -0.4 -0.2 0.2 Im@q1D , 1.2 1.3 1.4 1.5 1.6 1.7 Re@q1D -1.0 -0.8 -0.6 -0.4 -0.2 0.2 Im@q1D , 1.21.31.41.51.61.7 Re@q1D -1.0 -0.8 -0.6 -0.4 -0.2 0.2 0.4 Im@q1D , 1.0 1.2 1.4 1.6 1.8Re@q1D -1.0 -0.5 0.5 Im@q1D , 1.0 1.2 1.4 1.6 1.8Re@q1D -1.0 -0.5 0.5 Im@q1D , 1.0 1.2 1.4 1.6 1.8Re@q1D -1.0 -0.5 Im@q1D , 1.0 1.2 1.4 1.6 1.8Re@q1D -1.0 -0.5 0.5 Im@q1D , 1.2 1.4 1.6 1.8 Re@q1D -1.0 -0.5 0.5 Im@q1D , 1.0 1.2 1.4 1.6 1.8 Re@q1D -1.0 -0.5 Im@q1D , 1.0 1.2 1.4 1.6 1.8 Re@q1D -1.0 -0.5 0.5 Im@q1D > 8 Untitled-3

(21)

B

Random chain

The list below is all random values used in the simulation. Always using the first value in the list (−0.485786) as V1, the second value (−0.66244) as V2 and so on.

-0.485786, -0.66244, -0.730974, 0.733407, -0.516704, 0.225218, 0.169968, -0.395415, -0.297793, -0.677518, -0.161885, -0.644241, -0.398672, 0.228784, 0.250667, 0.830526, 0.517702, -0.0741628, 0.213087, 0.364442, -0.616038, 0.361357, -0.568729, 0.877715, 0.374256, -0.0789293, 0.462926, -0.653213, 0.398076, 0.553102, 0.934611, -0.102473, -0.472244, 0.54991, 0.320828, -0.218046, -0.760284, 0.772081, -0.830088, 0.566398, -0.405706, -0.200326, -0.604278, 0.792555, 0.99921, -0.973142, 0.603608, -0.185167, 0.639737, -0.687197, 0.611642, -0.57459, -0.936145, -0.502369, -0.0567764, -0.247019, -0.493122, -0.977279, -0.195978, 0.580145, 0.737844, -0.923921, 0.983789, -0.437439, 0.407859, 0.971148, 0.138714, -0.422082, -0.390664, -0.513364, -0.876501, -0.850353, -0.495019, -0.0304584, 0.564389, -0.225581, 0.754904, 0.956463, 0.29576, -0.63506, 0.748796, 0.994665, -0.334891, -0.108553, -0.250815, 0.79525, 0.819362, 0.244421, 0.563012, -0.199316, 0.806498, 0.917213, 0.383193, 0.51981, 0.874409, -0.233505, 0.490223, -0.288882, 0.26837, -0.587346, -0.908257, -0.9587, 0.776594, 0.969251, 0.812041, -0.0999044, 0.920006, 0.219313, -0.780989, -0.411296, 0.0983765, -0.723132, -0.485891, 0.347257, -0.440121, 0.750995, -0.0425068, 0.613119, -0.152976, -0.304898, -0.82869, -0.831647, -0.76359, -0.754775, 0.612815, 0.280301, -0.740633, 0.697339, 0.0490698, 0.578603, 0.997044, -0.122689, -0.746765, -0.981923, 0.206192, -0.502883, 0.401822, 0.0904307, -0.752999, -0.652357, 0.265451, 0.915246, -0.859325, -0.313632, -0.597015, 0.585123, 0.274431, 0.391394, 0.857614, 0.57084, 0.0958453, 0.714885, 0.768803, -0.231486, -0.40044, -0.541678, 0.687465, 0.515846, 0.135539, 0.409874, 0.195364, -0.392189, 0.271666, -0.827333, -0.649693, 0.192318, -0.868866, -0.212508, -0.658037, -0.135848, -0.402126, -0.736202, -0.821544, -0.502946, 0.396093, 0.747395, 0.820428, 0.923992, 0.896479, 0.459592, -0.124916, -0.88087, 0.705875, 0.861678, -0.837992, 0.665797, 0.00186692

(22)

(23)

H*Parameter values*L ClearAll@"Global`*"D

H*I denna cell defineras konstanter för systemet*L H*d skalar avståndet mellan vågledare när den plottas*L d = 100;

H*B anger graden av oårdning i systemet*L B = 1.9;

H*A@tD anger drivkraften i första randsajten*L At = 1;

A@tD = At * H1 - Exp@-t • 1000DL;

H*g anger hur mycket "gain" systemet har*L g = 0.1;

w = 0.2;

H*h anger dämpning i systemet,h är beroende av f, f är kommer att fixeras till olika värden när

ekvationsystemet löses för att få en lista med

en lista amed lösningar för olika värden på h*L

h =Hf - 1L * fdif + hmin; H*n anger antalet vågledare

Hutom de som sitter i randen av systemetL,

totala antalet vågledare kommer vara n+2*L

n = 100;

H*Systemet kommer att simuleras mellan tiden 0 och time*L time = 10 000; fdif = 0.0003; hmin = 0.003; fmax = 120; line = 1; stör = 10 ^ - 5; LyMin = 0.95; steglengd = 1; steglengdlyd = 2; steglengdly = 1;

H*randomlist en lista med slumtpal som

används för att ge en oårdning i systemet*L

root = "D:\\Dropbox\\viktor\\";

randomlist = <<Hroot <> "randomlist.txt"L;

ProgressIndicator@Dynamic@iD, 80, time - steglengd<D

ProgressIndicator@Dynamic@kD, Dynamic@80, 5<DD

H*I denna cell löses ekvationsystemet*L nonlin = Evaluate@-g • H1 + Abs@uk@tDD^2L + hD;

H*U@tD är lista med de funktioner som ska lösas*L U@tD = Table@uk@tD, 8k, 1, n + 2<D;

H*bound1 beskriver ekvationen

för första vågledaren för systemet*L

bound1 =8I * u1'@tD + u2@tD + Abs@u1@tDD^2 * u1@tD

-HB * randomlist@@1DD + wL * u1@tD +

I * nonlin * u1@tD Š A@tD< •. k ® 1;

H*bound1 beskriver ekvationen

för första vågledaren för systemet*L

bound2 =8I * un+2'@tD + un+1@tD + Abs@un+2@tDD^2 * un+2@tD

-H B * randomlist@@n + 2DD + wL * un+2@tD +

I * nonlin * un+2@tD Š 0< •. k ® n + 2 ;

H*V@tD är en lista med ekvationer som beskriver alla vågledare utom de som är i

kanternaHsom beskrivs av bound1 och bound2L*L

V@tD = Table@I * uk'@tD + Huk+1@tD + uk-1@tDL +

Abs@uk@tDD^2 * uk@tD - HB * randomlist@@kDD + wL * uk@tD +

I * nonlin * uk@tD Š 0, 8k, 2, 1 + n<D;

H*The main simulation code*L

at0 = AbsoluteTime@D;

For@i = 0, i £ time - 20, i += 20,

2 Sim1.nb

Clear@KD;

If@i Š 0, K = Table@uk@iD Š 0, 8f, 1, fmax<, 8k, 1, n + 2<D,

K = Table@uk@iD Š lines@@f, kDD •. t ® i,

8f, 1, fmax<, 8k, 1, n + 2<DE;

Clear@linesD;

lines = Table@

NDSolve@8V@tD, bound1, bound2, K@@fDD<, U@tD, 8t, i, i + 20<, MaxSteps ® 900 000, InterpolationOrder ® 1, PrecisionGoal ® 12D, 8f, 1, fmax<D@@All, 1, All, 2DD ; If@ i ³ LyMin * time, at1 = AbsoluteTime@D; If@i Š time - 20,

slist = Table@k, 8k, i, i + 20 , steglengd<D,

slist = Table@k, 8k, i, i + 20 - steglengd, steglengd<DD; H*gLinesPart =

Table@lines@@f,kDD•.t®slist,8f,fmin,fmax<,8k,1,n+2<D;*L

gLinesPart = lines•. t ® slist;

tjoin1 = AbsoluteTime@D;

If@i Š LyMin * time, gLines = gLinesPart,

gLines = Table@Join@gLines@@f, kDD, gLinesPart@@f, kDD, 1D,

8f, 1, fmax<, 8k, 1, n + 2<DD;

tjoin2 = AbsoluteTime@D

H*Print@"köring: ",

i•100 +1,"; sampling: ",at2-at1," sekunder; join: ",tjoin2-tjoin1," sekunder" D;*L

D; D;

Sim1.nb 3

SetAccuracy@Hat3 - at0L • 60, 3D, " minuter."D; H*Intglines =

Table@ListInterpolation@Abs@gLines@@f,kDDD^2,80,time<,

InterpolationOrder®0D@tD,8f ,1,fmax<,8k,1,n+2<D;*L

H*This section plots the imaginary

part against the real part of theta*L

dLines = Table@gLines@@All, k, AllDD, 8k, 1, 1<D; Dlines = Table@gLines@@All, k, jDD, 8k, 1, 1<, 8j, 1, 500<D;

Dimensions@dLinesD

Dimensions@DlinesD

ImReInt = ListInterpolation@Table@

Max@Abs@Dlines@@1, All, jDDDD - Min@Abs@Dlines@@1, All, jDDDD,

8j, 1, fmax<D, 88h •. f ® 1, h •. f ® fmax<<D

Table@ListPlot@Table@

8Re@dLines@@1, j, nDDD, Im@dLines@@1, j, nDDD<, 8n, 1, 500<D, AspectRatio ® 1, AxesLabel ®8"Re@q1D", "Im@q1D"<D, 8j, 1, fmax<D

Plot@ImReInt@f D, 8f, h •. f ® 1, h •. f ® fmax<D

(24)

H*The two sections below is used to plot the system,

find the participation number and the total norm*L

P = Table@ Sum@Abs@gLines @@f, j, -1DDD^2 ,8j, 1, n + 2<D^2 • Sum@Abs@gLines @@f, j, -1DDD^4 ,8j, 1, n + 2<D ,8f, 1, fmax<D; norm = Table@ Sum@Abs@gLines @@f, j, -1DDD^2 ,8j, 1, n + 2<D, 8f, 1, fmax<D; IntP = ListInterpolation@P, 8h •. f ® 1, h •. f ® fmax<, InterpolationOrder ® 1D@fD; IntN = ListInterpolation@norm, 8h •. f ® 1, h •. f ® fmax<,

InterpolationOrder ® 1D@fD;

Print@"h=", h , " g= ", g, " w= ", w ,

" V= " , B, " n= ", n, " Res= ", fdif, " A=", AtD "P"

Plot@IntP, 8f, h •. f ® 1, h •. f ® fmax<,

PlotRange ®880, All<, 80, All<<, AxesLabel ® 8"h", "P"<D "Norm"

Plot@IntN, 8f, h •. f ® 1, h •. f ® fmax<,

PlotRange ®880, All<, 80, All<<, AxesLabel ® 8"h", "h"<D Sim1.nb 5

Intglines =

Table@ListInterpolation@Abs@gLines@@f, kDDD^2, 80, time<,

InterpolationOrder ® 0D@tD, 8f , 1, fmax<, 8k, 1, n + 2<D;

P =HSum@Abs@Intglines@@1, iDDD^2, 8i, 1, n + 2<DL^2 •

Sum@Abs@Intglines@@1, iDDD^4, 8i, 1, n + 2<D; tot = Sum@Abs@Intglines@@1, iDDD^2, 8i, 1, n + 2<D;

Dimensions@IntglinesD

Print@"h=", h, " g= ", g, " w= ", w ,

" V= " , B, " n= ", n, " Res= ", fdif, " A=", AtD

ParametricPlot3D@

Evaluate@Table@8i, t, Intglines@@1, iDD<, 8i, 1, n + 2<DD, 8t, 0, time<, PlotRange ® All,

AxesLabel ®8"n", "z", "ÈqÈ^2"<, PlotStyle ® Black, BoxRatios ®83, 5, 1<, PerformanceGoal ® "Quality"D

ParametricPlot3D@Evaluate@

Table@8hmin + fdif * i, t, Intglines@@i, 1DD<, 8i, 1, fmax<DD, 8t, 0, time<, PlotRange ® All, AxesLabel ® 8"h", "z", "ÈqÈ^2"<, PlotStyle ® Black, BoxRatios ®83, 5, 1<,

PerformanceGoal ® "Quality"D

ParametricPlot3D@Evaluate@

Table@8hmin + fdif * i, t, Intglines@@i, 1DD<, 8i, 1, 1<DD, 8t, 0, time<, PlotRange ® All, AxesLabel ® 8"h", "z", "ÈqÈ^2"<, PlotStyle ® Black, BoxRatios ®83, 5, 1<,

PerformanceGoal ® "Quality"D

Plot@tot, 8t, 0, time<, PlotRange ® 880, All<, 80, All<<, AxesLabel ®8"z", "N"<D

Plot@P, 8t, 0, time<, PlotRange ® 880, All<, 80, All<<, AxesLabel®8"z", "P"<D

6 Sim1.nb

H*This section needs to run in

order to calculate the lyapunov exponent*L

For@i = time * LyMin, i £ time - steglengdly, i += steglengdly,

Clear@KD;

K = Table@

uk@iD Š

KroneckerDelta@line, kD * Exp@I * RandomReal@80, 2 Pi<DD * stör + gLines@@f, k, Hi - LyMin * timeL • steglengd + 1DD ,8f, 1, fmax<, 8k, 1, n + 2<D;

Clear@linesD;

lines = Table@

NDSolve@8V@tD, bound1, bound2, K@@fDD<, U@tD, 8t, i, i + steglengdlyd <, MaxSteps ® 10 000,

InterpolationOrder ® 1, PrecisionGoal ® 12D

,8f, 1, fmax<D@@All, 1, line, 2DD ;

If@i Š time - steglengdly,

slist = Table@k, 8k, i, i + steglengdly, steglengdly<D, slist = Table@k, 8k, i, i + steglengdly, steglengdly<DD;

gLinesPart = lines•. t ® slist;

If@i Š time * LyMin, LygLines = gLinesPart, LygLines = Table@

Join@LygLines@@fDD, gLinesPart@@fDD, 1D, 8f, 1, fmax<DD;

D;

Print@"hela simuleringen tog ",

Sim1.nb 7

H*This section calculates the lyapunov exponent*L Print@"h=", h, " g= ", g, " w= ", w ,

" V= " , B, " n= ", n, " Res= ", fdif, " A=", AtD k = 0;

LyDif1 = Table@

LygLines@@f, 2 * j 1DD

-gLines@@f, line, HHj - 1L * steglengdlyL • steglengd + 1DD ,8f, 1, fmax<, 8j, 1, time * H1 - LyMinL • steglengdly - 1<D; k = 1;

LyDif2 = Table@ LygLines@@f, 2 * jDD

-gLines@@f, line, Hj * steglengdlyL • steglengd + 1DD

,8f, 1, fmax<, 8j, 1, time * H1 - LyMinL • steglengdly - 1<D; k = 2;

Lyaexp = 1• Htime * H1 - LyMinLL * Table@ Sum@

Log@Abs@LyDif2@@f, jDD • LyDif1@@f, jDDDD

,8j, 1, time H1 - LyMinL • steglengdly - 2<D, 8f, 1, fmax<D; k = 3;

IntLy = ListInterpolation@Lyaexp,

8h •. f ® 1, h •. f ® fmax<, InterpolationOrder ® 1D@fD; k = 4;

Plot@IntLy, 8f, h •. f ® 1, h •. f ® fmax<, PlotRange ® All, AxesLabel®8"h", "l"<D

Plot@IntLy, 8f, h •. f ® 1, h •. f ® fmax<, AxesLabel®8"h", "l"<D

k = 5;

H*The code below is used to plott

the manifold of the stationary case*L

ClearAll@"Global`*"D g = 0.1; h = 1; V = 1.9; w = 0.2; n = 100; M = 600; 8 Sim1.nb

(25)

dexp = - 1; Kint = 500; stegl = 10 ^ - 4; min = 0.00; max = 0.2; root = "D:\\Dropbox\\viktor\\"; SetDirectory@root <> "Saved\\"D; randomlist = <<Hroot <> "randomlist.txt"L;

slimrand = Table@randomlist@@jDD, 8j, 1, n<D;

G = -HV * slimrand@@1DD + w Abs@xD^2 L x -IHg • H1 + Abs@xD^2L - h * jL x + y;

Dimensions@FD;

ProgressIndicator@Dynamic@jD, 8min, max<D

Dynamic@kD

Dynamic@dexpD

Dynamic@ListPlot@A, Joined ® True,

AxesLabel ®8"x", "A"<, PlotRange ® AllDD

For@j = min, j £ max, j += stegl, a1 = 0; a2 = 10 ^ - 30; A1 = AbsoluteTime@D; For@k = 1, a1 < a2, k ++ , x = 2 * 10 ^ - 13 * k * 2 ^ dexp; y = 0 ; For@ i = n , i ³ 2, i += -1, t = x; x =HV * slimrand@@iDD + w - Abs@xD^2 L x + IHg • H1 + Abs@xD^2L - h * jL x - y; y = t; D; Sim1.nb 9 D; a1 = a2; a2 = Abs@GD; kt = k D;

If@kt > 2 * Kint, dexp = dexp + 1D; If@kt £ Kint, dexp = dexp - 1D;

A2 = AbsoluteTime@D;

If@a1 < 2 * 10^100,

z =8h * Hj - stegl • 2L, Ha1 + a2L • 2<; If@j == min, A =8z<, A = Join@A, 8z<DD; D D; Print@"h=", h, " g= ", g, " w= ", w , " V= " , V, " n= ", n, " stegl= ", steglD ListPlot@A, Joined ® True, AxesLabel ® 8"h", "A"<D

ListPlot@A, Joined ® True,

AxesLabel ®8"h", "A"<, PlotRange ® 80, All<D

EmitSound@Sound@SoundNote@DDD

H*The code below is used to plott the first turning

point of the manifold in the stationary case*L

ClearAll@"Global`*"D g = 0; h = 0.2; V = 1.9; w = - 1.00 ; n = 100; M = 500; d = 1; root = "D:\\Dropbox\\viktor\\"; 10 Sim1.nb

SetDirectory@root <> "Saved\\"D; randomlist = <<Hroot <> "randomlist.txt"L;

slimrand = Table@randomlist@@jDD, 8j, 1, n<D;

F =HV * slimrand + w - Abs@xiD^2 L xi+ IHg • H1 + Abs@xiD^2L - hL xi- xi+1; xn= Table@8 * 10^-15 * HjL, 8j, 1, M, d<D; xn+1= 0 ; Dimensions@FD ProgressIndicator@Dynamic@iD, 8n, 1<D For@ i = n , i ³ 2, i += -1, xi-1= F@@iDD D; a = -HV * slimrand@@1DD + w - Abs@x1D^2 L x1 -IHg • H1 + Abs@x1D^2L - hL x1+ x2; X = Table@8a@@jDD, x1@@jDD<, 8j, 1, HM - 1L • d + 1<D;

Z = Table@8a@@jDD, Sum@ Abs@xk@@jDDD^2, 8k, 1, n<D<,

8j, 1, HM - 1L • d + 1<D;

Y = Table@8Re@x1@@jDDD, Im@x1@@jDDD<, 8j, 1, HM - 1L • d + 1<D;

ListPlot@Abs@XD, Joined ® True,

PlotRange ®880, 1<, 80, 1<<, AxesLabel ® 8"ÈAÈ", "Èq1È"<D

ListPlot@Abs@XD, Joined ® True, PlotRange ® 880, 0.4<, 80, 0.2<<, AxesLabel ®8"ÈAÈ", "Èq1È"<D

ListPlot@Abs@x1D, Joined ® TrueD

ListPlot@Abs@ZD, Joined ® True, AxesLabel ® 8"a", "N"<D

Dimensions@XD

Dimensions@x1D

ListPlot@Y, Joined ® TrueD

U = Table@Dot@Abs@X@@iDDD - Abs@X@@i + 1DDD,

Abs@X@@i + 2DDD - Abs@X@@i + 1DDDD, 8i, 1, HM - 1L • d - 1<D;

ListPlot@U, Joined ® TrueD

H*ListPlot@

Sim1.nb 11

Abs@ZD

Characterization of attractors in a model for boundary-driven nonlinear optical waveguide arrays with disorder, gain and damping