Citation
Modeling the response of Rayleigh - Van der Pol oscillators to stochastic excitation near the Hopf bifurcation

Material Information

Title:
Modeling the response of Rayleigh - Van der Pol oscillators to stochastic excitation near the Hopf bifurcation
Creator:
Lu, Bingshen
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English

Thesis/Dissertation Information

Degree:
Master's ( Master of integrated sciences)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
College of Liberal Arts and Sciences, CU Denver
Degree Disciplines:
Integrated sciences
Committee Chair:
Asadi-Zeydabadi, Masoud
Committee Members:
Tagg, Randall
Simon, Burt

Notes

Abstract:
In this thesis we investigate the Rayleigh-van der Pol oscillator with additive noise. This dynamical system is a fundamental model that can be tested experimentally by the Wien-bridge electrical circuit. The analysis of the Rayleigh-van der Pol oscillator in the presence of noise provides a good understanding of the stochastic nonlinear dynamical system. The Rayleigh-van der Pol oscillator is a two-variable nonlinear system with one control parameter. As the control parameter varies at the bifurcation (transition) point the state of system changes from steady to self-sustained oscillation. This transition is known as Hopf bifurcation. The focus of this project is on the behavior of the system near bifurcation point in the presence of noise. We found that the bifurcation point widens to a “bifurcation region” and the ”bifurcation region” shifts as the level of noise increases.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
Copyright Bingshen Lu. Permission granted to University of Colorado Denver to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.

Downloads

This item has the following downloads:


Full Text
MODELING THE RESPONSE OF RAYLEIGH - VANDER POL OSCILLATORS
TO STOCHASTIC EXCITATION NEAR THE HOPF BIFURCATION
by
BINGSHEN LU
B.S., University of Colorado Denver, 2015
A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Integrated Sciences Integrated Sciences
2019


This thesis for the Master of Integrated Sciences degree by Bingshen Lu
has been approved for the Integrated Sciences Program by
Masoud Asadi-Zeydabadi, Chair Randall Tagg Burt Simon
Date: May 18, 2019
m


Lu, Bingshen (M.S., Integrated Sciences)
Modeling the Response of Rayleigh - Van der Pol Oscillators to Stochastic Excitation Near the Hopf Bifurcation
Thesis directed by Associate Professor Masoud Asadi-Zeydabadi
ABSTRACT
In this thesis we investigate the Rayleigh-van der Pol oscillator with additive noise. This dynamical system is a fundamental model that can be tested experimentally by the Wien-bridge electrical circuit. The analysis of the Rayleigh-van der Pol oscillator in the presence of noise provides a good understanding of the stochastic nonlinear dynamical system. The Rayleigh-van der Pol oscillator is a two-variable nonlinear system with one control parameter. As the control parameter varies at the bifurcation (transition) point the state of system changes from steady to self-sustained oscillation. This transition is known as Hopf bifurcation. The focus of this project is on the behavior of the system near bifurcation point in the presence of noise. We found that the bifurcation point widens to a “bifurcation region” and the ’’bifurcation region” shifts as the level of noise increases.
The form and content of this abstract are approved. I recommend its publication.
Approved: Masoud Asadi-Zeydabadi
IV


TABLE OF CONTENTS
I INTRODUCTION......................................................... 1
II ESSENTIAL IDEAS FROM NONLINEAR DYNAMICS.............................. 7
11.1 Linear Systems................................................... 7
11.11 Simple Harmonic Motion with Damping.............................. 9
11.111 Limit Cycle and Hopf Bifurcation................................ 11
II.IV A Nonlinear System: the Classic Van der Pol Model............... 11
II.V Modified Van der Pol Model...................................... 18
II. VI Rayleigh-van der Pol Model...................................... 21
III ESSENTIAL IDEAS FROM STOCHASTIC DIFFERENTIAL EQUATION
(SDE).............................................................. 27
111.1 Classical Example of Stochastic Process: Brownian Motion........ 27
111.11 Generating a Stochastic Process Numerically..................... 27
111.111 Embedding Stochastic Processes in Differential Equations....... 28
III. IV Numerical Methods for SDE...................................... 28
III.IV.I Euler-Maruyama method......................................... 29
III.V Solving Stochastic Integrals..................................... 29
III.V.I Ito’s Rule.................................................... 29
III.V.II Ito’s Chain Rule.............................................. 30
III. V.Ill Analytical examples......................................... 30
IV STOCHASTIC DYNAMICS................................................ 37
IV. I Stochastic Van der Pol Model.................................... 37
IV. II Stochastic Rayleigh-van der Pol................................. 37
V BEHAVIOR NEAR THE BIFURCATION...................................... 42
V. I Numerical Simulation............................................. 42
V.I.I Poincare Section.............................................. 42
v


V.I.II Behavior of Mean and Variance with Bifurcation Parameter for Fixed
Noise Level................................ 46
VI CONCLUSION AND FUTURE WORK............................. 60
VI.I Future Work.......................................... 60
VI.I.I Coupled System..................................... 60
REFERENCES................................................ 62
APPENDIX
A. ANALYSIS AND SIMULATION FOR COUPLED SYSTEM............. 64
vi


CHAPTER I
INTRODUCTION
Bifurcation, a sudden change of the dynamical system, lies in the center of nonlinear dynamics[29][19]. It has plenty of applications in daily activity and fields of science such as the quick change of medical condition from stable to severe, financial market change from balanced to crisis[ 19], or neural network change from stationary to generates motion. The needs to do a proper prediction of such phenomena rise as time goes by. However, reality has uncertainties so it is important to understand how randomness effects the nonlinear systems.
The general theory of random excitation in physical systems was first studied by Albert Einstein in 1905 and was generalized by Marian Smoluchowski in 1916[8]. Later on, physicists and mathematicians Andrey Kolmogorov, Adriaan Fokker, Max Planck, Leonard Ornstein, George Eugene Uhlenbeck, Norbert Wiener, Joseph L. Doob and Kiyoshi Ito laid the foundation for solving linear SDE between 1930s and 1950s[8], However, those earlier works are for solving linear SDE. Difficulties arise for nonlinear system because of the complexity nature of nonlinear systems[8].
Back in 1957, Garstens started studying noise in nonlinear oscillators[12], Garstens used a nonlinear Van der Pol equation with driving noise term to model nonlinear oscillators that ’’the output is fed back into the system thus modifying in a complicated manner the effective noise input”[12] here the nonlinear terms in such model are multiplicative (i.e. Xi X2, where ^ =X2). Garstens presented a method
1


for estimating the effect of nonlinear noise in an oscillator at low levels of oscillation. This required an approximate solution of a non-homogeneous van der Pol type equation with a noisy driven force by reducing the equation to a linear form. The power spectrum of the output was used to study the effect of noise in the nonlinear system.
Caugheys paper in 1959[6] showed that there was interest in the how nonlinear systems response to random excitation, especially in the systems with multiplicative nonlinear terms. He studied the self-excited system like Van der Pol Oscillator and its response to external stochastic terms acting as external random excitation.
Caughey found that the solution has a stochastic/noise term because of the stochastic term in the self-excitation system, and that the strength of the stochasticity decreases as the amplitude of oscillation increase. Interestingly enough, so does the band width of the stochasticity/noise which prove the early work done by Garstens incorrect.
In Caugheys paper in 1967[7], the use of a Fokker-Planck equation had been introduced to study the stationary response of self-excited oscillators. Using such technique, Caughey obtained the exact solution for the stationary joint probability density function, as well as the spectral density for oscillator with weak damping and weak excitation.
Following his previous work, Caughey developed a nonlinear theory of random vibrations[8]. As Caughey pointed out, the field’s importance had rise rapidly because of the developments in fields like high speed flight and the needs of dealing the random fluctuation cause by turbulence back at that time. Caughey summarized and developed some solution techniques for random vibrations in nonlinear systems[8].
During the similar time in Japan (1974), Kawakubo and Kabashima also used the Fokker-Planck approach to study stochastic processes in self-excited oscillation using a Van der Pol stochastic nonlinear system in a more experimental way[17], Kawakubo and Kabashima focused on how the fluctuations, namely the amplitude and phase
2


fluctuations in the system, behave in the experiment of Wien-Bridge oscillators.
They showed how the voltage fluctuates with the feedback factor and how the variance fluctuates with the parameter controlling the nonlinearity. Kawakubo and Kabashima also included the probability distribution function obtained with the Fokker-Planck equation, examining how it changed with the parameter controlling the nonlinearity.
In agreement with Caughey and others find. Spanos showed that in 1978 such problems were drawing more and more attention[28]. This included several more applications of stochastic phenomena in physics and engineering, such as control, aerospace, natural disaster, ocean wave and transportation engineering.
Spanos pointed out that even though linear random vibration has be well studied, the theory of non-linear random vibrations was not. This is still the case: there are still a lot of development needed in this topic.
The Fokker-Planck approach was used in Spanos work to obtain solutions, Spanos also showed the application of Fokker-Planck approach in Van der Pol and Rayleigh oscillators.
Later on in 1995, Leung did a systematic study of the behavior of Van der Pol oscillator with additive and multiplicative stochastic terms using different nonlinearity parameters[20], comparing the bifurcation behaviors of the system with different stochastic terms. By such methods Leung stated that we could get a more realistic understanding of the bifurcation of a nonlinear system.
Leung found that all level of stochastic terms / noise chosen disturb the usual behavior of the Van der Pol oscillator. Moreover, he found that noise could postpone the bifurcation between fixed point and limit cycle and that noise could induce a transition from a stable fixed point to a stable limit cycle, and vice versa. Leung suggested that noisy nonlinear systems could have more and different phenomena than the deterministic one, such as stochastic resonance.
3


In Leung’s later work in 1998[21], he pointed out studies on stochastic bifurcation focus on the shift of bifurcation point, where a shifting of bifurcation point is an indication of the noise effect on the two competing mechanisms, namely self-sustained oscillation and nonlinear dissipation.
Yuan and others made further developments in this direction in 2013 and found that zero-mean multiplicative noise can induce the bifurcation points and sometimes even change the topology of the system[31], or in another words, the structure does not correspond to the deterministic counterpart, which make the analysis of such system more difficult. They also used the steady-state distribution of the process to approach such case, in where correspond means the probability density function does not change over time and have a maximum at the stable attractors, which later on define as the steady state distribution multiplying negative one being a Lyapunov function of the deterministic counterpart of the stochastic dynamics. Yuans group apply a new stochastic interpretation rather than the common Ito and Stratonovichs integration to obtain a correspondence between stochastic and deterministic dynamics.
Just last year (2018), Doan and others were trying to merge bifurcation theory for stochastic dynamical systems[10], Doan and others work shows that the earlier concepts developed in the 1990s, namely theory of phenomenological (P) bifurcation and dynamical (D) bifurcation, cannot comprehensively capture the complexity of bifurcation in random dynamical system. Doan points out that despite its relevance for many applications of topical interest, a bifurcation theory of random dynamical systems is still in its infancy[10],
With the developments made in mid 20th century on time dependent probability theory[2][16][26] and stochastic differential equations[25], and we began to have tools to solve the stochastic system in the form of stochastic differential equation[13][25]. Yet The bifurcation of stochastic chaotic system is still something we do not sufficiently understand yet. But the importance of bifurcation in such system is becoming
4


increasingly unavoidable[l][10].
In recent years, neuroscience (as well as machine learning and artificial intelligent) draws more and more attention[22] and funding, so do the nonlinear and stochastic components of this field[l][22][15]. There are some breakthroughs on the applications of such topics in crucial areas such as neuroscience[14][l 1][17], However, most of the developments about stochastic differential equations are in the field of stock market forecasting and related topics in quantitative/financial engineering[25], i.e. the early detection of economic bubbles[9]. How to analyze the stochastic nonlinear dynamical model is still a difficult task which been rarely discussed in the context of stochastic differential equation but mostly done by deterministic approximation, especially for nonlinear system with multiplicative terms in neural networks (i.e. System contains Xi X2, where ^ = X2). Most of the analytical and numerical methods cannot be applied to such terms[5][18][25]. Yet the interest in self-excited oscillation and its application in neural network pushes us forward to study such system[3][17], especially the response of a class of oscillators to stochastic excitation[28]. Now we might have tools and methods to study it in the content of stochastic differential equation instead of deterministic approximation methods.
In this thesis we are building a modified Rayleigh Van der Pol model with additive stochastic components. We hope to test if the system behaves as others observed, when we put the system in the proper context of stochastic differential equation instead of deterministic methods. We also hope to lay a foundation for future work and potential applications of the Rayleigh Van der Pol model. The applications include neural networks for the Rayleigh Van der Pol[6][7] based central pattern generator (CPG) model[3], as well as other potential applications in physics, neuroscience, finance and other related fields.
This work is also offered as a guide for some methods of how to simulate and analyze models of nonlinear system with stochastic components in the context of
5


stochastic differential equations. This is primarily focused on how to solve the system with terms containing random variable products Xi • X2, where = X2, which we
do not have many good options to solve yet.
6


CHAPTER II
ESSENTIAL IDEAS FROM NONLINEAR
DYNAMICS
This chapter organizes a conceptual framework for analyzing two-dimensional nonlinear dynamics. This provides a broad and powerful foundation for analyzing many applications especially after stochastic forcing is added in later chapters.
II.I Linear Systems
In general a two-dimensional dynamical linear system has the form of
X-\ — dy\X-\ + d-\2X2 (II. la)
X2 = 921-^1 + 922-^2 (II.lb)
where x 1 = —, x 2 = —, and an, 9i2, 9i3, 3u are constant parameters.
dt dt
We can write the system in matrix form as
where A
x = Ax
[?n ai2D ^1 â–¡
and x =
â–¡ â–¡
a2i a22 â–¡ â–¡
x2
We seek solution of the form [29],
(II 2)
7


x=eAv. (II. 3)
Substituting (2.3) into (2.2), we get
Av = Av,
(II 4)
We say the solution in (2.3) is an eigensolution if v is an eigenvector of A with corresponding eigenvalue A [29],
The characteristic equation of matrix A is given by
Equation (2.5) yields
di 1 ~A a-i2 det
C&21 ^22 —A
o,
(II. 5)
where
A2 -tA + A = 0
(II. 6)
t = trace(A) = a-n+ 322 A = det(A) = 311 322 -312321
(II. 7a) (II. 7b)
Solving equation (2.6), we get
M =
A2 —
V
y+ 7-2 _4A
V 2
T - T2 - 4A
2
(II. 8 a) (II. 8b)
where A1 and ^2 are the eigenvalues of matrix A. We find the eigenvectors by substituting eigenvalues into one of the equations in (2.4). We get two linearly independent eigenvectors Vi and V2. Now we can get the general solution for x,
x=c-\eMtv-\ +C2eA2tv2.
(II 9)
8


II.II Simple Harmonic Motion with Damping
One good example for the linear system is simple harmonic motion with damping. The F in this case is the sum of the spring and damping forces:
dx
F=-kx-b~ (II. 10)
where k and b are the spring and damping constants. Applying Newton’s second law for (2.10), we get
(fix dx
m_____= -kx - b—
dt2 dt
We can rewrite equation (2.11) as
or
where (3=^, U)0 =
(fix dx
m_____+ b____+ kx = 0
dt2 dt
x + /3x + (fiQx = 0
The solution (2.12) and (2.13) has the form
x = Ae*(
Substituting (2.14) into (2.12), then we have
(11.11)
(11.12)
(11.13)
(11.14)
(,mA2 + bA + k)AeAt = 0
Equation (2.15) for the nontrivial solution of (2.12) implies
>2 b i k
A2 +___A +____= o.
m m
Equation (2.16) is the characteristic equation of (2.12).
(11.15)
(11.16)
9


The roots of equation (2.16), eigenvalues of (2.12), for underdamped or oscillatory
case, (|L < J-), is given by
b
k b2
A = — ±i - 9
O m m A
(11.17)
The general solution of (2.12) or (2.13) for weak damping, 4hs< is given by
b j K x = Aie^te m - b2 f A - b~ f Am- -h A2e2nl m Am~ (11.18)
or V V
x = A e^e' -fPXtt a gjSfg-/ w2-p2/4t (11.19)
1 2
For real valued solutions, A-\ and A2 will be complex conjugates.
We call x = 0, x = 0 a fixed point solution. For (3 < 0 all solutions tend to this fixed point solution and it is therefore said to be stable, for /3 > 0 the fixed point is unstable.
Figure (2.1) shows time series and phase plot for damped harmonic motion.
Figure II. 1: Underdamped Harmonic Motion
10


II.III Limit Cycle and Hopf Bifurcation
A limit cycle occurs when the phase plane (x v.s. contains an isolated closed trajectory. Isolated means the neighboring trajectories either spiral towards or away from the limit cycle. [29]
A Hopf bifurcation happens where a fixed point of a system changes to a limit cycle.
II.IV A Nonlinear System: the Classic Van der Pol Model
The Van der Pol equation has been used as a model for neuron networks with feedback controls. It is named after Dutch physicist Balthasar van der Pol[29][30], who applied the equation to electronic circuit oscillators.
CpX 2 Cfx 2
dt2 ~(P )Wo dt + W0X = (11.20)
Provided that [J /= 0 we can write:
CpX ( 2 cfx 2
dt2 ~^1 ~fjx )'w° dt + = 0 (II.21)
Using the rescaled variables, t =U)ot and x = ^ equation (2.21) becomes:
cfx
dt2 ~/J( 1 -x )
dx
2 —
dt
+ x = 0
(11.22)
This is the standard, or classical form of the van der Pol equation. However, the scale factor ^becomes singular at jt/ = 0 at the same point that the actual output, x, has an amplitude that is just about to rise from a value of zero[30]. This is one of the reasons that we modified the classical Van der Pol model as shown in the later section.(see section 2.3 and 2.4)
An alternate dimensionless form of the dynamics is the single scaled nonlinear
11


equation in x:
12


2
dx
cPx
dp-^-x)d;+x=0'
(11.23)
Here, the scaled time and dynamic variable are t = U)ot and x = - where
x* =V -ipor notational clarity, we will drop the tilde’s in the following but assume the equations are still in the scaled form.
We can applying method of averaging to the Van der Pol equation to find approximate solutions to the nonlinear system.
The dynamics, equation (2.23), can be re-written as two first-order differential equations
dx = w dt dw 2 (II. 24a)
dt=~x + Qj -x )w. (II. 24b)
We are going to assume that the the solution will evolve from the simple sinsoidal solution when the nonlinear term is turned on. Thus we think of time variations in amplitude and phase around the limit cycle such that the phase-space trajectory can be written in the form:
x(0 = /\(0cos[f+ w(t) = -A(t) sin[f+ Our goal is to find equations for the amplitude function A(t) and the phase function In order to do so, we insert (2.25a) and (2.25b) for x and w into the dynamical equations (2.24a) and (2.24b) to obtain dA(t) dw(t}
df cos[t+cp(t)]~A(t) 1+ df sin[t+cp(t)] = -A(t)sin[t+cp(t)] (11.26)
_ cW(f)sin[f+cp(f)] ~A(t) 1+ cos[t+cp(f)] =
dt dt
1
-A(f)cos[t+(p(t)]-fJ-A2(t)cos2[t+cp(t)] /\(f)sin[f+ 13


This simplifies to
dA(t) dt
cos [t+cp(t)]-A(t)
d(p(t)
dt
sin[f + (p(t)] = 0
(11.28)
dP^l. sin[t+ cp(t)] - A(tcos[t+ cp(t)] = dt dt
1
- fj ~A (£) cos [t + (p(t)] A(t) sin[f + (p(t)]. (11.29)
Use (2.28) to write
dcp(t) = ^cos[f+ cp(t)]
dt A(t) sin[f+ cp(t)]
Substitute (2.30) into (2.29) and multiply through by sin[f + (pit)]:
(11.30)
dA{t)
dA{t)
dt
sin [t + (p(t)]
dt
cos [t + cp(t)] =
1
-E-A (f)cos [?+cp(t)] A(t)sin [?+ Collect terms and use the fact that cos2[f + + sin2[f + dm 22 1 2
dt = pi-A (t) cos [t + (Pit)] A(t) sin [t+(p(t)]. (11.32)
We can also get
d(P(f) o o 1
dt
= fj -A (t)cos [t+(p(f)] cos[f+^(f)]sin[f+^o(f)].
(11.33)
So far this has been an exact change of coordinates. We now make an approxi-
mation. Anticipating that A2(f) ~ 0(p), we take both time derivatives
mu
dt
and
dW)
dt
to likewise be 0(p) and thus the amplitude and phase are slowly varying deviations
from the circular limit cycle. We will then treat them as effectively constant over one
cycle and replace the right-hand side with averages. Taking G =t + (p(tX dA(t) i 2 n
0 fj-A^t) cos P A(t) sin £)dd (II.34a)
dt 2 TT
1
277
1
277
14


cos 0sin ddd
(II. 3 4b)
VA(t)
2 TT
sin ddd ~A (t)
2 TT
2fjA(t) - g/\ (t).
(II.34c)
15


So we have
dA(t) 1 2
dt - 8 4a/ -/\ (0 /\(0- (II.35)
This is an example of an *amplitude equation*. It maybe formally integrated starting with
to arrive at
At long times
8-
dA
[4/7 -A2] A
= dt
(11.36)
V-
V
A+ 2 fj
V = dJt.
A2{2 jt/ -A)
(11.37)
V—
A —>2 y/j
(11.38)
(justifying thies assumption that A2 ~ 0(jL/)), while the next order correction is
V -L -j
- /j 2 -e ^
(11.39)
Returning to dimensioned coordinates,
( ) i/
x~ 2-e pWof x* fj cos(aJot), and the predicted final steady oscillation is given by
(11.40)
X —>2X* jt/ cos(coof).
where we recall that x* = ^
(11.41)
For future purpose, we re-write the classcial Van der Pol system as two first-order equations,
dx
— =y (11.42)
dt
16


dx
dt
2
(11.43)
= \j <1 -x ) y-x
As we mentioned the Hopf bifurcation happens where a fixed point changes to a limited cycle. The fixed point of classic Van der Pol model is (.x*, y*) = (0, 0). In this system, the Hopf bifurcation happens at jL/ = 0, when the system from steady state figure (2.3) turn into a self-sustain system (2.2).
Figure II. 2: Classical Van der Pol model, Stable Fixed Point, X versus t for /j = -0.1.
17


Figure II. 3: Classical Van der Pol model, Stable Fixed Point, V versus t for /j = -0.1.
Figure II. 4: Classical Van der Pol model, Stable Fixed Point, Phase Portrait for /j = -0.1.
We can see from both the figures of X and Y versus t that the model started oscillating and damping to zero after some time. The phase portrait shows the red dot, which is our initial point, spiral in and become a fixed point, that indicates the model reaches a fixed point.
18


Figure II. 5: Classical Van der Pol model, Limit Cycle, X versus t for /j = 0.1.
0 50 100 150 200 250 300
t
Figure II. 6: Classical Van der Pol model, Limit Cycle, V versus t for /j = 0.1.
19


> 0 •
-2
-1
2 â– 
1 â– 
-3 -2
-1
0
X
2
3
Figure IF 7: Classical Van der Pol model, Limit Cycle, Phase Portrait for/j = 0.1.
We can see from both the figures of X and Y versus t that the model oscillates at an fixed amplitude after some time. The phase portrait shows the red dot, which is our initial point, spiral out and become a circular trajectory, that also indicate the model reaches a limit cycle.
II.V Modified Van der Pol Model
The classical Van der Pol is a good simple model for the current goal. However, during the numerical simulation, the simulation breaks when [J is close to zero, which agree with the discussion in section 2.2. Therefore, we are using one version of modified Van der Pol model,
Comparing the second equation in classical with modified Van der Pol equations, (comparing equations (2.43) with and (2.45)) we find [J is outside the parenthesis in (2.43) but is inside in (2.45).That makes the nonlinear term more smooth (with no singularity) and more close to real world phenomena. For this reason it is easier to
dx
(11.44)
dy
2
(11.45)
20


observe the model’s behavior near bifurcation. The behavior of the model is more
consistent and does not break during simulation when jt/ is close to zero. For modified Yan der Pol Model, the bifurcation still happens at fj = 0, as figures 2.8-2.13 indicate.
Figure II. 8: Modified Van der Pol model, Stable Fixed Point, X versus t for /j = -0.1.
Figure II. 9: Modified Van der Pol model, Stable Fixed Point, V versus t for /j = -0.1.
21


0 0075 •
0.0050 ■ 00025 • 0.0000 ■ -0 0025 • -0 0050 ■ -0 0075 • -0 0100
-0010 -0005
0000
X
0005
0010
0015
Figure II. 10: Modified Van der Pol model, Stable Fixed Point, Phase Portrait for /j
= -0.1.
Given the same initial condition and parameters, we can see that the modified Van der Pol behaves almost the same with the classical counterpart. However, the modified Van der Pol model reaches fixed point faster than the classical model because the nonlinear term [J is less destructive.
Figure II. 11: Modified Van der Pol model, Limit Cycle, X versus t for [j = 0.1.
22


t
Figure II. 12: Modified Van der Pol model, Limit Cycle, V versus t for /j = 0.1.
06
04
02
>- 00 -0 2 -0 4 -0.6
Figure II. 13: Modified Van der Pol model, Limit Cycle, Phase Portrait for/j = 0.1.
Given the same initial condition and parameters, we can see as well that the modified Van der Pol behaves almost the same with the classical counterpart. However, the modified Van der Pol model reaches limit cycle faster than the classical model because the nonlinear term [j is less destructive.
II.VI Rayleigh-van der Pol Model
We use (Modified) Rayleigh-van der Pol equations as our standard system. One of the characteristics worth noting about Rayleigh-van der Pol equations is that its
23


phase portrait is circular, even for larger value of jL/ which make our future work on the phase plane a lot easier.
___= y (11.46)
dt
dy_ 22
dt = (ju-y -x) y-x (11.47)
Rayleigh-van der Pol is also a good representation of the Wein bridge circuit which can be conduct into experiment for testing [30],
Now we change the system to polar coordinate,
x = rcos 6 y = rsin 6
(II.48a) (II.48b)
X________
r= x2 + y2 (II. 49a)
0 = arctan — (II. 49b)
x
24


1
dr
dx
= ✓ (2x +2v ) dt 2 X2 + V2 dt dt (II. 50a)
1 dx dy
Is + >y + CNI X II (II. 50b)
1 2 2
= xy + y[Qj -y -x )y-x] (II. 50c)
1 = r(xy+ijy2 -y4 -xV -xy) (II.50d)
CM X 1 % 1 3 -H| >*. II (II.50e)
= - X2 -y2) (II.50f)
l2 sin2 0 9 [A/-(X2+y2)] r (II.50g)
= (jj-r frsin (II. 5 Oh)
dO 1 c/y c/x — = —(x— - y—) dt rf dt dt J. 2 2 2 (11.51 a)
= ?xt0u~y -x )y -x] -y _I 2 2 2 2 (II. 5 lb)
= r2 [W ~y -x )xy -(x + y )] (II. 51c)
= 2 [(/J -r )r sin 0 cos 0 -r ] (II.5 Id)
= (jt/ -r f sin 0 cos 0-1 (II.51e)
The constant-amplitude solution is
II (II. 52a)
0 = -1. (II.52b)
It predicts a pure sinusoidal oscillation at all values of jt/, that is, one that contains no harmonics - unlike the van der Pol model. Moreover the limit cycle is a circle of
V _
fixed radius proportional to jt/ in (x, y) phase space. Figures figures 2.6 and 2.7
25


show the bifurcation at jL/ = 0.
26


Figure II. 14: Rayleigh-van der Pol model, Stable Fixed Point, X versus t for /j 0.1.
0.0075 00050 00025 0.0000 -0 0025 -0 0050 -0 0075 -0 0100
0 50 100 150 200 250 300
t
Y
fr
1------------r------------i-----------i------------r
Figure II. 15: Rayleigh-van der Pol model, Stable Fixed Point, Y versus t for g 0.1.
27


00075
0 0050 ■ 00025 • 0.0000 ■ -0 0025 ■ -0 0050 ■ -0 0075 • -0 0100
-0010 -0005
0000
X
0005
0010
0015
Figure 11.16: Rayleigh-van der Pol model, Stable Fixed Point, Phase Portrait for/j =
-0.1.
Given the same initial condition and parameters, the figures shows that the Rayleigh-van der Pol model reaches the fixed point as fast as the modified Van der Pol model with a different oscillation trajectory and a more circular phase portrait even for the case of fixed point.
Figure 11.17: Rayleigh-van der Pol model, Limit Cycle, X versus t for g = 0.1.
28


03
02 0.1 00 -0 1 -0 2 -0 3
t
Figure II. 18: Rayleigh-van der Pol model, Limit Cycle, Y versus t for g = 0.1.
Figure 11.19: Rayleigh-van der Pol model, Limit Cycle, Phase Portrait for g = 0.1.
Given the same initial condition and parameters, the figures shows that the Rayleigh-van der Pol model reaches the limit cycle as fast as the modified Van der Pol model with a different oscillation trajectory and a more circular phase portrait. The black circle in the phase portrait indicate the constant-amplitude solution where
v_
r = ij.
29


CHAPTER III
ESSENTIAL IDEAS FROM STOCHASTIC DIFFERENTIAL EQUATION (SDE)
III.I Classical Example of Stochastic Process: Brownian Motion
A Brownian motion Bt (Wiener process Wi), is a sequence of random variables satisfies the following conditions [25]:
1) Bt(0) = 0 with probability 1,
2) for 0 V_____
distributed with mean 0 and variance t~S, or Bt~Bs ~ t ~sN(0, 1), where N(0, 1) is normal distribution with zero mean and unit variance.
3) in 0 III.II Generating a Stochastic Process Numerically
For the purpose of the numerical simulation, we discretize Brownian Motion as
V-
dWt ~ tN (0, 1). A Sample path of Brownian motion is shown in figure 3.1, it is a random path start at 0 with mean 0 and variance dt.
30


Figure III. 1: sample path of Brownian Motion
With the help of measure theory, it has been proved that Brownian Motion like showed in figure 3.1 is nowhere differentiable [2], we can consider it as a special case of fractal to understood its non-differentiablity.
III.III Embedding Stochastic Processes in Differential Equations
A stochastic process [25] has a general form of
dXt = f{Xt, t) dt+ g(Xt, t) dBt (III. 1)
where Bt stand for Brownian motion/Wiener process(Wf ), f(Xt, t)dt is the drift term which carry the deterministic characteristics and g(Xt, t)Bt is the diffusion term which carry the stochastic characteristics.
It is worthy noted that because Brownian motion is nowhere differentiable [2], which lead to a result that stochastic process is nowhere differentiable as well, therefore we need write it in difference equation form before we try to solve it analytically and numerically.
III.IV Numerical Methods for SDE III.IV.I Euler—Maruyama method
31


A stochastic differential equation in general in the integral form is given by
t t
X(t) =X0 + f(X(s))ds+ g(X(s))dW(s), 0 0 0
In general we use ltd version of integral since most of the methods in stochastic differential equation are based on ltd's chain rule numerically and analytically.
Let’s write this integral solution in differential form:
dX(t) = f(X(f))dt + g(X(t))dW(f), X(0) = X0, 0 To solve this equation numerically, we let At = T/N where N is the total steps andXj=X{tj). The Euler-Maruyama (EM) method is similar to Euler’s method for ODE with additional stochastic term.
Xj — Xj-1 + ftXj-^At + g(X/-i)Al/l/, / — 1, 2...A/.
For the case of additive noise, we write the diffusion terms
V____
AW ~A AtN (0, 1) (III.2a)
= AN(0, At) (III.2b)
where A is a constant parameter that influence the strength of noise. The level of stochasticity is controlled by this parameter.
III.V Solving Stochastic Integrals III.V.I Ito’s Rule
Ito’s rule prove that for stochastic process dXt = f{Xt, t)dt+ g(Xt, t)dBt,
df=0,dtdBt=0 (III. 3 a)
dBfdt = 0, dBf = dt (III. 3b)
This is the one of the major breakthrough leading the method to solve the stochastic integrals.
32


III.V.II Ito’s Chain Rule
III.V.II.I Single Variable Ito’s Chain Rule
Ito’s Chain Rule says, for Yt = g(t, Xt), where Xt is stochastic process, we have
dg
dg
dYt = Xt)dt + (t, Xt)dXt + -—
dt dx 2 dx2
III.V.II.II Multivariate Ito’s Chain Rule
i d2g 2
(t, Xt)(dXt)
(HI 4)
for Yj= g[t,Xi), where Xt is stochastic process, we have
2 2 2 2
da. da. I d_g_
dY = Ti(t,Xt)dt+ f*(t,xt)dx,+ L (tXDdXjdx, (III 5)
dt 1=1 dxi 21=1 Sx
1=1 '
for i = 1,2. Note that aa 9' = 0 Vj, k.
7 dxjdXk J
i= 1 k=*\
III.V.III Analytical examples III.V.III.I Simple Brownian Motion
We can see the main difference between Stochastic differential equation (SDE) and ordinary differential equation (ODE) by integrating over Brownian Motion B{.
t
Bs dBs.
(HI 6)
first we guess the form of the solution from deterministic model, in this case, it
is y = g(t, x) = yX2, then we try the form for Yt.
Y=g(t,X)=l_B2
i i 2 f
Using Ito’s chain rule we have
(III. 7)
1 2
dYt=0 dt+B t dBt+ â–  1 (dB)t
(III. 8)
d(Is2) = BdB + l_dt
^ f where (dBt)2 = dt by Ito’s Rule,
30
(III. 9)


t t
2
31


1 o2 _ r 1 *
-°t~ Bs dBs + - l ds. (III. 10)
2 Q 2 o
Which give us the solution,
BsdBs= l-B2(-lt (HI 11)
0 z z
III.V.III.II RLC circuit
For a simple RLC circuit, we have system
1
LQt + RQt + £ Qt = Ft (III. 12)
where Ft = Gt+aWt,
Wt : White noise (Wtdt = dBt),
Bt : one dimensional Brownian Motion,
Qt : the change at time t at a fixed point in an electric circuit,
Ft : the potential source at time t,
Gt: the potential source at time t without noise, R : the inductance, L: the resistance, C: the resistance.
We rename the variable as following:
□*1 □ [£h
X = x,= c C= C II (III 13)
x2 Qt
X2 Qt
We write equation (3.12) as a system of the first order SED:
X1 =X2 (III. 14a)
1
Z_X2 = -RX2 - Xi+ Gt+ aWt (III. 14b)
o
Now we write the system of SED, equations (3.14a) and (.14b), in the from of difference equation as following
32


1
1
a
(III. 15 a)
dXi =X2dt R
dX2 = ~LX2dt - Xi dt+ LGtdt+ Wtdt
(III.15b)
where Wtdt= dBt.
In order to solve this two dimensional SDE, we represent the system, (3.15a) and (3.15b), in the matrix form equation,
dX
dX
dXt = AXdt + Hdt + KdBf
â–¡
0 1

CL L
â–¡ â–¡ â–¡ c
dX2 -i—fi-
ll
â–¡
â–¡ â–¡ 0 L
â–¡ â–¡
E.
(III.16a) (III. 16b)
Now solving the matrix equation (Note that e At = °^=Q AKkfK converges.), we
move AXtdt to the left side of equation (3.16a) and multiply both sides by e~At.
dXt -AXtdt = H(t)dt + KdBt (III. 17a)
e~At[dX-AXdt] = e~At[H(t)dt+ KdBt] (III. 17b)
Now we want to show if we can write the left hand side of equation 3.14b as d(e~AtX). In attempt to figure this out, we try to apply multivariate Ito’s chain rule
x2
to a(t xi, xo) = e At x'
dY = d(e~AtX) = -Ae~AtXdt+ e~AtdX (III. 18a)
= e~At[dX-AXdt].
(III.18b)
(Note that Ae At= e AtA.)
33


d(e~AtX) = e~At[H(s)dt+ KdBt\.
(III. 19)
Integrate both sides,
t t
e~AtX -Xo = e~AsH(s)ds+ e~AsKdBs o o
(III. 20)
Recall integration by parts we have
t t
f(s)dBs = f(f)Bt - ds(f(s))Bs (III.21)
0 0
Use integration by parts for the last term on the right hand side in equation (3.20),
t t
e~AsKdBs= e~AtKBt - (-A)e~AsKBsds
o o
(III.22)
Substituting (3.22) into (3.20), we get
t t
e~AtXt = X0+ e~AsH(s)ds + e~AtKBt + Ae~AsKBsds o o
(III. 23)
Which give us the solution
Xt = eAtX0 + KBt + eAT e~As[H(s)+AKBs]ds (III.24)
*At „-As
III.V.IILIII Damped Harmonic Oscillator
The equation of motion for Damped Harmonic Oscillator with Brownian Noise is given by
34


dx Bt d2Xt
-/cx -6 “+ a “ = m ~ (HI 25)
Bt: Brownian noise (Wtdt= dBt where Wt is white noise), k : spring constant, b : damping constant, a : constant, m : mass.
Now let us to write it in the content of stochastic differential equation.
d2X dXt ,v_ dB< 171 dt2 +b dt+ f a dt (III.26)
Equation (3.26) can be written as a system of two first order SDE as following
X, =X2 (III.27a)
mX2 + bX2 + AXi = aBt (III.27b)
Writing the system of SDE in the from of difference equation, we have
dXi =X2dt (III.28a)
b k a
dX2 = - ~X2dt -~X1 dt+ ~dBt m m m (III.28b)
Now we rewrite the system in the form of matrix equation,
â–¡ II G' >< ii X (III. 29)
x2
x2
35


dX = AXdt+ CdBt
(III.30a)
:dX] â–¡
dX =ur u, A
â–¡ C dX2

(III. 3 Ob)
c =
Moving AXdt to the left side of (3.30a) and multiplying both sides by e At, we
get
dX -AXdt = CdBt (III. 31 a)
e~At[dX-AXdt] = e~AtCdBt (III. 3 lb)
We already know from the previous example that we can write the left hand side of equation 3.31b as d(e~AtX).
Which give us
d(e~AtX) = e~AtCdBt (III.32)
Integrate both sides,
t
e~AtXt -Xq = e~AsCdBs (III.33)
o
Use integrate by parts for the last term in the right hand side in equation 3.33,
t t
e~AsCdBs = e~AtCBt - (-A)e-AsCBsds (III.34)
o o
Bring the term in we have
t
e~AtXt = X0+e~AtCBt Ae~AsCBsds (III.35)
o
Which give us the solution
36


t
Xt = eAtXo + CBt + eAt Ae~AsCBsds (III. 3 6)
o
It is worth noting that the damped harmonic oscillator and RLC circuit in the previous section have similar form of solution. In the RLC circuit, in addition to the noise, the external source is also included in the system. The stochastic RLC in previous section mathematically is similar to the stochastic damped driven harmonic oscillator. Figure 3.2 show the result for the case of damped harmonic motion with weak noise (A = 0.01) for weak damping < ^), we can see that for small A, the stochastic term did not influence the damping that such, but started to take over after. The phase portrait indicate that the stochastic term keep kicking the system out of the fix point.
Figure III.2: Damped Harmonic Motion with A = 0.01, m = 5, b = 0.01, k= 1.
37


CHAPTER IV
STOCHASTIC DYNAMICS
IV.I Stochastic Van der Pol Model
In general, the system always have some level of noise. We take the external additive white noise into consideration and treated it as a forcing term.
The stochastic Van der Pol model with additive noise has the form
dX=Y dt + dW^t) (IV.la)
dY =(jj -X f Y dt -X dt + dW2(f) (IV. lb)
IV.II Stochastic Rayleigh-van der Pol
Stochastic Rayleigh-van der Pol equation is given by
dX = Y â–  dt + dW^f) (IV. 2a)
dY =Qj-Y 2-X )2V dt -X dt + dW2(f) (IV.2b)
We found out that the analytical method for linear stochastic differential equations does not work for multiplicative variables. We cannot rewrite stochastic nonlinear system into matrix equations when there are multiplicative variables. Moreover, we found no theory about analytical method for stochastic nonlinear systems with
38


multiplicative variables. However, we can solve such system using Euler-Maruyama method.
In the case of circuit, we set it up in a more realistic way such that dW-\(f) = 0.
dX= Y dt (IV. 3 a)
dY =Qj-Y 2-X )2 Y dt -X dt + dW(t) (IV.3b)
In figures 4.1 to 4.6 the time sires and phase portraits are shown for fj = ~0.1 and fj = 0.1, below and above deterministic bifurcation value fj = 0, for A = 0.01.
Figure IV. 1: Rayleigh-van der Pol model, Stable Fixed Point, X versus t for jj = -0.1 and A = 0.01.
39


300
Figure IV.2: Rayleigh-van der Pol model, Stable Fixed Point, Y versus t for jj = -0.1 and A = 0.01.
Figure IV.3: Rayleigh-van der Pol model, Stable Fixed Point, Phase Portrait for/j = -0.1 and A = 0.01.
Figure 4.1 to 4.3 show that in the case of Stable Fixed Point, the figures for X and Y versus t suggest that the oscillation damps away and the stochasticity takes over. The phase portrait is now a fuzzy circular trajectory: the trajectory fluctuate while reaching the stable fixed point. Right after the system reach the stable fixed point, the stochasticity/noise kick the system out of it and that cause the system to restart damping.
40


Figure IV.4: Stochastic Rayleigh-van derPolmodel, Limit Cycle, X versus tfor/j = 0.1 and A =0.01.
04
0.3
02
0.1
>- 0.0 -0.1 -0 2 -0 3
Figure IV. 5: Stochastic Rayleigh-van derPol model, Limit Cycle, Y versus tfor[j = 0.1 and A =0.01.
41


04
0.3 â– 
02 -01 •
> 00 â– 
-0.1 â– 
-0.2 •
-0.3 •
-0.6 -04 -0.2 0.0 0.2 04 0.6
X
Figure IV.6: Stochastic Rayleigh-van der Pol model, Limit Cycle, Phase Portrait for IJ = 0.1 and A = 0.01.
Figure 4.4 to 4.6 showed that in the case of limited cycle, the figures for X and Y versus t suggest that the model have a fuzzy limit cycle instead, the fuzzy trajectory and widened cycling region of the phase portrait agree with it. The black circle in
_ V
the phase portrait indicate the constant-amplitude solution where r = jJ
42


CHAPTER V
BEHAVIOR NEAR THE BIFURCATION
V.I Numerical Simulation
In this chapter the numerical Euler-Maruyama method is used to solved the stochastic Rayleigh-van der Pol equation.
V.I.I Poincare Section
Study and statistical analysis of a stochastic system provides useful tools for understanding the behavior and properties of the system. In order to obtain and measure the statistical behavior of the nonlinear stochastic system, we chose to measure the statistical properties of the Poincare section.
Numerical simulation of Poincare section for deterministic system requires Jacob-nian which need take the (partial) derivatives of the system. That bring the difficulty since stochastic differential equation is nowhere differentiable.
In order to obtain the Poincare section for our system, we set up a (partial) hyper plane orthogonal to the phase plane which also goes through the phase plane’s x axis. Every time the trajectory in phase plane cross the hyper plane, we mark it. This can give a good estimate of system’s statistical properties. We only need to collect data for any one of the axis(i.e. the positive x axis) that the trajectory cross thanks to the circular phase portrait of Rayleigh-van der Pol.
43


0 4
0 3 02
01
>â–  00 -0.1 -0 2 -0 3
-0.6 -04 -0.2 0.0 02 04 06
X
Figure V. 1: Poincare Section for Phase Portrait of Stochastic Rayleigh-van der Pol model under different noise level, /j = 0.1, A = 0.01, t > 1300, fort e[0, 2000].
As figure 5.1 shows, under the influence of relatively weak noise (A = 0.01), the system starting became slightly unstable. The stochastic term causes the fluctuation of the limit cycle, which showed as the fuzzy path and the widened circular region of the phase portrait. As the strength (the level) of noise gets larger the fluctuation also increases.
44


0.4
02
>â–  00 -0.2 -0.4
Figure V.2: Poincare Section for Phase Portrait of Stochastic Rayleigh-van der Pol model under different noise level, y = 0.1, A = 0.05, t > 1300, fort e[0, 2000].
Under the influence of relatively stronger noise (A = 0.05), figure 5.2 suggests that the system became way more unstable. We can also tell by the increasing of the width of Poincare section.
45


02fl 0.29 0.30 031 0.32 0.33 0.34 0.35
r
Figure V.3: Histogram of Poincare Section for Phase Portrait of Stochastic Rayleigh-van der Pol model under different noise level, g = 0. A = 0.01.
We obtained the histogram of Poincare section to observe the statistical property of the system. It is a equivalent to view the cross section of the Poincare section. Figure 5.3 showed that under the influence of relatively weak noise (A = 0.01), we can see an almost-normal distribution. The red line indicate the constant-amplitude
y_
solution where r= [J.
46


40
0.15 0.20 0.25 0.30 0.35 0.40
r
Figure V.4: Histogram of Poincare Section for Phase Portrait of Stochastic Rayleigh-van der Pol model under different noise level, g = 0. A = 0.05.
Figure 5.4 showed that in the case of relatively stronger noise (A = 0.05), the distribution had increased the bell width and shift to the left while be less normal distributed.
V.I.II Behavior of Mean and Variance with Bifurcation Parameter for Fixed Noise Level
We found that the variance does not influence by dt for our system with mu /=
Oven though for a stochastic process [ V ar] = c2 dt where C2 = A for our case.
After the data collecting, we computed the mean and variance of the Poincare
section versus jt/ to see if we can get any information for the bifurcation, as shown
from figure 11.12 to figure 11.17. As we can tell from the figures, for small noise level
(i.eA = 0.01), we can easily tell where the bifurcation happens (around jt/ = 0). From
the sudden increase of mean in mean versus jt/ and the peak in variance versus jt/, at
47


this noise level we can also observe the phase portrait for indications of bifurcation by changing [J to see the difference in phase portraits for limit cycle and fixed point. But for raised noise level (i.e. A = 0.05), it is difficult to find out where the bifurcation happens by the phase portrait, but the change of mean in mean versus jt/ and the peak in variance versus /J. For intensified noise level (i.e.A > 0.5), the mean versus jt/ cannot give us much indication about where the bifurcation happens, however the variance versus [J can still give us a region where the bifurcation happens.
The theory of Hopf bifurcation defines where bifurcation happens as a critical point. However, after observing the experimental data, it seems like for the case of bifurcation under the influence of the noise, the critical point disappears, a “critical region” or “bifurcation region” appears instead. The “bifurcation region” widen with a higher amplitude while shifting to the right instead of stay at [J = 0 as the stochastic term increases.
We have obtained the mean and variance of the Poincare section versus /J. It is worth noting that while we ran the simulation, we seeded the pseudo-random number generator to ensure we were studying the same sequence of random variables.
48


mean v.s. mu, 1 seed,lambda=0.01

1 n nnnnnnnn
. UU UUUU UU onnnnnnn • • •
u n . oUUUUUUU finnnnnnn 4 • •
u n .DUUUUUUU yinnnnnnn • •
u n .M-UUUUUUU onnnnnnn /
u . zuuuuuuu ’.oftoSoSo^ • 1 ft '
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(a) Mean v.s. n for A = 0.01.
Var v.s. mu, 1 seed,lambda=0.01
0 n .00070000
U n »
u n •
n • • •
u n .UUUJUUUU ^ • • 4
n • • • •
4 u .uuuxww .00000000 •*|
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(b) Variance v.s. /j for A = 0.01.
Figure V.5: Stochastic Rayleigh-van der Pol model, behavior of mean and variance,
for A =0.01, t> 1300.
49


For .A = 0.01, we can tell from figure 5.5 where the bifurcation happened by the sudden increase in mean versus jt/ as well as the peak in variance versus /J. The peak of variance is around jt/ = 0 where the bifurcation happens for the deterministic model.
50


mean v.s. mu, 1 seed,lambda=0.02

1 nnnnnnnn
I ri .UuUWJVUU Qnpuwinn • • •
U Pi .BUUlAluUU AnrwiAnn 4 • •
u ( .DUUUlAlUu Annnrwirwi • •
n innnnnnn /
u .ZUUOwUu IcficDoSoo
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(a) Mean v.s. n for A = 0.02.
Var v.s. mu, 1 seed,lambda=0.02
—f f .00250000 •
u PI »
u PI a • • • * «
u f 1 • •
U —t •t/Uu Jt/UtF® • * .00000000 • • 1 * * • • 1
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(b) Variance v.s. /j for A =0.02.
Figure V.6: Stochastic Rayleigh-van der Pol model, behavior of mean and variance,
for A =0.02, t> 1300.
51


For A = 0.02, we are still be able to tell from figure 5.6 where the bifurcation happened by the sudden increase in mean versus jt/ as well as the peak in variance versus /J. However, the “sudden change” in the mean versus jt/ is a little more smooth than in the care of A = 0.01. The variance also had a higher value at the peak and wider peak region. The peak of variance slightly shifts to the right side of jt/ = 0.
52


mean v.s. mu, 1 seed,lambda=0.05
—t 1 20000000 nnnnnnnn h
± r\ uuuuuuuu onnnnnnn • • •
U n oUUUUUUU cnnnnnnn • • •
u n OUUUUUUU nnnnnnnn • •
u n H-UUUUUUU onnnnnnn 1
4 ►#••••••* 0.00000000 '
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(a) Mean v.s. n for A = 0.05.
Var v.s. mu, 1 seed,lambda=0.05
0, A ,01000000 OAQflAAOA #
u, n UU-7UUUUV flAQAAAAA • •
u, UUoUUUUU

A AAfiAAAflA 4
U i A UUOUUUUU ' nntnnnnn * * 1
u, A UUjUUUUU ,004000001 nr\ jnnnnAB
U, A •
Ui ri UUOUUUUU* t nmnnnnn
Ui n uu zuuuuu • # fiAi nnnnn •*4 • •
4 ► • • • 0, AJ^r J.UUUUaJ 00000000 **.4 ►
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(b) Variance v.s. /j for A =0.05.
Figure V.7: Stochastic Rayleigh-van der Pol model, behavior of mean and variance,
for A =0.05, t> 1300.
53


In figure 5.7, for A = 0.05, we are still be able to tell where the bifurcation happened by the sudden increase in mean versus jt/ as well as the peak in variance versus /J. However, the “sudden change” in the mean versus jt/ is more smooth than in the care of A = 0.2. The variance also had a higher value at the peak and wider peak region compare with the case for A = 0.02. The peak of variance shifts to the right side of jt/ = 0 even more.
54


mean v.s. mu, 1 seed,lambda=0.1
1.20000000

i onrvwi/vi • • •
J.uUwUUw i fin/wtnnn i • •
J. OUi/UUUUU i ^nrvwwvi •
JihUUUUUUu 3.20000000 j 3.00000000 • »—

-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(a) Mean v.s. /j tor A = 0.1.
Var v.s. mu, 1 seed,lambda=0.1
1.03000000 •
i j.uz jUUU UU •
#
UJ. uUUU UU 1 m •
1 7 • •
• 3.00000000 •
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
(b) Variance v.s. /j for A = 0.1.
Figure V.8: Stochastic Rayleigh-van der Pol model, behavior of mean and variance, for A = 0.1, t> 1300.
For ^ = 0.1, as showed in figure 5.8, we had difficulty to tell where the bifurca-
55


tion happened by mean versus jL/, but we can tell the region where the bifurcation happened. We can still be able to tell where the bifurcation happened by the peak in variance versus /J. The variance also had a higher value at the peak and wider peak region compare with the previous cases. The peak of variance shifts to the right side of jL/ = 0 more significantly.
56


mean v.s. mu, 1 seed,lambda=0.7
2 00000000 ••
1 1 60000000 40000000 >T
• •**
•
1 1 ZUUUUUUU 00000000
J • • (•*
••
0* 20000000 00000000
—6 —0

-2.0000 -1.0000 0.0000 1.0000 2.0000 3.0000 4.0000
(a) Mean v.s. n for A = 0.7.
Var v.s. mu, 1 seed,lambda=0.7
0.20000000
0.18000000 •

0.16000000 •

0.14000000 ♦ • •

0.12000000 • 1 .*
• • •
0.10000000*« • o.osoooWJb Tpr n^nfinnnnnn
• ••
•

0.02000000 • • y • •*

—0.00000000
-2.0000 -1.0000 0.0000 1.0000 2.0000 3.0000 4.0000
(b) Variance v.s. n for A = 0.7.
Figure V.9: Stochastic Rayleigh-van der Pol model, behavior of mean and variance,
for A =0.7, t> 1300.
57


Figure 5.9 shows mean and variance versus jL/ for significantly bigger noise/s-tochastic term (A = 0.7). We cannot tell where bifurcation region is just by mean versusjt/. The peak region of the variance, which indicates the bifurcation region, continued to widen with increased peak value compare with the previous cases. There seems to be a “second” peak in the bifurcation region at fj = 0.7 as well. The surface plots in appendix A shows a better visualization for coupled system.
58


v
H *
mean v.s. mu, 1 seed,lambda=0.8
2.00000000 I
1.80000000 *•
*
1.60000000 «,*
*
1.40000000 _•
1.20000000
1.00000000 •*
0.30000000 -•* J- -
0,60000000 *****
ortJflflfiooo*]
0,20000000 0.00000000 I
-2.0000 -1.0000 0.0000 1.0000 2.0000 3.0000 4.0000
(a) Mean v.s. /j tor A = 0.8.
Varv.s. mu, 1 seed,lambda=0.8
0.25000000 ,
0.20000000
0.15000000
«
0.10000^
0,05000000
^%sv *
•*v.*
0.00000000
-2.0000 -1.0000 0.0000 1.0000 2.0000 3.0000 4.0000
(b) Variance v.s. /j tor A = 0.8.
Figure V. 10: Stochastic Rayleigh-van der Pol model, behavior of mean and variance,
for A =0.8, t> 1300.
59


For intensified noise/stochastic term (A = 0.8) showed in figure 5.10. We are still unable to tell where bifurcation region is just by mean versus jL/. The bifurcation region can still be found by the peak region of the variance, which is continued to widen with increased peak value like in previous case. It is more clear to see the “second peak” in the bifurcation which is around jL/ = 1 here.
60


CHAPTER VI
CONCLUSION AND FUTURE WORK
For the system of (modified) Rayleigh Van der Pol with additive stochastic term, the critical point disappears, a “critical region” or “bifurcation region” appears instead. The “bifurcation region” widens with a higher variance amplitude while shifting to the right instead of staying at jt/ = 0 as the stochastic term increases.
VI.I Future Work
One of the applications of Rayleigh Van der Pol Oscillator is in the simulation of the Wein Bridge circuit[30],
VI.I.I Coupled System
One of the applications of the coupled nonlinear system is in the simulation of the animal gaits [3] [4] [27],
VI.I.I.I Two Synaptic Coupled Stochastic Rayleigh-van der Pol Model
In the content of neuroscience, the neuron network can be represented in the form of coupled oscillators. There are different kind of coupling method such as synaptic coupling and diffusive coupling. We choose synaptic coupling since it is a more realistic way to represents how different neurons connect each other.
61


System of Two Synaptic Coupled Stochastic Rayleigh-van der Pol Equations:
dX'i = Y-i dt+ Y\ X2 dt+ dW'i(t), (VI.la)
c/Yi = ([j-Y ?-X )fYi dt -Xi dt + y2Y2 dt + dW2(f), (VI.lb)
dX2 = Y2 dt+ Y3 X\ dt + d\A/3(t), (VI. lc)
dY2 = ([j-Y i-X )\Y2 dt -X2 dt + y4 -Yi dt + dW4(f), (VI. 1 d)
Where dW^(t) = dW3(t) = 0, yy = y, / = 1, 2,3,4.
See appendix A for analysis and simulation for two coupled Rayleigh-van der Pol system.
62


REFERENCES
[1] Anishchenko Vadim S., Astakhov Vladimir, Neiman Alexander, Vadivasova Tatjana., Schimansky-Geier Lutz, Nonlinear Dynamics of Chaotic and Stochastic Systems: Tutorial and Modern Developments, Springer, 2007.
[2] Patrick Billingsley, Probability and Measure, Wiley Series in Probability and Statistics, 2012.
[3] Pietro-Luciano Buono, Martin Golubitsky , Models of central pattern generators for quadruped locomotion I. Primary gaits, J. Math. Biol, 2001.
[4] Pietro-Luciano Buono, Models of central pattern generators for quadruped locomotion II. Secondary gaits, J. Math. Biol, 2001.
[5] Kevin Burrage, Ian Lenane, Grant Lythe, Numerical Methods for Second-Order Stochastic Differential Equations, SIAM J. SCI. COMPUT Vol. 29, No. 1, pp. 245264,2007.
[6] T. K. Caughey, Response of Van der Pol’s Oscillator to Random Excitation, lournal of Applied Mechanics, 1959.
[7] T. K. Caughey, On the Respones of a Class of Self-excited Oscillators to Stochastic Excitation, Int. I. Non-Linear Mechanics, 1967.
[8] T. K. Caughey, Nonlinear Theory of Random Vibration, Advances in Applied Mechanics Volume 11, 1971.
[9] Andrey Dmitrieva, Victor Dmitrieva, Oleg Sagaydakb, Olga T sukanovaa, The Application of Stochastic Bifurcation Theory to the Early Detection of Economic Bubbles, Procedia Computer Science, 2017.
[10] Thai Son Doan, Maximilian Engel, Jeroen S W Lamb, Martin Rasmussen, Hopf bifurcation with additive noise, Nonlinearity 31 4567,2018.
[11] Arinin Fuchs, Nonlinear Dynamics in Complex Systems, Theory and Applications for the Life- Neuro- and Natural Sciences, Springer, 2013.
[12] M. A. Garstens, Noise in Non-Linear Oscillators, J. Appl. Phys., 1957.
[13] Kiyosi Ito, Stochastic Integral, Proc. Imp. Acad., 1944.
[14] Eugene M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, The MIT Press, 2007.
63


[15] Kurt Jacobs, Stochastic Processes for Physicists: Understanding Noisy Systems, Cambridge University Press, 2010
[16] Samuel Karlin, Howard M. Taylor, A First Course in Stochastic Processes, Academic Press, 1975.
[17] Tatsuyuki Kawakubo, Shigeharu Kabashima, Stochastic Process in Self-Excited Oscillation, Journal of Physical Society Japan, 1974.
[18] Peter Eris Kloeden, Eckhard Platen, Henri Schurz, Numerical Solution of SDE Through Computer Experiments, Springer, 2003.
[19] Christian Kuehn, A Mathematical Framework for Critical Transitions: Bifurcations, Fast-Slow Systems and Stochastic Dynamics, Physica D: Nonlinear Phenomena, Vol. 240, No. 12, pp. 1020-1035, 2011.
[20] H.K. Leung, Stochastic Transient of a noisy Van der Pol Oscillator, Physica A, 1995.
[21] H.K. Leung, Stochastic Hopf Bifurcation in a biased Van der Pol Model, Physica A, 1998.
[22] Carlo Liang, Gabriel J. Lord, Stochastic Methods in Neuroscience, Oxford University Press, 2010.
[23] Marc Mendler, Johannes Falk, Barbara Drosse , Analysis of stochastic bifurcations with phase portraits, arXiv, 2018.
[24] John F. Moxnes, Kjell Hausken, Introducing Randomness into First-Order and Second-Order Deterministic Differential Equations, Advances in Mathematical Physics Volume 2010, 2010.
[25] Bernt Oksendal, Stochastic Differential Equations (6th Ed.): An Introduction with Applications, Springer, 2010.
[26] Sheldon M. Ross, Stochastic Processes, John Wiley and Sons, Inc., 1996.
[27] Mozhdeh Saffari-Parizi, Coupled Biological Oscillators and Animals Gaits,
MIS Thesis, University of Colorado Denver, 2004.
[28] P-T. D. Spanos, Stochastic Analysis of Oscillators with Non-Linear Damping,
Int. J. Non-Linear Mechanics, 1978.
[29] Steven H. Strogatz, Nonlinear Dynamics and Chaos (2nd Ed.): With Applications to Physics, Biology, Chemistry, and Engineering, Westview Press, 2014.
[30] Randall Tagg, UC Denver Physics Department, Private Communication.
[31] Ruoshi Yuan, Xinan Wang, Yian Ma, Bo Yuan, Ping Ao , Explor-
ing a noisy van der Pol type oscillator with a stochastic approach, PHYSICAL REVIEW E 87, 062109, 2013.
63


APPENDIX A
ANALYSIS AND SIMULATION FOR COUPLED SYSTEM
This appendix contains the analysis and simulation for two coupled Rayleigh-van der Pol system.
For the coupled case, we are studying
- coupling strength y v.s. jt/ for different noise level A, find where the bifurcation happens;
- mean and variance of the Poincare section, find where the bifurcation happens; -surface plot of mean and variance of each oscillator v.s. y and jt/ for fixed noise
level (i.e. ^2,4 = 0.05).
The results for the first oscillator in the system shows in the figures below.
64


Gamma v.s. mu, Iambda2,4=0
0.15
--------------------- •----04—•--
--1--- -#--- m • ojos i • •• •
• Fix Point
--•-------•-----------•—
-0.4 -a3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 • Limited Cycle
m • -0.05 • • • • • •
---•--- --------•------•----0.1 •---•---------•- ■> HIM---------•--
-0.15
Figure A. 1: Coupled System of Stochastic Rayleigh-van der Pol, y v.s.p, A2A = 0.
For the deterministic twp-coupled Rayleigh-van der Pol equations, we can see from the Figure A.l that the bifurcation happened at y = ~2p \4hich is one of the eigenvalues for the system.
65


Gamma v.s. mu, Iambda2,4=0.01
---OAS—|--
• •••••••••••« • 0.1 4 • • m
• ---•-- «••••• 0.05 0 --•------•-----
• Fix Point
--•-------•---IH i > t i------•--------------------------- • Limited Cyde
-a4 -0.3 -a2 -ai o o.i 0,2 oj
• • • • -0.05 • • ••••••
• • • • • • • • -0A-* • • •-------------------•
-0 I 5
Figure A.2: Coupled System of Stochastic Rayleigh-van derPol, y v.s.p, K2a = 0.01, the gray points indicate the bifurcation region.
The intriguing part happened once we introduce stochastic term into the coupled system as showed in figure A.2. Even for A = 0.01, the bifurcation transit from a line to a region, the gray points marked as “?” are the points in the bifurcation region. The bifurcation region is a combination of positive part of y = -yjL/ and negative part of y = yV both shift to the right of axis. Moreover, the bifurcation region widened more for positive y.
66


• Fix Point
• Limit Cycle
Figure A.3: Coupled System of Stochastic Rayleigh-van derPol, y v.s.p, A2A = 0.05.
For .A = 0.05 which showed in figure A.3, the bifurcation region is still a combination of positive part of y = - and negative part of y = The bifurcation region keep shifting to the right of axis while widened more for positive y than negative y.
Gamma v.s. mu, Iambda2,4=0.05
0.15 ,
67


Gamma v.s. mu, Iambda2,4=0.1
0.1S
-0.4
• Fix Points
• Limited Cycle
-0.15
Figure A.4: Coupled System of Stochastic Rayleigh-van der Pol, y v.s.p, A2A =0.1.
For .A = 0.05 which showed in figure A.4, the bifurcation region is still a combination of positive part of y = - and negative part of y = The bifurcation region keep shifting to the right of axis while widened more for positive y than negative y.
68


meanl v.s.mu, Iambda2,4=0.05
• gamma = -0.1
• gamma = -0.05
• gamma = 0
• gamma = 0.05
• gamma = 0.1
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
4 * < • Ah » ►
( • • • • • • * • •
0 60000000 • • • •4 • • • ► • »

O 40000000 • • • 1*.
0.20000000 1 •

• • « ,*8 •* 0.00000000 »
Figure A. 5: Coupled System of Stochastic Rayleigh-van der Pol, mean of the first oscillator v.s. p for different y, Ma = 0.05.
Now we fixed the stochastic term(i.e. A = 0.05), we want to see how the statistical properties behave for different coupling strength y.
Figure A. 5 suggests that only the magnitude of y will influence the mean of the system. Increasing the magnitude of y, has the opposite effect as increasing A. The increase of the magnitude of y make the “sudden change” of mean sharper and the value of mean bigger while shifting the bifurcation point/region to the left of the axis.
69


Varl v.s.mu, Iambda2,4=0.05
â–  . 0.01600000
0.01400000 •
0.01200000
o 01000000
0.00800000
« i 0.00600iDO » >
• • n mvinnnnn 1

• 1 • ! *0.00200000 ^ • *: • • ♦
« ! « * nnnnnnnnn • : • 1 • • • ! t
' 1 1 W—
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
• gamma = -0.1
• gamma = -0.05
• gamma = 0
• gamma = 0.05
• gamma = 0.1
Figure A. 6: Coupled System of Stochastic Rayleigh-van derPol, variance of the first oscillator v.s. p for different y, Ma = 0.05.
Meanwhile, Figure A.6 suggests similar tendency, that only the magnitude of Y will influence the mean of the system. Increasing the magnitude of y, has the opposite effect as increasing A. The increase of the magnitude of y make the decrease the amplitude of the peak region of bifurcation region while narrowing and shifting the bifurcation region to the left of the axis.
70


meanl
-1.0
0.10 005 0.00 -0.05 -0.10
gamma
Figure A. 7: Coupled System of Stochastic Rayleigh-van der Pol, surface plot of the first oscillator for mean, p and y, A2A = 0.05.
meanl
1.0
mu
Figure A. 8: Coupled System of Stochastic Rayleigh-van der Pol, surface plot of the first oscillator for mean, p and yata different angle, A2A = 0.05.
Figure A. 7 and A. 8 shows the surface plots for mean, p and y of the first oscillator at different angles. The surface plots give a better visualization about mean, p, y and bifurcation at fixed stochastic terms (A2,a = 0.05): mean increases as the magnitude of y increase, the “sudden change” in mean indicate where the bifurcation region is.
71


Varl
0014
0012
0010
-0 008 â– I- 0 006
i C 034 â–  C
Figure A. 9: Coupled System of Stochastic Rayleigh-van der Pol, surface plot of the first oscillator for variance, p and y, A2A = 0.05.
Varl
-0 008
10006 0004 0002
Figure A. 10: Coupled System of Stochastic Rayleigh-van der Pol, surface plot of the first oscillator for variance, p and y at a different angle, A2A = 0.05.
Figure A.9 and A. 10 shows the surface plots for variance, p and y of the first oscillator at different angles. The surface plots give a better visualization about variance, p y and bifurcation at fixed stochastic terms (A2,a = 0.05): variance decreases as the
72


magnitude of y increase, the peaks indicate the bifurcation region. There seems to have a “main bifurcation” and “sub bifurcation” inside the bifurcation region.
mean2
-l.o
010 005 0.00 -0.05 -0.10
gamma
Figure A. 11: Coupled System of Stochastic Rayleigh-van derPol, surface plot of the second oscillator for mean, p and y, A2a = 0.05.
mean2
mu
-100 075 - c sc
Figure A. 12: Coupled System of Stochastic Rayleigh-van derPol, surface plot of the second oscillator for mean, p and y at a different angle, A2a = 0.05.
Figure A.l 1 and A. 12 shows the surface plots for mean, p and y of the second oscillator at different angles. The surface plots give a better visualization about
73


mean, y and bifurcation at fixed stochastic terms (A2,4 = 0.05): mean increases as the magnitude of y increase, the “sudden change” in mean indicate where the bifurcation region is.
Var2
0010
0008
0006
0004
0002
-0 10
-0 5
U U3
mu
C 1C
- 0 006
I 0;:4 A- C
Figure A. 13: Coupled System of Stochastic Rayleigh-van derPol, surface plot of the second oscillator for variance, p and y, A2A = 0.05.
Var2
- 0 006 - C 034
! I C
Figure A. 14: Coupled System of Stochastic Rayleigh-van derPol, surface plot of the second oscillator for variance, p and y at a different angle, A2a =0.05.
Figure A. 13 and A. 14 shows the surface plots for variance, p and y of the second
74


oscillator at different angles. The surface plots give a better visualization about variance y and bifurcation at fixed stochastic terms (^2,4 = 0.05): variance decreases as the magnitude of y increase, the peaks indicate the bifurcation region. There seems to have a “main bifurcation” and “sub bifurcation” inside the bifurcation region and both bifurcation regions are a little different compare to the first oscillator.
The coupled system behaves a little different than we anticipated, the nature of the coupled stochastic nonlinear system and its behavior is not well understand yet. As we can tell from the yv.s.mu plots, the bifurcation region has an excepted shift once the randomness has been introduced to the system, then the bifurcation region widens and shift to the right as the stochasticity increases. For the fixed A, both positive and negative coupling shift the bifurcation region to the left, the bifurcation narrowed as the absolute value of y increases.
75


Full Text

PAGE 1

MODELING THE RESPONSE OF RAYLEIGH VAN DER POL OSCILLATORS TO STOCHASTIC EXCITATION NEAR THE HOPF BIFURCATION by BINGSHEN LU B.S., University of Colorado Denver, 2015 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Integrated Sciences Integrated Sciences 2019

PAGE 2

iii This th esis for the Master of Integrated Sciences degree by Bingshen Lu has been approved for the Integrated Sciences Program by Masoud Asadi Zeydabadi, Chair Randall Tagg Burt Simon Date: May 18, 2019

PAGE 3

iv Lu, Bingshen (M.S., Integrated Sciences) Modeling the Response of Rayleigh Van der Pol Oscillators to Stochastic Excitation Near the Hopf Bifurcation Thesis directed by Associate Professor Masoud Asadi Zeydabadi ABSTRACT In this thesi s we investigate the Rayleigh van der Pol oscillator with additive noise. This dynamical system is a fundamental model that can be tested experimentally by the Wien bridge electrical circuit. The analysis of the Rayleigh van der Pol oscillator in the presence of noise provides a good understanding of the stochastic nonlinear dynamical system. The Rayleigh van der Pol oscillator is a two variable nonlinear system with one control parameter. As the control parameter varies at the bifurc ation (transition) point the state of system changes from steady to self sustained oscillation. This transition is known as Hopf bifurcation. The focus of this project is on the behavior of the system near bifurcation point in the presence of noise. We fou nd that the bifurcation point widens to a shifts as the level of noise increases. The form and content of this abstract are approved. I recommend its publication. Approved: Masoud Asadi Zeydabadi

PAGE 4

TABLEOFCONTENTS IINTRODUCTION...............................1 IIESSENTIALIDEASFROMNONLINEARDYNAMICS..........7 II.ILinearSystems...............................7 II.IISimpleHarmonicMotionwithDamping.................9 II.IIILimitCycleandHopfBifurcation.....................11 II.IVANonlinearSystem:theClassicVanderPolModel...........11 II.VModiedVanderPolModel........................18 II.VIRayleigh-vanderPolModel........................21 IIIESSENTIALIDEASFROMSTOCHASTICDIFFERENTIALEQUATION SDE......................................27 III.IClassicalExampleofStochasticProcess:BrownianMotion.......27 III.IIGeneratingaStochasticProcessNumerically...............27 III.IIIEmbeddingStochasticProcessesinDierentialEquations........28 III.IVNumericalMethodsforSDE........................28 III.IV.IEuler{Maruyamamethod........................29 III.VSolvingStochasticIntegrals........................29 III.V.IIt^o'sRule.................................29 III.V.IIIt^o'sChainRule.............................30 III.V.IIIAnalyticalexamples...........................30 IVSTOCHASTICDYNAMICS..........................37 IV.IStochasticVanderPolModel.......................37 IV.IIStochasticRayleigh-vanderPol......................37 VBEHAVIORNEARTHEBIFURCATION..................42 V.INumericalSimulation............................42 V.I.IPoincareSection.............................42 v

PAGE 5

V.I.IIBehaviorofMeanandVariancewithBifurcationParameterforFixed NoiseLevel..........................46 VICONCLUSIONANDFUTUREWORK...................60 VI.IFutureWork.................................60 VI.I.ICoupledSystem.............................60 REFERENCES ...................................62 APPENDIX A.ANALYSISANDSIMULATIONFORCOUPLEDSYSTEM........64 vi

PAGE 6

1 dX 1 CHAPTER I INTRODUCTION Bifurcation, a sudden change of the dynamical system, lies in the center of nonlinear dynamics[29][19]. It has plenty of applications in daily activity and fields of science such as the quick change of medical condition from stable to severe, financial mar ket change from balanced to crisis[19], or neural network change from stationary to generates motion. The needs to do a proper prediction of such phenomena rise as time goes by. However, reality has uncertainties so it is important to understand how ran domness effects the nonlinear systems. The general theory of random excitation in physical systems was first studied by Albert Einstein in 1905 and was generalized by Marian Smoluchowski in 1916[8]. Later on, physicists and mathematicians Andrey Kolmogorov, Adriaan Fokker, Max Planck, Leonard Ornstein, George Eugene Uhlenbeck, Norbert Wiener, Joseph L. D o ob and Ki y oshi I t o laid the foundation for solving linear SDE b e t w een 1930s and 1950s[8]. However, those earlier works are for solving linear SDE. Difficulties arise for nonlinear system because of the complexity nature of nonlinear systems[8]. Back in 1957, Garstens started studying noise in nonlinear oscillators[12]. Garstens used a nonlinear Van der Pol equation with driving noise term to model nonlinear s in such model are multiplicative (i.e. X 1 · X 2 , where dt = X 2 ). Garstens presented a method

PAGE 7

2 for estimating the effect of nonlinear noise in an oscillator at low levels of oscillation. This required an approximate solution of a non homogeneous van der Pol type equa tion with a noisy driven force by reducing the equation to a linear form. The powe r spectrum of the output was used to study the effect of noise in the nonlinear system. Caugheys paper in 1959[6] showed that there was interest in the how nonlinear systems response to random excitation, especially in the systems with multiplicative nonli near terms. He studied the self excited system like Van der Pol Oscillator and its response to external stochastic terms acting as external random excitation. Caughey found that the solution has a stochastic/noise term because of the stochastic term i n the self excitation system, and that the strength of the stochasticity decreases as the amplitude of oscillation increase. Interestingly enough, so does the band width of the stochasticity/noise which prove the early work done by Garstens incorrect. In C augheys paper in 1967[7], the use of a Fokker Planck equation had been introduced to study the stationary response of self excited oscillators. Using such technique, Caughey obtained the exact solution for the stationary joint probability density function, as well as the spectral density for oscillator with weak damping and weak excitation. Following his previous work, Caughey developed a nonlinear theory of random cause of t he developments in fields like high speed flight and the needs of dealing the random fluctuation cause by turbulence back at that time. Caughey summarized and developed some solution techniques for random vibrations in nonlinear systems[8]. During the simi lar time in Japan (1974), Kawakubo and Kabashima also used the Fokker Planck approach to study stochastic processes in self excited oscillation using a Van der Pol stochastic nonlinear system in a more experimental way[17]. Kawakubo and Kabashima focused on how the fluctuations, namely the amplitude and phase

PAGE 8

fluctuations in the system, behave in the experiment of Wien Bridge oscillators. They showed how the voltage fluctuates with the feedback factor and how the variance fluctuates with the parameter controlling the nonlinearity. Kawakubo and Kabashima also included the probability distribution function obtained with the Fokker Planck equation, examining how it changed with the parameter controlling the nonlinearity. In agreement with Caughey and others find. Spanos showed that in 1978 such problems were drawing more and more attention[28]. This included several more applications of stochastic phenomena i n physics and engineering, such as control, aerospace, natural disaster, ocean wave and transportation engineering. Spanos pointed out that even though linear random vibration has be well studied, the theory of non linear random vibrations was not. This is still the case: there are still a lot of development needed in this topic. The Fokker Planck approach was used in Spanos work to obtain solutions, Spanos also showed the application of Fokker Planck approach in Van der Pol and Rayleigh oscillators. Later on in 1995, Leung did a systematic study of the behavior of Van der Pol oscillator with additive and multiplicative stochastic terms using different nonlinear ity parameters[20], comparing the bifurcation behaviors of the system with di fferent stochastic terms. By such methods Leung stated that we could get a more realistic understanding of the bifurcation of a nonlinear system. Leung found that all level of stochastic terms / noise chosen disturb the usual behavior of the Van der Pol o scillator. Moreover, he found that noise could postpone the bifurcation between fixed point and limit cycle and that noise could induce a transition from a stable fixed point to a stable limit cycle, and vice versa. Leung suggested that noisy nonlinear sys tems could have more and different phenomena than the deterministic one, such as stochastic resonance. 3

PAGE 9

In later work in 1998[21], he pointed out studies on stochastic bifurcation focus on the shift of bifurcation point, where a s hifting of bifurcation point is an indication of the noise effect on the two competing mechanisms, namely self sustained oscillation and nonlinear dissipation. Yuan and others made further developments in this direction in 2013 and found that zero mean multiplicative noise can induce the bifurcation points and sometimes even change the topology of the system[31], or in another words, the structure does not correspond to the deterministic counterpart, which make the analysis of such system more difficult. They also used the steady state distribution of the process to approach such case, in where correspond means the probability density function does not change over time and have a maximum at the stable attractors, which later on define as t he steady state distribution multiplying negative one being a Lyapunov function of the deterministic counterpart of the stochastic dynamics. Yuans group apply a new stochastic interpretation rather than the common Ito and Stratonovichs integration to ob tain a correspondence between stochastic and deterministic dynamics. Just last year (2018), Doan and others were trying to merge bifurcation theory for stochastic dynamical systems[10]. Doan and others work shows that the earlier concepts developed in t he 1990s, namely theory of phenomenological (P) bifurcation and dynamical (D) bifurcation, cannot comprehensively capture the complexity of bifurcation in random dynamical system. Doan points out that despite its relevance for many applications of topical interest, a bifurcation theory of random dynamical systems is still in its infancy[10]. With the developments made in mid 20th century on time dependent probability theory[2][16][26] and stochastic differential equations[25], and we began to have tools to solve the stochastic system in the form of stochastic differential equation[13][25]. Yet The bifurcation of stochastic chaotic system is still something we do not suffi ciently understand yet. But the importance of bifurcation in such system is becoming 4

PAGE 10

5 dX 1 increasingly unavoidable[1][10]. In recent years, neuroscience (as well as machine learning and artificial intelligent) draws more and more attention[22] and funding, so do the nonlinear and stochastic components of this field[1][22][15]. There are some breakthroughs on the applications of such topics in crucial areas such as neuroscience[14][11][17]. However, most of the developments about stochastic differential equations are in the field of stock market forecasting and related topics in quantitative/financial engineering[25], i.e. the early detection of economic bubbles[9]. How to analyze the stochastic nonlinear dyn amical model is still a difficult task which been rarely discussed in the context of stochastic differential equation but mostly done by deterministic approximation, especially for nonlinear system with multiplicative terms in neural networks (i.e. System contains X 1 · X 2 , where dt = X 2 ). Most of the analytical and numerical methods cannot be applied to such terms[5][18][25]. Yet the interest in self excited oscillation and its application in neural network pushes us forward to study such system[3][17], especially the response of a class of oscillators to stochastic excitation[28]. Now we might have tools and methods to study it in the content of stochastic di fferential equation instead of deterministic approximation methods. In this thesis we are building a modified Rayleigh Van der Pol model with additive stochastic components. We hope to test if the system behaves as others observed, when we put the system in the proper context of stochastic differential equation in stead of deterministic methods. We also hope to lay a foundation for future work and potential applications of the Rayleigh Van der Pol model. The applications include neural networks for the Rayleigh Van der Pol[6][7] based central pattern generator (CPG) model[3], as well as other potential applications i n physics, neuroscience, fi nance and other related fields. This work is also offered as a guide for some methods of how to simulate and analyze models of nonlinear system with stochastic components in the context of

PAGE 11

6 dX 1 stochastic differential eq uations. This is primarily focused on how to solve the system with terms containing random variable products X 1 · X 2 , where dt = X 2 , which we do not have many good options to solve yet.

PAGE 12

7 CHAPTER II ESSENTIAL IDEAS FROM NONLINEAR DYNAMICS This chapter organizes a conceptual framework for analyzing two dimensional non linear dynamics. This provides a broad and powerful foundation for analyzing many applications especially after stochastic forcing is added in later chapters. II.I Linear Systems In general a two dimensional dynamical linear system has the form of x 1 = a 11 x 1 + a 12 x 2 (II.1a) x 2 = a 21 x 1 + a 22 x 2 (II.1b) where x 1 = dx , x 2 = dx 2 , and a 11 , a 12 , a 13 , a 14 are constant parameters. dt dt We can write the system in matrix form as x = A x (II.2) where A = a 11 a 1 2 and x = x 1 . a 21 a 22 x 2 We seek solution of the form [29],

PAGE 13

8 x = e v. (II.3) Substituting (2.3) into (2.2), we get Av = (II.4) We say the solution in (2.3) is an eigensolution if v is an eigenvector of A with corresponding eigenvalue [29]. The characteristic equation of matrix A is given by det a 11 a 12 = 0 , (II.5) Equation (2.5) yields a 21 a 22 where 2 0 (II.6) = trace ( A ) = a 11 + a 22 (II.7a) = det ( A ) = a 11 a 22 a 12 a 21 (II.7b) Solving equation (2.6), we get 1 = 2 = + 2 , (II.8a) 2 2 . (II.8b) 2 where 1 and 2 are the eigenvalues of matrix A . We find the eigenvectors by substituting eigenvalues into one of the equations in (2.4). We get two linearly inde pendent eigenvectors v 1 and v 2 . Now we can get the general solution for x , x = c 1 e 1 t v 1 + c 2 e 2 t v 2 . (II.9)

PAGE 14

9 0 m m II.II Simple Harmonic Motion with Damping One good example for the linear system is simple harmonic motion with damping. The F in this case is the sum of the spring and damping forces: dx F = kx b dt , (II.10) where k and b are the spring and damping for (2.10), we get d 2 x dx m = kx b (II.11) dt 2 dt We can rewrite equation (2.11) as d 2 x dx m + b + kx = 0 (II.12) dt 2 dt or x ¨ + x + 2 x = 0 ( I I.13) where = b 0 = l k . The solution (2.12) and (2.13) has the form x = Ae (II.14) Substituting (2.14) into (2.12), then we have ( 2 + + k ) Ae = 0 (II.15) Equation (2.15) for the nontrivial solution of (2.12) implies 2 + b + k = 0 . (II.16) m m Equation (2.16) is the characteristic equation of (2.12).

PAGE 15

10 4 m 2 m = 2 m ± i m 4 m 2 (II.17) 4 m 2 m x = A 1 e 2 m te m 4 m 2 + A 2 e 2 m m 4 m 2 The roots of equation (2.16), eigenvalues of (2.12), for underdamped or oscillatory case, ( b 2 < k ), is given by b I k b 2 The general solution of (2.12) or (2.13) for weak damping, b 2 < k , is given by b i I k b 2 t b te i I k b 2 t or x = A e t e i w 2 2 / 4 t + A e t e i w 2 2 / 4 t ( I I.19) 1 0 2 0 For real valued solutions, A 1 and A 2 will be complex conjugates. We call x = 0 , x = 0 a fixed point solution. For 0 all solutions tend to this fixed point solution and it is therefore said to be stable, for 0 the fixed point is unstable. Figure (2.1) shows time series and phase plot for damped harmonic motion. Figure II.1: Underdamped Harmonic Motion (II.18)

PAGE 16

11 dt µ 2 This is the standard, or classical form of the van der Pol equation. However, the µ II.III Limit Cycle and Hopf Bifurcation A limit cycle occurs when the phase plane ( x v.s. dx ) contains an isolated closed trajectory. Isolated means the neighboring trajectories either spiral towards or away from the limit cycle.[29] A Hopf bifurcation h appens where a fixed point of a system changes to a limit cycle. II.IV A Nonlinear System: the Classic Van der Pol Model The Van der Pol equation has been used as a model for neuron networks with feedback controls. It is named after Dutch physicist Balthasar va n der Pol[29][30], who applied the equation to electronic circuit oscillators. d 2 x 2 dx 2 dt 2 ( µ ) 0 dt + 0 x = (II.20) Provided that µ / = 0 we can write: d 2 x 2 dx 2 dt 2 µ (1 µ x ) 0 dt + 0 x = 0 (II.21) Using the rescaled variables, t 0 t and x l x , equation (2.21) becomes: d 2 x d t 2 dx µ (1 x ) d t + x = 0 (II.22) scale factor l becomes singular at µ = 0 at the same point that the actual output, x , has an amplitude that is just about to rise from a value of zero[30]. This is one of the reasons that we modified the classical Van der Pol model as shown in the later section.(see section 2.3 and 2.4) An alternate dimensionless form of the dynamics is the single scaled nonlinear

PAGE 17

12 equation in x :

PAGE 18

13 x 1 l dt l d 2 x 2 dx dt 2 ( µ x ) dt + x = 0 . (II.23) Here, the scaled time and dynamic variable are t 0 t and x x where x . For notational clarity, we the equations are still in the scaled form. We can applying method of averaging to the Van der Pol equation to find ap proximate solutions to the nonlinear system. The dynamics, equation (2.23), can be re written as two first order differential equations dx = w (II.24a) dt dw 2 dt = x + ( µ x ) w. (II.24b) We are going to assume that the the solution will evolve from the simple sinsoidal solution when the nonlinear term is turned on. Thus we think of time variations in amplitude and phase around the limit cycle such that the phase space trajectory can be written in the form: x ( t ) = A ( t ) cos[ t + ( t )] (II.25a) w ( t ) = A ( t ) sin[ t + ( t )] (II.25b) Our goal is to find equations for the amplitude function A ( t ) and the phase function ( t ). In order to do so, we insert (2.25a) and (2.25b) for x and w into the dynamical equations (2.24a) and (2.24b) to obtain dA ( t ) ( t ) cos[ t + ( t )] A ( t ) 1 + sin[ t + ( t )] = A ( t ) sin[ t + ( t )] (II.26) dA ( t ) ( t ) sin[ t + ( t )] A ( t ) 1 + cos[ t + ( t )] = dt dt A ( t ) cos[ t + ( t )] µ A 2 ( t ) cos 2 [ t + ( t )] 1 A ( t ) sin[ t + ( t )] . (II.27) dt

PAGE 19

14 2 2 2 2 2 2 2 1 2 2 1 2 dA ( t ) ( t ) dt dt 2 2 2 2 3 This simplifies to dA ( t ) cos[ t + ( t )] A ( t ) dt ( t ) sin[ t + ( t )] = 0 (II.28) dt dA ( t ) sin[ t + ( t )] A ( t ) dt ( t ) cos[ t + ( t )] = dt µ A ( t ) cos [ t + ( t )] 1 A ( t ) sin[ t + ( t )] . (II.29) Use (2.28) to write ( t ) dt dA ( t ) cos[ t + ( t )] = dt . (II.30) A ( t ) sin[ t + ( t )] Substitute (2.30) into (2.29) and multiply through by sin[ t + ( t )]: dA ( t ) 2 dA ( t ) 2 sin [ t + ( t )] dt cos [ t + ( t )] = dt E A ( t ) cos [ t + ( t )] 1 A ( t ) sin [ t + ( t )] . (II.31) Collect terms and use the fact that cos 2 [ t + ( t )] + sin 2 [ t + ( t )] = 1 to obtain dA ( t ) dt = µ A ( t ) cos [ t + ( t )] A ( t ) sin [ t + ( t )] . (II.32) We can also get ( t ) dt = µ A ( t ) cos [ t + ( t )] cos[ t + ( t )]sin[ t + ( t )] . (II.33) So far this has been an exact change of coordinates. We now make an approxi mation. Anticipating that A 2 ( t ) O ( µ ), we take both time derivatives and to likewise be O ( µ ) and thus the amplitude and phase are slowly varying deviations from the circular lim it cycle. We will then treat them as effectively constant over one cycle and replace the right hand side with averages. Taking t + ( t ), dA ( t ) 1 2 dt 2 µ A ( t ) cos ( t ) sin (II.34a) 1 2 1 2 0 0 0 2 2

PAGE 20

15 2 8 µA ( t ) 2 sin A ( t ) 2 cos sin (II.34b) 1 1 3 = µA ( t ) A ( t ) . (II.34c)

PAGE 21

16 1 So we have dA ( t ) 1 2 dt 8 4 µ A ( t ) A ( t ) . (II.35) This is an example of an *amplitude equation*. It may be formally integrated starting with to arrive at dA 8 = dt (II.36) [4 µ A 2 ] A A + 2 µ µt µ = e . (II.37) At long times A 2 (2 µ A ) A 2 µ (II.38) (justifying thies assumption that A 2 O ( µ )), while the next order correction is µ ( 2 e µt ) . (II.39) Returning to dimensioned coordinates, x ( 2 e 0 t ) x µ cos( 0 t ) , (II.40) and the predicted final steady oscillation is given by x 2 x µ cos( 0 t ) . (II.41) where we recall that x . For future purpose, we re write the classcial Van der Pol system as two first order equations, dx = y (II.42) dt

PAGE 22

17 dy 2 dt = µ · (1 x ) · y x (II.43) As we mentioned the Hopf bifurcation happens where a fixed point changes to a limited cycle. The fixed point of classic Van der Pol model is ( x , y ) = (0 , 0). In this system, the Hopf bifurcation happens at µ = 0, when the system from steady state figure (2.3) turn into a self sustain system (2.2). Figure II.2: Classical Van der Pol model, Stable Fixed Point, X versus t for µ = 0.1.

PAGE 23

18 Figure II.3: Classical Van der Pol model, Stable Fixed Poi nt, Y versus t for µ = 0.1. Figure II.4: Classical Van der Pol model, Stable Fixed Point, Phase Portrait for µ = 0.1. We can see from both the figures of X and Y versus t that the model started oscillating and damping to zero after some time. The phase portrait shows the red dot, which is our initial point, spiral in and become a fixed point, that indicates the model reaches a fixed point.

PAGE 24

19 Figure II.5: Classical Van der Pol model, Limit Cycle, X versus t for µ = 0.1. Figure II.6: Classical Van der Pol model, Limit Cycle, Y versus t for µ = 0.1.

PAGE 25

20 Figure II.7: Classical Van der Pol model, Limit Cycle, Phase Portrait for µ = 0.1. We can see from both the figures of X and Y versus t that the model oscillates at an fixed amplitude after some time. The phase portrait shows the red dot, which is our initial point, spiral out and become a circular trajectory, that also indicate the model reaches a limit cycle. II.V Modified Van der Pol Model The classical Van der Pol is a good simple model for the current goal. However, during the numerical simulation, the simulation breaks when µ is close to zero, which agree with the discussion in section 2.2. Therefore, we are using one version of modified Van der Pol model, dx = y (II.44) dt dy 2 dt = ( µ x ) · y x (II.45) Comparing the second equation in classical with modified Van der Pol equations, (comparing equations (2.43) with and (2.45)) we find µ is outside the parenthesis in (2.43) but is inside in (2.45).That makes the nonlinear term more smooth (with no singularity) and more close to real world phenomena. For this reason it i s easier to

PAGE 26

21 more consistent and does not break during simulation when µ is close to zero. For modified Van der Pol Model, the bifurcation still happens at µ = 0, as figures 2.8 2.13 indicate. Figure II.8: Modified Van der Pol model, Stable Fixed Point, X versus t for µ = 0.1. Figure II.9: Modified Van der Pol model, Stable Fixed Point, Y versus t for µ = 0.1.

PAGE 27

22 Figure II.10: Modified Van der Pol model, Stable Fixed Point, Phase Portrait for µ = 0.1. Given the same initial condition and parameters, we can see that the modified Van der Pol behaves almost the same with the classical counterpart. However, the modified Van der Pol model reaches fixed point faster than the classical model because the nonlinear term µ is less destructive. Figure II.11: Modified Van der Pol model, Limit Cycle, X versus t for µ = 0.1.

PAGE 28

23 Figure II.12: Modified Van der Pol model, Limit Cycle, Y versus t for µ = 0.1. Figure II.13: Modified Van der Pol model, Limit Cycle, Phase Portrait for µ = 0.1. Given the same initial condition and parameters, we can see as well that the mod ified Van der Pol behaves almost the same with the classical counterpart. However, the modified Van der Pol model reaches limit cycle faster than the classical model because the nonlinear term µ is less destructive. II.VI Rayleigh van der Pol Model We use (Modified) Rayleigh van der Pol equations as our standard system. One of the characteristics worth noting about Rayleigh van der Pol equations is that its

PAGE 29

24 / phase portrait is circular, even for larger value of µ which make our future work on the phase plane a lot easier. dx = y (II.46) dt dy 2 2 dt = ( µ y x ) · y x (II.47) Rayleigh van der Pol is also a good representation of the Wein bridge circuit which can be conduct into experiment for testing [30]. Now we change the system to polar coordinate, x = r cos (II.48a) y = r sin (II.48b) r = x 2 + y 2 (II.49a) y = arctan x (II.49b)

PAGE 30

25 / dt dt / r r 2 2 2 2 dr 1 dx dy = (2 x + 2 y ) (II.50a) dt 2 x 2 + y 2 = x 2 + y 2 ( x dt + y dt ) (II.50b) 1 2 2 = r xy + y [( µ y x ) y x ] (II.50c) 1 = ( xy + µy 2 y 4 x 2 y 2 xy ) (II.50d) 1 = ( µy 2 y 4 x 2 y 2 ) (II.50e) y 2 = ( µ x 2 r y ) (II.50f) r 2 sin 2 = [ µ ( x 2 r + y 2 )] (II.50g) = ( µ r ) r sin (II.50h) 1 dy dx = ( x y ) (II.51a) dt r 2 dt dt 1 2 2 2 = r 2 x [( µ y x ) y x ] y (II.51b) 1 2 2 2 2 = r 2 [( µ y x ) xy ( x + y )] (II.51c) 1 2 2 2 = r 2 [( µ r ) r sin cos r ] (II.51d) = ( µ r ) sin cos 1 (II.51e) The constant amplitude solution is r = µ (II.52a) = 1 . (II.52b) It predicts a pure sinusoidal oscillation at all values of µ , that is, one that contains no harmonics unlike the van der Pol model. Moreover the limit cycle is a circle of fixed radius proportional to µ in ( x, y ) phase space. Figures figures 2.6 and 2.7 1 dx dy

PAGE 31

26 show the bifurcation at µ = 0.

PAGE 32

27 Figure II.14: Rayleigh van der Pol model, Stable Fixed Point, X versus t for µ = 0.1. Figure II.15: Rayleigh van der Pol model, Stable Fixed Point, Y versus t for µ = 0.1.

PAGE 33

28 Figure II.16: Rayleigh van der Pol model, Stable Fixed Point, Phase Portrait for µ = 0.1. Given the same initial condition and parameters, the figures shows that the Rayleigh van der Pol model reaches the fixed point as fast as the modified Van der Po l model with a different oscillation trajectory and a more circular phase portrait even for the case of fixed point. Figure II.17: Rayleigh van der Pol model, Limit Cycle, X versus t for µ = 0.1.

PAGE 34

29 Figure II.18: Rayleigh van der Pol model, Limit Cycle, Y versus t for µ = 0.1. Figure II.19: Rayleigh van der Pol model, Limit Cycle, Phase Portrait for µ = 0.1. Given the same initial condition and parameters, the figures shows that the Rayleigh van der Pol model reaches the limit cycle as fast as the modified Van der Pol model with a different oscillation trajectory and a more circular phase portrait. The black circle in the phase portrait indicate the constant amplitude solution where r = µ .

PAGE 35

30 CHAPTER III ESSENTIAL IDEAS FROM STOCHASTIC DIFFERENTIAL EQUATION (SDE) III.I Classical Example of Stochastic Process: Brownian Motion A Brownian motion B t (Wiener process W t ), is a sequence of random variables satisfies the following conditions [25]: 1) B t (0) = 0 with probability 1, 2) for 0 s < t T the random variable given by difference B t B s is normally distributed with mean 0 and variance t s , or B t B s t sN (0 , 1), where N (0 , 1) is normal distribution with zero mean and unit variance. 3) in 0 s < t u < v T , the increments B t B s and B v B u are independent. III.II Generating a Stochastic Process Numerically For the purpose of the numerical simulation, we discretize Brownian Motion as dW t tN (0 , 1). A Sample path of Brownian motion is shown in figure 3.1, it is a random path start at 0 with mean 0 and variance dt.

PAGE 36

31 Figure III.1: sample path of Brownian Motion With the help of measure theory, it has been proved that Brownian Motion like showed in figure 3.1 is nowhere differentiable [2], we can consider it as a special case of fractal to understood its non differen tiablity. III.III Embedding Stochastic Processes in Differential Equations A stochastic process [25] has a general form of dX t = f ( X t , t ) · dt + g ( X t , t ) · dB t (III.1) where B t stand for Brownian motion/Wiener process( W t ), f ( X t , t ) dt is the drift term which carry the deterministic characteristics and g ( X t , t ) B t is the diffusion term which carry the stochastic characteristics. It is worthy noted that because Brownian motion is nowhere differentiable [2], which lead to a result that stochastic process is nowhere differentiable as well, there fore we need write it in difference equation form before we try to solve it analytically and numerically. III.IV Numerical Methods for SDE III.IV.I Euler Maruyama method

PAGE 37

32 t A stochastic differential equation in general in the integral form is given by X ( t ) = X 0 + t f ( X ( s )) ds + 0 t g ( X ( s )) dW ( s ) , 0 t T In general w e use I t o v ersion of i n tegral since most of the meth o ds in st o c hastic di ff ere n tial equation are b ased on I t o c hain rule n u merically and analyticall y . integral solution in differential form: dX ( t ) = f ( X ( t )) dt + g ( X ( t )) dW ( t ) , X (0) = X 0 , 0 t T To solve this equation numerically, we t = T/N where N is the total steps and X i = X ( t i ). The Euler Maruyama (EM) method is similar to method for ODE with additional stochastic term. X i = X i 1 + f ( X i 1 t + g ( X i 1 W, i = 1 , 2 ...N. For the case of additive noise, we write the diffusion terms W tN (0 , 1) (III.2a) = (0 , t ) (III.2b) where is a constant parameter that influence the strength of noise. The level of stochasticity is controlled by this parameter. III.V Solving Stochastic Integrals III.V.I I t o Rule I t rule pr o v e that for st o c hastic pr o cess dX t = f ( X t , t ) dt + g ( X t , t ) dB t , dt 2 = 0 , dtdB t = 0 (III.3a) dB t dt = 0 , dB 2 = dt (III.3b) This is the one of the major breakthrough leading the method to solve the stochas tic integrals. 0

PAGE 38

30 i i i i t t j 2 x x j t j k for i = 1,2. Note that g 2 t t 2 t III.V.II I t o Chain Rule III.V.II.I Single V ariable I t o Chain Rule I t Chain Rule s a ys, for Y t = g ( t, X t ), where X t is st o c hastic pr o cess, w e h a v e 1 2 g 2 dY t = ( t, X t ) dt + ( t, X t ) dX t + ( t, X t )( dX t ) 2 2 (III.4) III.V.II.II Multi v ariate I t o Chain Rule for Y i = g i ( t, X t ), where X t is stochastic process, we have 2 2 2 2 dY = ( t, X ) dt + ( t, X ) dX 1 g + ( t, X ) dX dX (III.5) 2 i x j x k = 0 j, k . III.V.III Analytical examples III.V.III.I Simple Brownian Motion We can see the main difference between Stochastic differential equation (SDE) and ordinary differential equation (ODE) by integrating over Brownian Motion B t : t B s dB s . (III.6) 0 first we guess the form of the solution from deterministic model, in this case, it is y = g ( t, x ) = 1 x 2 , then we try the form for Y t . 1 Y = g ( t, X ) = B 2 (III.7) t t 2 t Using It c hain rule w e h a v e dY = 0 · dt + B · dB 1 2 + · 1 · ( dB ) (III.8) 1 d ( B 2 ) = B dB 1 + dt (III.9) 2 t where ( dB t ) 2 = dt b y It Rule, k j =1 k =1 j j =1 k j =1 k =1 j j =1 t

PAGE 39

31 t t 2

PAGE 40

32 t t 1 B 2 = 2 t 1 B s dB s + 0 2 t 1 ds. (III.10) 0 Which give us the solution, t B s dB s = 0 1 1 B 2 t (III.11) 2 2 III.V.III.II RLC circuit For a simple RLC circuit, we have system where F t = G t + t , 1 LQ t + RQ t + C Q t = F t (III.12) W t : White noise ( W t dt = dB t ), B t : one dimensional Brownian Motion, Q t : the change at time t at a fixed point in an electric circuit, F t : the potential source at time t, G t : the potential source at time t without noise, R : the inductance, L : the resistance, C : the resistance. We rename the variable as following: X = X = X 1 = Q t (III.13) We write equation (3.12) as a system of the first order SED: X 1 = X 2 (III.14a) 1 LX 2 = RX 2 C X 1 + G t + t (III.14b) Now we write the system of SED, equations (3.14a) and (.14b), in the from of difference equation as following Q t X 2 Q t X 2 t

PAGE 41

33 L 1 2 x 2 dX 1 = X 2 dt (III.15a) R 1 1 dX 2 = L X 2 dt CL X 1 dt + L G t dt + L W t dt (III.15b) where W t dt = dB t . In order to solve this two dimensional SDE, we represent the system, (3.15a) and (3.15b), in the matrix form equation, dX = dX t = AXdt + Hdt + KdB t (III.16a) d X = dX 1 , A = 0 1 , H ( t ) = 0 , K = 0 (III.16b) dX 2 1 R 1 G t N o w solving the ma t rix equation (Note that e At = k =0 A K t K k ! converges.), we move AX t dt to the left side of equation (3.16a) and multiply both sides by e At . dX t AX t dt = H ( t ) dt + KdB t (III.17a) e At [ dX AXdt ] = e At [ H ( t ) dt + KdB t ] (III.17b) Now we want to show if we can write the left hand side of equation 3.14b as d ( e At X ). I n attempt to fi gure this out, w e try to apply m ulti v ariate It o c hain rule to g ( t, x , x ) = e At x 1 , dY = d ( e At X ) = Ae At Xdt + e At dX (III.18a) = e At [ dX AXdt ] . (III.18b) (Note that Ae At = e At A .) L CL L L

PAGE 42

34 d ( e At X ) = e At [ H ( s ) dt + KdB t ] . (III.19) Integrate both sides, e At X X 0 = t e As H ( s ) ds + 0 t e As KdB s (III.20) 0 Recall integration by parts we have t f ( s ) dB s = f ( t ) B t 0 t d s ( f ( s )) B s (III.21) 0 Use integration by parts for the last term on the right hand side in equation (3.20), t e As KdB s = e At KB t 0 Substituting (3.22) into (3.20), we get t ( A ) e As KB s ds (III.22) 0 e At X t = X 0 + t e As H ( s ) ds + e At KB t + 0 t Ae As KB s ds (III.23) 0 Which give us the solution X t = e At X 0 + KB t + e At t e As [ H ( s ) + AKB s ] ds (III.24) 0 III.V.III.III Damped Harmonic Oscillator The equation of motion for Damped Harmonic Oscillator with Brownian Noise is given by

PAGE 43

35 dx B t 2 X t kx b dt + a dt = m 2 (III.25) B t : Brownian noise ( W t dt = dB t where W t is white noise), k : spring constant, b : damping constant, a : constant, m : mass. Now let us to write it in the content of stochastic differential equation. 2 X t m 2 + b dX t + kX = a dt t dB t (III.26) dt Equation (3.26) can be written as a system of two first order SDE as following X 1 = X 2 (III.27a) mX 2 + bX 2 + kX 1 = aB t (III.27b) Writing the system of SDE in the from of difference equation, we have dX 1 = X 2 dt (III.28a) b k a dX 2 = m X 2 dt m X 1 dt + m dB t (III.28b) Now we rewrite the system in the form of matrix equation, X = X ( t ) = X 1 (III.29) X 2 X 2

PAGE 44

36 dX 2 k m b m a m dX = AXdt + CdB t (III.30a) d X = dX 1 , A = 0 1 , c = 0 (III.30b) get Moving AXdt to the left side of (3.30a) and multiplying both sides by e At , we dX AXdt = CdB t (III.31a) e At [ dX AXdt ] = e At CdB t (III.31b) We already know from the previous example that we can write the left hand side of equation 3.31b as d ( e At X ). Which give us d ( e At X ) = e At CdB t (III.32) Integrate both sides, e At X t X 0 = t e As CdB s (III.33) 0 Use integrate by parts for the last term in the right hand side in equation 3.33, t e As CdB s = e At CB t 0 Bring the term in we have e At X t = X 0 + e At CB t Which give us the solution t ( A ) e As CB s d s (III.34) 0 t Ae As CB s d s (III.35) 0

PAGE 45

37 4 m 2 m X t = e At X 0 + CB t + e At t Ae As CB s ds (III.36) 0 It is worth noting that the damped harmonic oscillator and RLC circuit in the previous section have similar form of solution. In the RLC circuit, in addition to the noise, the external source is also included in the system. The stochastic RLC in previous section mathematically is similar to the stochastic damped driven harmonic oscillator. Figure 3.2 show the result for the case of damped harmonic motion with weak noise ( = 0 . 01) for weak damping ( b 2 < k ), we can see that for small , the stochastic term did not influence the damping that such, but started to take over after. The phase portrait indicate that the stochastic term keep kicking the system out of the fix point. Figure III.2: 1, m = 5, b = 0.01, k = 1.

PAGE 46

38 2 2 2 CHAPTER IV STOCHASTIC DYNAMICS IV.I Stochastic Van der Pol Model In general, the system always have some level of noise. We take the external additive white noise into consideration and treated it as a forcing term. The stochastic Van der Pol model with additive noise has the form dX = Y · dt + dW 1 ( t ) (IV.1a) dY = ( µ X ) · Y · dt X · dt + dW 2 ( t ) (IV.1b) IV.II Stochastic Rayleigh van der Pol Stochastic Rayleigh van der Pol equation is given by dX = Y · dt + dW 1 ( t ) (IV.2a) dY = ( µ Y X ) · Y · dt X · dt + dW 2 ( t ) (IV.2b) We found out that the analytical method for linear stochastic differential equa tions does not work for multiplicative variables. We cannot rewrite stochastic nonlin ear system into matrix equations when there are multiplicative variables. Moreover, we found no theory about analytical method for stochastic nonlinear systems with

PAGE 47

39 2 2 multiplicative variables. However, we can solve such system using Euler Maruyama method. In the case of circuit, we set it up in a more realistic way such that dW 1 ( t ) = 0. dX = Y · dt (IV.3a) dY = ( µ Y X ) · Y · dt X · dt + dW ( t ) (IV.3b) In figures 4.1 to 4.6 the time sires and phase portraits are shown for µ = 0 . 1 and µ = 0 . 1, below and above deterministic bifurcation value µ = 0, for = 0 . 01. Figure IV.1: Rayleigh van der Pol model, Stable Fixed Point, X versus t for µ = 0.01.

PAGE 48

40 Figure IV.2: Rayleigh van der Pol model, Stable Fixed Point, Y versus t for µ = Figure IV.3: Rayleigh van der Pol model, Stable Fixed Point, Phase Portrait for µ = Figure 4.1 to 4.3 show that in the case of Stable Fixed Point, the figures for X and Y versus t suggest that the oscillation damps away and the stochasticity takes over. The phase portrait is now a fuzzy circular trajectory: the trajectory fluctuate while reaching the stable fixed point. Right after the system reach the stable fixed point, the stochasticity/noise kick the system out of it and that cause the system to restart damping.

PAGE 49

41 Figure IV.4: Stochastic Rayleigh van der Pol model, Limit Cycle, X versus t for µ = 0.01. Figure IV.5: Stochastic Rayleigh van der Pol model, Limit Cycle, Y versus t for µ = 0.01.

PAGE 50

42 Figure IV.6: Stochastic Rayleigh van der Pol model, Limit Cycle, Phase Portrait for Figure 4.4 to 4.6 showed that in the case of limited cycle, the figures for X and Y versus t suggest that the model have a fuzzy limit cycle instead, the fuzzy trajectory and widened cycling region of the phase portrait agree with it. The black circle in the phase portrait indicate the constant amplitude solution where r = µ .

PAGE 51

43 CHAPTER V BEHAVIOR NEAR THE BIFURCATION V.I Numerical Simulation In this chapter the numerical Euler Maruyama method is used to solved the stochastic Rayleigh van der Pol equation. V.I.I P oinca r ´ e Section Study and statistical analysis of a stochastic system provides useful tools for un derstanding the behavior and propertie s of the system. In order to obtain and measure the statistical behavior of the nonlinear stochastic system, we chose to measure the statistical pro p erties of the P oinca r ´ e section. Numerical si m ulation of P oinca r ´ e section for deterministic system requires Jacob nian which need take the (partial) derivatives of the system. That bring the difficulty since stochastic differential equation is nowhere differentiable. In o r der to obtain the P oinca r ´ e section for o ur system, w e set up a (partial) h y p er Every time the trajectory in phase plane cross the hyper plane, we mark it. This can give a good We only need to collect data for any one of the axis(i.e. the positive x axis) that the trajectory cross thanks to the circular phase portrait of Rayleigh van der Pol.

PAGE 52

44 Figure V.1: Poin c a r ´ e S e ction for Phase Port r ait of St o chastic R ayleigh van der Pol model under different noise level, µ = 0 . 1 = 0 . 01 , t > 1300 , for t [0 , 2000] . As figure 5.1 shows, under the influence of relatively weak noise ( = 0 . 01), the system starting became slightly unstable. The stochastic term causes the fluctuation of the limit cycle, which showed as the fuzzy path and the widened circular region of the phase portrait. As the strength (the level) of noise gets larger the fluctuation also increases.

PAGE 53

45 Figure V.2: Poin c a r ´ e S e ction for Phase Port r ait of St o chastic R ayleigh van der Pol model under different noise level, µ = 0 . 1 1300 , for t [0 , 2000] . Under the influence of relatively stronger noise ( = 0 . 05), figure 5.2 suggests that the system became way more unstable. We can also tell by the increasing of the width of P oinca r ´ e section.

PAGE 54

46 Figure V.3: Hist o g r am of Poin c a r ´ e S e ction for Phase Port r ait of St o chastic R aylei gh van der Pol model under different noise level, µ = 0 . 1 = 0.01. W e obtained the histogram of P oinca r ´ e section to obser v e the statistical pro p er t y of the system. It is a equi v ale n t to vi e w the cross section of th e P oinca r ´ e section. Figure 5.3 showed that under the influence of relatively weak noise ( = 0 . 01), we can see an almost normal distribution. The red line indicate the constant amplitude solution where r = µ .

PAGE 55

47 Figure V.4: Hist o g r am of Poin c a r ´ e S e ction for Phase Port r ait of St o chastic R aylei gh van der Pol model under different noise level, µ = 0 . 1 = 0.05. Figure 5.4 showed that in the case of relatively stronger noise ( = 0 . 05), the distribution had increased the bell width and shift to the left while be less normal distributed. V.I.II Behavior of Mean and Variance with Bifurcation Parameter for Fixed Noise Level We found that the variance does not influence by dt for our system with mu / = 0 even though for a stochastic process [ V ar ] = c 2 · dt where c 2 = for our case. After the data collecting, w e computed the mean and v a r iance of the P oinca r ´ e section versus µ to see if we can get any information for the bifurcation, as shown from figure II.12 to figure II.17. As we can tell from the figures, for small noise level (i.e. = 0 . 01), we can easily tell where the bifurcation happens (around µ = 0). From the sudden increase of mean in mean versus µ and the peak in variance versus µ , at

PAGE 56

48 this noise level we can also observe the phase portrait for indications of bifurcation by changing µ to see the difference in phase portraits for limit cycle and fixed point. But for raised noise level (i.e. = 0 . 05), it is difficult to find out where the bifurcation happens by the phase portrait, but the change of mean in mean versus µ and the peak in variance versus µ . For intensified noise level ( 0 . 5), the mean versus µ cannot give us much indication about where the bifurcation happens, however the variance versus µ can still give us a region where the bifurcation happens. The theory of Hopf bifurcation defines where bifurcation happens as a critical point. However, after observing the experimental data, it seems like for the case of bifurcation under the influence o a higher amplitude while shifting to the right instead of stay at µ = 0 as the stochastic term increases. W e h a v e o btained the mean and v ariance of the P oinca r ´ e section v ersus µ . It is worth noting that while we ran the simulation, we seeded the pseudo random number generator to ensure we were studying the same sequence of random variables.

PAGE 57

49 (a) Mean v.s. µ for = 0.01. (b) Variance v.s. µ for = 0.01. Figure V.5: Stochastic Rayleigh van der Pol model, behavior of mean and variance, 1300 .

PAGE 58

50 For = 0 . 01, we can tell from figure 5.5 where the bifurcation happened by the sudden increase in mean versus µ as well as the peak in variance versus µ . The peak of variance is around µ = 0 where the bifurcation happens for the deterministic model.

PAGE 59

51 (a) Mean v.s. µ for = 0.02. (b) Variance v.s. µ for = 0.02. Figure V.6: Stochastic Rayleigh van der Pol model, behavior of mean and variance, 1300 .

PAGE 60

52 For = 0 . 02, we are still be able to tell from figure 5.6 where the bifurcation happened by the sudden increase in mean versus µ as well as the peak in variance versus µ µ is a little more smooth than in the care of = 0 . 01. The variance also had a higher value at the peak and wider peak region. The peak of variance slightly shifts to the right side of µ = 0.

PAGE 61

53 (a) Mean v.s. µ for = 0.05. (b) Variance v.s. µ for = 0.05. Figure V.7: Stochastic Rayleigh van der Pol model, behavior of mean and variance, 1300 .

PAGE 62

54 In figure 5.7, for = 0 . 05, we are still be able to tell where the bifurcation happened by the sudden increase in mean versus µ as well as the peak in variance versus µ . However, µ is more smooth than in the care of = 0 . 2. The variance also had a higher value at the peak and wider peak region compare with the case for = 0 . 02. The peak of variance shifts to the right side of µ = 0 even more.

PAGE 63

55 (a) Mean v.s. µ for = 0.1. (b) Variance v.s. µ for = 0.1. Figure V.8: Stochastic Rayleigh van der Pol model, behavior of mean and variance, 1300 . For = 0 . 1, as showed in figure 5.8, we had difficulty to tell where the bifurca

PAGE 64

56 tion happened by mean versus µ , but we can tell the region where the bifurcation happened. We can still be able to tell where the bifurcation happened by the peak in variance versus µ . The variance also had a higher value at the peak and wider peak region compare with the previous cases. The peak of variance shifts to the right side of µ = 0 more significantly.

PAGE 65

57 (a) Mean v.s. µ for = 0.7. (b) Variance v.s. µ for = 0.7. Figure V.9: Stochastic Rayleigh van der Pol model, behavior of mean and variance, 1300 .

PAGE 66

58 Figure 5.9 shows mean and variance versus µ for significantly bigger noise/s tochastic term ( = 0 . 7). We cannot tell where bifurcation region is just by mean versus µ . The peak region of the variance, which indicates the bifurcation region, con tinued to widen with increased peak value compare with the previous cases. There seems to be a peak in the bifurcation region at µ = 0 . 7 as well. The surface plots in appendix A shows a better visualization for coupled system.

PAGE 67

59 (a) Mean v.s. µ for = 0.8. (b) Variance v.s. µ for = 0.8. Figure V.10: Stochastic Rayleigh van der Pol model, behavior of mean and variance, 1300 .

PAGE 68

60 For intensified noise/stochastic term ( = 0 . 8) showed in figure 5.10. We are still unable to tell where bifurcation region is just by mean versus µ . The bifurcation region can still be found by the peak region of the variance, which is continued to widen with increased peak value like in previous case. It is more clear to see the µ = 1 here.

PAGE 69

61 CHAPTER VI CONCLUSION AND FUTURE WORK For the system of (modified) Rayleigh Van der Pol with additive stochastic term, the the right instead of staying at µ = 0 as the stochastic term increases. VI.I Future Work One of the applications of Rayleigh Van der Pol Oscillator is in the simulation of the Wein Bridge circuit[30]. VI.I.I Coupled System One of the applications of the coupled nonlinear system is in the simulation of the animal gaits[3][4][27]. VI.I.I.I Two Synaptic Coupled Stochastic Rayleigh van der Pol Model In the content of neuroscience, the neuron network can be represented in the form of coupled oscillators. There are different kind of coupling method such as synaptic coupling and diffusive coupling. We choose synaptic coupling since it is a more realisti c way to represents how different neurons connect each other.

PAGE 70

62 2 2 1 1 2 2 2 2 System of Two Synaptic Coupled Stochastic Rayleigh van der Pol Equations: dX 1 = Y 1 · dt + 1 · X 2 · dt + dW 1 ( t ) , (VI.1a) dY 1 = ( µ Y X ) · Y 1 · dt X 1 · dt + 2 · Y 2 · dt + dW 2 ( t ) , (VI.1b) dX 2 = Y 2 · dt + 3 · X 1 · dt + dW 3 ( t ) , (VI.1c) dY 2 = ( µ Y X ) · Y 2 · dt X 2 · dt + 4 · Y 1 · dt + dW 4 ( t ) , (VI.1d) Where dW 1 ( t ) = dW 3 ( t ) = 0, i = i = 1 , 2 , 3 , 4. See appendix A for analysis and simulation for two coupled Rayleigh van der Pol system .

PAGE 71

63 REFERENCES [1] Anishchenko Vadim S., Astakhov Vladimir, Neiman Alexander, Vadivasova Tatjana., Schimansky Geier Lutz, Nonlinear Dynamics of Chaotic and Stochastic Systems: Tutorial and Modern Developments , Springer, 2007. [2] Patrick Billingsley, Probability and Measure , Wiley Series in Probability and Statistics, 2012. [3] Pietro Luciano Buono, Martin Golubitsky , Models of central pattern generators for quadruped locomotion I. Primary gaits , J. Math. Biol, 2001. [4] Pietro Luciano Buono, Models of central pattern generators for quadruped locomotion II. Secondary gaits , J. Math. Biol, 2001. [5] Kevin Burrage, Ian Lenane, Grant Lythe, Numerical Methods for Second Order Stochastic Differential Equations , SIAM J. SCI. COMPUT Vol. 29, No. 1, pp. 245264, 2007. [6] T. K. Caughey, Response of Van , Journal of Applied Mechanics, 1959. [7] T. K. Caughey, On the Respones of a Class of Self excited Oscillators to Stochastic Excitation , Int. J. Non Linear Mechanics, 1967. [8] T. K. Caughey, Nonlinear Theory of Random Vibration , Advances in Applied Mechanics Volume 11, 1971. [9] Andrey Dmitrieva, Victor Dmitrieva, Oleg Sagaydakb, Olga Tsukanovaa, The Application of Stochastic Bifurcation Theory to the Early Detection of Economic Bubbles , Procedia Computer Science, 2017. [10] Thai Son Doan, Maximilian Engel, Jeroen S W Lamb, Martin Ras mussen, Hopf bifurcation with additive noise , Nonlinearity 31 4567, 2018. [11] Armin Fuchs, Nonlinear Dynamics in Complex Systems, Theory and Applica tions for the Life , Neuro and Natural Sciences , Springer, 2013. [12] M.A. Garstens, Noise in Non Linear Oscillators , J. Appl. Phys., 1957. [13] Kiyosi It o , St o chastic Int e g r al , Pr o c . Imp. Acad., 1 944. [14] Eugene M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting , The MIT Press, 2007.

PAGE 72

63 [15] Kurt Jacobs, Stochastic Processes for Physicists: Understanding Noisy Sys tems , Cambridge University Press, 2010 [16] Samuel Karlin, Howard M. Taylor, A First Course in Stochastic Pro cesses , Academic Press, 1975. [17] Tatsuyuki Kawakubo, Shigeharu Kabashima, Stochastic Process in Self Excited Oscillation , Journal of Physical Society Japan, 1974. [18] Peter Eris Kloeden, Eckhard Platen, Henri Schurz, Numerical Solu tion of SDE Through Computer Experiments , Springer, 2003. [19] Christian Kuehn, A Mathematical Framework for Critical Transitions: Bi furcations, Fast Slow Systems and Stochastic Dynamics , Physica D: Nonlinear Phenomena, Vol. 240, No. 12, pp. 1020 1035, 2011. [20] H.K. Leung, Stochastic Transient of a noisy Van der Pol Oscillator , Physica A, 1995. [21] H.K. Leung, Stochastic Hopf Bifurcation in a biased Van der Pol Model , Phys ica A, 1998. [22] Carlo Liang, Gabriel J. Lord, Stochastic Methods in Neuroscience , Oxford University Press, 2010. [23] Marc Mendler, Johannes Falk, Barbara D rosse , Analysis of stochastic bifurcations with phase portraits , arXiv, 2018. [24] John F. Moxnes, Kjell Hausken, Introducing Randomness into First Order and Second Order Deterministic Differential Equations , Advances in Mathemat ical Physics Volume 2010, 2010. [25] Bernt Oksendal, Stochastic Differential Equations (6th Ed.): An Introduction with Applications , Springer, 2010. [26] Sheldon M. Ross, Stochastic Processes , John Wiley and Sons, Inc., 1996. [27] Mozhdeh Saffari Parizi, Coupled Biological Oscillators and Animals Gaits , MIS Thesis, University of Colorado Denver, 2004. [28] P T. D. Spanos, Stochastic Analysis of Oscillators with Non Linear Damping , Int. J. Non Linear Mechanics, 1978. [29] Steven H. Strogatz, Nonlinear Dynamics and Chaos (2nd Ed.): With Appli cations to Physics, Biology, Chemistry, and Engineering , Westview Press, 2014. [30] Randall Tagg, UC Denver Physics Department, Private Communication. [31] Ruoshi Yuan, Xinan Wang, Yian Ma, Bo Yuan, Ping Ao , Explor ing a noisy van der Pol type oscillator with a stochastic approach , PHYSICAL REVIEW E 87, 062109, 2013.

PAGE 73

64 APPENDIX A ANALYSIS AND SIMULATION FOR COUPLED SYSTEM This appendix contains the analysis and simulation for two coupled Rayleigh van der Pol system. For the coupled case, we are studying coupling strength v.s. µ for different noise level , find where the bifurcation happens; mean and v a riance of the P oinca r ´ e section, fi nd where the bifurcation hap p ens; surface plot of mean and variance of each oscillator v.s. and µ for fixed noise level (i.e. 2 , 4 = 0 . 05). The results for the first oscillator in the system shows in the figures below.

PAGE 74

65 1 Figure A.1: Coupled System of Stochastic Rayleigh 2 , 4 = 0 . For the deterministic twp coupled Rayleigh van der Pol equations, we can see from the Figure A.1 that the bifurcation happened at = 2 µ which is one of the eigenvalues for the system.

PAGE 75

66 2 2 Figure A.2: Coupled System of Stochastic Rayleigh van 2 , 4 = 0 . 01 , the gray points indicate the bifurcation region. The intriguing part happened once we introduce stochastic term into the coupled system as showed in figure A.2. Even for = 0 . 01, the bifurcation transit from a line The bifurcation region is a combination of positive part of = 1 µ and negative part of = 1 µ both shift to the right of axis. Moreover, the bifurcation region widened more for positive .

PAGE 76

67 1 1 2 2 Figure A.3: Coupled System of Stochastic Rayleigh 2 , 4 = 0 . 05 . For = 0 . 05 which showed in figure A.3, the bifurcation region is still a combina tion of positive part of = µ and negative part of = µ . The bifurcation region keep shifting to the right of axis while widened more for positive than negative .

PAGE 77

68 1 1 2 2 Figure A.4: Coupled System of Stochastic Rayleigh 2 , 4 = 0 . 1 . For = 0 . 05 which showed in figure A.4, the bifurcation region is still a combina tion of positive part of = µ and negative part of = µ . The bifurcation region keep shifting to the right of axis while widened more for positive than negative .

PAGE 78

69 Figure A.5: Coupled System of Stochastic Rayleigh van der Pol, mean of the first oscillator v.s. µ for different 2 , 4 = 0 . 05 . Now we fixed the stochastic term(i.e. = 0 . 05), we want to see how the statistical properties behave for different coupling strength . Figure A.5 suggests that only the magnitude of will influence the mean of the system. Increasing the magnitude of , has the opposite effect as increasing . The increase of the magnitude of value of mean bigger while shifting the bifurcation point/region to the left of the axis.

PAGE 79

70 Figure A.6: Coupled System of Sto chastic Rayleigh van der Pol, variance of the first oscillator v.s. µ for different 2 , 4 = 0 . 05 . Meanwhile, Figure A.6 suggests similar tendency, that only the magnitude of will influence the mean of the system. Increasing the magnitude of , has the opposite effect as increasing . The increase of the magnitude of make the decrease the amplitude of the peak region of bifurcation region while narrowing and shifting the bifu rcation region to the left of the axis.

PAGE 80

71 Figure A.7: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the 2 , 4 = 0 . 05 . Figure A.8: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the 2 , 4 = 0 . 05 . Figure A.7 and A.8 shows the surface plots for mean, µ and of the first oscillator at different angles. The surface plots give a better visualization about mean, µ , and bifurcation at fixed stochastic terms ( 2 , 4 = 0 . 05): mean increases as the magnitude of ion region is.

PAGE 81

72 Figure A.9: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the 2 , 4 = 0 . 05 . Figure A.10: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the 2 , 4 = 0 . 05 . Figure A.9 and A.10 shows the surface plots for variance, µ and of the first oscil lator at different angles. The surface plots give a better visualization about variance, and bifurcation at fixed stochastic terms ( 2 , 4 = 0 . 05): variance decreases as the

PAGE 82

73 magnitude of increase, the peaks indicate the bifurcation region. There seems to Figure A.11: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the second oscillator for mean, µ and 2 , 4 = 0 . 05 . Figure A.12: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the second oscillator for mean, 2 , 4 = 0 . 05 . Figure A.11 and A.12 shows the surface plots for mean, µ and of the second oscillator at different angles. The surface plots give a better visualization about

PAGE 83

74 mean, and bifurcation at fixed stochastic terms ( 2 , 4 = 0 . 05): mean increases as the magnitude of bifurcation region is. Figure A.13: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the second oscillator for variance, µ and 2 , 4 = 0 . 05 . Figure A.14: Coupled System of Stochastic Rayleigh van der Pol, surface plot of the second 2 , 4 = 0 . 05 . Figure A.13 and A.14 shows the surface plots for variance, µ and of the second

PAGE 84

75 oscillator at different angles. The surface plots give a better visualization about variance and bifurcation at fixed stochastic terms ( 2 , 4 = 0 . 05): variance decreases as the magnitude of increase, the peaks indicate the bifurcation region. There seems both bifurcation regions are a little different compare to the first osc illator. The coupled system behaves a little different than we anticipated, the nature of the coupled stochastic nonlinear system and its behavior is not well understand yet. As we can tell from the plots, the bifurcation region has an excepted shift once the randomness has been introduced to the system, then the bifurcation region widens and shift to the right as the stochasticity increases. For the fixed , both positive and negative coupling shift the bifurcation region to the left, the bifurcation narrowed as the absolute value of increases.