The extended Malkus-Robbins dynamo [Moroz, 2003] reduces to the Lorenz equations when one of the key parameters, $\beta$, vanishes. In a recent study [Moroz, 2004] investigated what happened to the lowest order unstable periodic orbits of the Lorenz limit as $\beta$ was increased to the end of the chaotic regime, using the classic Lorenz parameter values of r = 28; $\sigma$ = 10 and b = 8=3. In this paper we return to the parameter choices of [Moroz, 2003], reporting on two of the cases discussed therein

Moroz, I. M.

English

Alarcón, E.

Repositorio da Universidade da Coruña

Métodos numéricos E. Alarcón ETSII. Universidad Politécnica de Madrid l. EL METODO DE LOS ELEMENTOS FINITOS (MEF) Sea un domino elástico O, de contorno an, sometido a unas fuerzas por unidad de volumen ~ y tensiones X = T en an de modo que - >:o:::- = an= ao u ao,. -Si las tensiones ?' en Q están en equilibrio con ~ y X y si además ~ y ~ son los desplazarientos y correspondientes defo~acio~es del cuerpo es sabido que es posible escribir el principio de los trabajos virtuales (PTV) como: Si la variación de los desplazamientos 8 ~ respeta las condiciones de contorno en aou debe ser nula allí y la última integral sólo se extiende a ao donde, como se ha dicho, se tiene X = T t ~ = la ecuación (1) puede escribirse: In 8 E: ·<!-In 8 ~T · ~- Ian 8 uT • T = O (1.2) 244 E. ALARCÓN El método de Rayleigh-Ritz consiste en aproximar el campo de des-plazamientos en la forma: u(x y z)"" N(x y z)ue (1.3) donde ue son valores incógnita en una serie de puntos seleccionados y N(x Y z) funciones de interpolación. Con ello se cumple E (x y z)=B(x y z)ue (1.4) donde B se obtiene de N utilizando la definición de las deformaciones ;::,; :::::; 1 E= (u.+u .. ) l] 2 l,j j,l (1.5) Las variaciones de ~ y ~ sólo afectan a ~e puesto que las "! son funciones fijas. Así o u= N oue OE=Boue La sustitución de (1.6) en (1.2) conduce a Puesto que las variaciones son arbitrarias ello implica f~t<!= f~T ~+ f~T·! Q - - Q - - ()Qf - -(1.6) (1.7) ( 1.8 ) MÉTODOS NUMÉRICOS 245 Por otro lado supóngase que existe un estado inicial de deformaciones :o y que el comportamiento sigue la ley de Hooke ( 1.9) donde J? es una matriz que contiene el modulo de Young E y el coeficiente de Poisson v en el caso de cuerpo isótropo o los correspondientes juegos de E y v para los casos anisótropos. La sustitución y reordenación de ( 9 ) y ( 4 ) en ( 6 ) conduce a (1.10) Obsérvese que todos los términos del segundo miembro son conoci-dos, por lo que (1.10 ) es un sistema lineal de ecuaciones cuya resolución permite obtener u~e . El método de los elementos finitos ( MEF ) sigue esta forma de proceder y su aportación se refiere a la elección de las funciones de interpolación. Para ello se subdivide el dominio en elementos ge definidos de forma convencional mediante nudos (vg. :triángulos con tres nudos, cua-driláteros con cuatro nudos, triángulos curvos de seis nudos, cuadriláteros de lados curvos de ocho nudos, tetraedros de cuatro nudos, etc. ) de forma que su unión recree el dominio total. A continuación se definen las funciones de forma que sean unidad en cada nudo y nulas en el resto. Por ejemplo, la figura 1 recoge la subdivisión de un dominio plano en triángulos y la forma de la función de interpolación correspondiente a uno de los nudos. Con ello se consigue que los coeficientes u~e sean precisamente los desplazamientos nodales que se van buscando, que bajo ciertas condiciones 246 E. ALARCÓN Figura 1 de continuidad dependientes del problema, las integrales pueden construir-se como suma de integrales y que la matriz BT D B , además de ser si-métrica por el propio método de interpolación, tenga estructura en banda lo que presenta ventajas tanto de tipo computacional como de almacenamiento. El proceso se organiza entonces de forma repetitiva: para cada ele-mento se obtiene la matriz de rigidez k e y el vector de cargas fe. :::::: ::::: k_e = f Jt J!J! - - - -n, fe= J~T ~+ J~T·J.:+ f-'!T I}eo ~ n, an, n, ~ (1.11) Estos cálculos se realizan en el sistema de ejes «natural» del elemento (figura 2) que se denominan «ejes locales». MÉTODOS NUMÉRICOS 247 X Figura2 ;-+ + -' +- ..¡.. Figura 3 La mayor parte de las veces se realiza una aplicación sobre un rec-tángulo unidad ( figura 3 ) para utilizar un método de integración numérica ( vg. : regla de Gauss ). A continuación se procede a realizar una rotación a ejes globales lo que permite un suma booleana de matrices y vectores ~ = L k~e e ~= Lfe ( 1.12) 248 E. ALARCÓN y de esta forma se tiene el sistema Ku=F (1.13) cuya resolución permite obtener los desplazamientos ~ a partir de los cua-les se calculan los esfuerzos y tensiones deshaciendo, elemento a elemento, las transformaciones anteriores. El esquema general de un programa de elementos finitos en elastici-dad es el siguiente: 1 PROGRAMA PRINCIPAL 1 1 Introducción de datos 1 Para cada elemento - Cálculo de ke - Cálculo de f e - Rotación a coordenadas globales - Colocación en matrices globales 1$. y !: 1 Aplicación de las condiciones de contorno 1 Resolución del sistema de ecuaciones 1 Cálculo de tensiones y esfuerzos 1 Presentación de resultados MÉTODOS NUMÉRICOS 249 2. VISCOELASTICIDAD En este apartado se recuerdan dos alternativas que son de interés en el tratamiento numérico de problemas viscoelásticos. Por un lado se resu-me la formulación de la ley de comportamiento mediante variables de esta-do y a continuación se presenta una formulación de tipo potencial. Final-mente se muestra un ejemplo sencillo. 2.1. VARIABLES DE ESTADO Si se utiliza la representación integral de un cuerpo viscoelástico para el caso monodimensional t e(t) = J D(t, r )da (2.1) 'ro donde D es la función especifica de fluencia. La integración de ( 2.1 ) por partes conduce a e(t)= a((t))+ r d(t,r)a(r)dr E t 'ro a d(t, r) =- ar D(t, r) 1 E( t) = D( t' t) (2.2) Si se trata de un material sin envejecimiento a(t) Jt e(t) =- - + d(t- r)a(r)dr E 'ro (2.3) 250 E. ALARCÓN En la mayoría de los casos prácticos se tiene la posibilidad de aproxi-mar d(t- r) mediante una serie de Dirichlet-Prony n 1 _1-T d(t-r) = L e 8; i=l 1]¡ ( 2.4) donde los valores 1J; ; 8; se obtienen mediante identificación de parámetros tras una serie de ensayos de laboratorio. Si se llama se tiene y, derivando (2.5) 1-T r~ 1 --q¡(t) = J, -e 8; a( r)dr '<o 1]¡ a(t) ~ e(t) = -+ ~q¡(t) E i=l • q¡(t) a(t) q.(t) + --··· = -----1 8¡ 1]¡ (2.5) (2.6) (2.7) Las ecuaciones (2.6) y (2.7) junto con las condiciones iniciales (2.8) forman un sistema equivalente a (2.3). Se ha pasado así de una representación integral a un sistema de n ecuaciones desacopladas de primer orden. MÉTODOS NUMÉRICOS 251 Puesto que, tan pronto como se conocen los q¡ (t) es posible calcular E (t), se las suele denominar variables de estado. En el caso lineal se puede proceder por intervalos de tiempo de longitud !1t con lo que (t+LH-T) ( ) t+At-T ( ) 1 --- (5 r r+t.r --- (5 r q(t+l1t)= re e, --dr+J e 7),(r) --dr=l +/,(2.9) ' Jo 17¡ ( r) t 17¡ ( r) t -pero f>t t-T ( ) --J~ -- (5 r 1 = e e, e e, -- dr = q (t )e e, t o 17¡ ( r) ' C.t (2.10) (5( r) Para obtener 12 se puede suponer que 11. ( r) permanece constante en el intervalo !1t con lo que ' (5 t ------¡- Jt+f>t 7i (5 t -~ 1 = · e ' e' dr = --(). 1- e ' ( ) t+f>t T ( ) ( f>t J 2 17¡ ( t) t 11/ t) l (2.11) Con ello se tiene el esquema incremental (2.12) La ecuación (2.6) se puede interpretar como (2.13) 252 E. ALARCÓN lo que permite obtener unas tensiones (2.14) La comparación de (2.14) con (1.9) permite interpretar ev como una deformación inicial a efectos de_ cálculo con elementos finitos y ello conduce al siguiente sistema operativo : a) Se realiza un análisis lineal para t =O con las propiedades elásticas de la estructura, lo que permite calcular la tensión en cada pieza. b) Suponiendo que en el intervalo t = O t = 11t se mantienen constantes las tensiones se calcula ev mediante (2.12) y (2.13). e) Se repite el cálculo elástico añadiendo la carga equivalente a la deform~ción visc~sa que, según (1.10) es J:/t l} (. Ello permite obte-ner tensiOnes medtante (2.14) - - -El proceso se repite hasta finalizar el tiempo prescrito. 2.1. EJERCICIO Determinar la matriz de rigidez de una barra a tracción de longitud L, área A y módulo de elasticidad E. Usando el formalismo de elementos finitos. X MÉTODOS NUMÉRICOS 253 El elemento se muestra en la figura junto con los ejes locales y globales. Se suponen fuerzas de volumen ~ = (b,O) y un estado de deformación inicial t:0 • Para simplificar se admite que tanto b como t:0 son constantes dentro del elemento. Tomando una función de interpolación lineal se tiene : o X 1--L X L o La relación entre desplazamientos y deformaciones es : La ley de comportamiento es : a= E(t:-t:0 ) Con ello 1 o -1 o . l T AE o o o o K'= B EBdv= · V L -1 o 1 o o o o o 254 E. ALARCÓN y 1 '-1 o 0 1 +AEe0 1 o 0 Puesto que K y P están referidas a ejes locales xL, yL es preciso realizar una rotación R con lo que : donde la matriz de rotación es : e S o o -S e o o R= ; e = cosa; s = sen a o o e S o o -s e Con ello en coordenadas globales se tiene : e2 2 se -e -se 2 2 EA S -e -s K= L e2 se s2 MÉTODOS NUMÉRICOS 255 e -e L S -S P = Jt + Fa = Ab 2 + AEE0 e e S S Una vez calculados los desplazamientos se obtienen los esfuerzos axiles mediante: Sea por ejemplo la estructura de la figura: Fy Si se aplican los resultados anteriores a cada una de las barras se tiene: EA K-~-"- 2J2i -1 -1 -1 EA ·K-~­-1 ' 22 - 2J2i o o o o -2.fi o o 2.fi 1 1 -1 -1 EA 1 -1 -1 ·K-~-' 33- 2J2i 256 Con lo que el sistema global de ecuaciones es : F_d -1 -1 o o -1 Fv~ 1 o o 1 F'.rz o o o o F~.z 1 o -1 EA ... ... ... Fx3 =üu 2 o F,., 2(1 +..Ji) Si se aplican las condiciones de contorno : se tiene el sistema : con lo que u1=v1=0 u2=v2=0 u4=v4=0 2L V - F 3 - EA(2 + J2) y 2 o o o o .. . o -1 o u, =0 o v, =0 o u, =O o v, =0 ... ... o u, -1 v, 1 u, =o 1 v, =o E. ALARCÓN MÉTODOS NUMÉRICOS 257 y para los esfuerzos 2 S = F 2 2 + J2 y2 En el caso viscoelástico -1 o 1 o EA{ S = L ( u2 - u1 ) cos a + ( v 2 - v 1 ) sen a - Ev L} Supóngase que la barra 2 está formada por un material standard con propiedades E1 = 100; E2 = 100; 1J = 1000 mientras que las barras 1 y 2 son elásticas con modulo E1 = 1 OO. Se toma Fx3 = O, Fy 3 = 1 Se comienza con un calculo elástico en t = O V,= -0.59 o S1 = S3 = o.29 S2 = o.59 258 E. ALARCÓN Suponiendo que los esfuerzos en las barras permanecen constantes en el intervalo t =0 y t = l se calcula la deformación viscoplástica en la barra 2 lo que conduce a: t: = 059 (1- e-o.1) = 0.00056 V 100 Ahora se hace un nuevo cálculo elástico en t = l tomando t:v como una deformación inicial. La carga nodal correspondiente es ( 1)T Pv = 0.056 O -1 O Si se incluye esta nueva carga en el cálculo se tiene: y por tanto Vo=-0.62 -' S1 = S3 = o.31 S2 = o.s6 El proceso se repite para t = 2, 3, etc. y se obtienen los resultados de la figura V 0.5 -------------------sj.S¡=$3 0.1 5 10 MÉTODOS NUMÉRICOS 259 2.2. PoTENCIAL DE DISIPACIÓN Como es sabido en elasticidad se definen los potenciales de energía elástica u(~) y energía complementaria w( ~) mediante aw t:=-= aa ( ) au aw d a:t: =adt:+t:da= dt:+ -da=d(U+W) = = = = = = at: = aa = es decir (2.15) que representa la llamada transformación de Legendre respecto a la que los dos potenciales son conjugados. La existencia de estos potenciales se sigue del concepto de reversibilidad y de función de estado termodinámico. En viscoelásticidad no existe una motivación semejante pero es costumbre utilizar una formula-ción similar con potenciales au(~) a= -. at: (2.16) 260 E. ALARCÓN Además se suele admitir que los potenciales son homogéneos. Por ejemplo si U es homogéneo de grado m y W homogéneo de grado k existe entre ellos la relación m k= m-1 En efecto, en virtud de la definición de potenciales aw cr: E = cr: - -= k W ~ ~ ~ acr y por tanto u(~) = kW- w = (k - l)w( ~) Así Pero debido a la homogeneidad lo que implica k= (k -l)m MÉTODOS NUMÉRICOS El producto es la función de disipación y puede verse que 1 dD 1 dD E= --= k aa 261 (2.17) lo que justifica el conocido resultado de la función de disipación como po-tencial viscoso (2.18) El conocimiento de W, U, o D permite determinar la ley de compor-tamiento (2.18) En general se establece una tensión equivalente cr y una deforma-• e ción equivalente E tales que la disipación coincida con la de ensayo monodimensional e Con ello el control de la deformación viscosa se reduce al caso monodimensional y se puede repetir el procedimiento allí indicado. Como ejemplo se va a tratar el caso en que que corresponde a una familia muy popular de leyes. 262 E. ALARCÓN El caso más sencillo es n = 1 con lo que Aplicando las definiciones anteriores se tiene y generalizando el resultado • 3 E¡ = 2(jz ( <1¡ - aocr) o (2.20) que puede escribirse como (2.21) Puede observarse que en este caso se tiene una deformación volumétrica nula MÉTODOS NUMÉRICOS 263 Por otro lado usando las relaciones (2.20) es posible obtener y despejando las tensiones (2.22) Ello permite calcular la disipación sustituyendo (2.22) en (2.19) lo que conduce a Obsérvese que en el caso monodimensional se tiene: 3 (je e =- --2 2a 2 3 o (je 1 =- e =e 2a2 2 t 3 o (2.23) 264 E.ALARCÓN En particular ello implica la ley de comportamiento lineal que da el significado de 0"0 : cuya sustitución en (2.21) conduce a la ley (2.24) En el caso monodimensional se tiene D ~ <J1 ¡1 ~u,[::]= 2W de modo que utilizando (2.19) se llega a (2.25) y también según (2.23) MÉTODOS NUMÉRICOS 265 Así pues • (2.26) Una expresión más familiar se obtiene observando que, debido a la condición de incompresibilidad se tiene y por tanto e = e (2.27) El proceso operativo en el caso monodimensional consiste en las si-guientes etapas: a) En el instante t = O se calcula la solución elástica del problema lo que permite calcular desplazamientos y tensiones. b) Se supone que las tensiones anteriores permanecen constantes en un intervalo /j.t lo que permite calcular la tensión equivalente a e. Con ella se obtiene, por el procedimiento unidimensional, la correspondiente defor-• mación viscosa e v así como e • • 3 • S .. ev =e .. = -ee ___'J_ lj 2 V (J e e) Con las deformaciones anteriores se calculan las cargas equivalentes 266 E. ALARCÓN y se resuelve de nuevo el problema para obtener nuevas tensiones El proceso se repite hasta agotar el tiempo de cálculo. Evidentemente el proceso se puede mejorar iterando dentro de cada intervalo de tiempo hasta conseguir un valor de las tensiones que responda mejor que el del intervalo anterior, tal como se realiza en los métodos de integración paso a paso. Obsérve~e también que es relativamente sencillo modificar un pro-grama de cálculo con elementos finitos: basta añadir un lazo de instantes de tiempo y una subrutina que calcule las deformaciones viscoplásticas. 

Métodos numéricos

Mathematical Institute Eprints Archive

ENOC-2005, Eindhoven, Netherlands, 7-12 August 2005
UNSTABLE PERIODIC ORBITS OF
PERTURBED LORENZ EQUATIONS
Irene M. Moroz
Mathematical Institute
24-29 St Giles’, Oxford OX1 3LB
UK
moroz@maths.ox.ac.uk
Abstract
The extended Malkus-Robbins dynamo [Moroz,
2003] reduces to the Lorenz equations when one of the
key parameters, β, vanishes. In a recent study [Moroz,
2004] investigated what happened to the lowest order
unstable periodic orbits of the Lorenz limit as β was
increased to the end of the chaotic regime, using the
classic Lorenz parameter values of r = 28, σ = 10
and b = 8/3. In this paper we return to the parameter
choices of [Moroz, 2003], reporting on two of the cases
discussed therein.
Key words
Dynamos, unstable periodic orbits, Lorenz equations.
1 Introduction
Self-exciting dynamos can be nonlinear electro-
mechanical devices or naturally occurring magnetohy-
drodynamic (MHD) uid systems, such as the Earth’s
geodynamo, operating within the liquid metallic outer
core. These dynamos convert mechanical energy into
magnetic energy through the action of motional induc-
tion. Since the governing equations are fully nonlin-
ear partial differential equations in four independent
variables, valuable insights into on-going physical pro-
cesses can be achieved by investigating much lower di-
mensional models, which contain the key ingredients
for producing dynamo action, yet render themselves ac-
cessible for systematic study.
One such model is the extended Malkus-Robbins
dynamo (EMR) [Moroz, 2003], a four-mode system
of nonlinear coupled ordinary differential equations,
which reduce to the Malkus-Robbins dynamo [Rob-
bins, 1977] when the parameter β, measuring the in-
verse moment of inertia of the armature of the dynamo
motor, vanishes. The EMR dynamo is of interest since
it becomes the Lorenz system in this limit.
There has been considerable interest in recent years in
characterising chaotic attractors of nonlinear dynami-
cal systems, such as the Lorenz equations, in terms of
the spectrum of unstable periodic orbits (upos), deter-
mined from time series. While a chaotic attractor pos-
sesses an innite number of upos, it is believed that
many properties may be determined from those orbits
of lowest period.
In a recent study [Moroz, 2004] investigated what
happened to the lowest order upos of the Lorenz limit
of the EMR dynamo as β was increased, using the clas-
sic Lorenz parameter choices of r = 28, σ = 10 and
b = 8/3, for two cases considered in [Moroz, 2003]. In
one case it was possible to trace the evolution of some
of the leading upos through to the onset of stable peri-
odic oscillations. In the other case, the effects of non-
zero β rapidly destroyed the simple characeter of the
Poincar·e return map.
In this paper, we extend the study of [Moroz, 2004]
and present the results of investigations of the EMR
dynamo in which the β = 0 limit does not yield chaotic
behaviour, by choosing different values for r, σ and b.
2 The Extended Malkus-Robbins Dynamo
The Malkus-Robbins dynamo comprises a conducting
disk, driven to rotate with dimensionless angular speed
W (t) by a constant torque, with a magnetic eld per-
pendicular to the disk inducing a radial e.m.f. A dimen-
sionless current Z(t) ows in the disk, is removed by
a ring of brushes and is fed to a coil and external load.
The coil is aligned in such a way that its dimension-
less current Y (t) produces a magnetic eld to reinforce
the original eld. If U(t) is the angular velocity of the
motor, then the EMR dynamo takes the form [Moroz,
2003; Moroz, 2004],
Y ′ = σ(Z − Y )− βˆU,
Z ′ = Y [(1 + βˆΛσ )−W ]− Z,
W ′ = ZY − ν(W + Wˆe),
U ′ = Y − ΛU,
 (1)
where
W = (1 + βˆΛσ )−W,
We = R/ν − (1 + βˆΛσ ),
}
(2)
and all parameters appearing in (1) are real and posi-
tive. R measures the applied couple; ν is the re-scaled
mechanical friction; σ depends upon the inductances
and resistances of the coil and brushes; Λ represents
the mechanical friction of the motor and βˆ measures
the inverse moment of inertia of the armature of the
motor.
2.1 The Lorenz Limit
When βˆ = 0, the evolution of U decouples from the
rst three equations of (1), which reduce to the Lorenz
equations under the identication
(Y, Z,W +We, U) 7−→ (x, y, z, 0),
(σ,R/ν, ν, βˆ) 7−→ (σ, r, b, 0).
}
(3)
Our choice for the Poincar·e section is motivated by
[Koga, 1986], who, for the classic parameter choice of
r = 28, b = 8/3 and σ = 10 dened his Poincar·e
section as
S = [(x, y) : zˆ = 0, ˙ˆz > 0, x > 0], (4)
where zˆ = z − (r − 1). He was able to construct two
one-dimensional return maps
Xn+1 = f1(Xn), Yn+1 = f2(Yn), (5)
whose xed points corresponded to the upos of the se-
ries Ln1 for n = 1, 2, .... Lnm refers to n rotations
around the left focus and m around the right focus,
where the foci are the nontrivial equilibrium solutions
of the Lorenz equations. Upos with Lnm for m > 1,
are obtained by iterating the return maps m times and
nding their xed points.
In anticipation of what is to come, we remind the
reader that the subcritical Hopf bifurcation for the non-
trivial equilibria in the Lorenz system occurs when
rH =
σ(σ + b+ 3)
σ − 1− b (6)
so that σ > 1 + b must be satised for its existence.
2.2 The Poincar·e section for the EMR dynamo
By analogy with the Lorenz equations, we have cho-
sen the Poincar·e section for the EMR dynamo to be:
Sˆ = [(Y, Z, U) :W = 0, W˙ > 0, Y > 0]. (7)
The nontrivial equilibrium states are then We =
(Ye, Ze, 0, Ue), where
Ye = ±
√
νWˆe
(1+ βˆΛσ )
,
Ze = Ye(1 +
βˆ
Λσ ),
Ue = Ye/Λ,
 (8)
while the trivial equilibrium state is
(Y0, Z0,W0, U0) = (0, 0,−Wˆe, 0). (9)
We remark that a double-zero Takens-Bogdanov bifur-
cation occurs when
(Rd, βˆd) = (
ν(Λ + σ)
σ(1− Λ) ,
Λ2(1 + σ)
1− Λ ), (10)
provided Λ < 1, while the subcritical Hopf bifurcation
condition in the Lorenz limit translates to σ > 1 + ν.
2.3 Parameter Regimes
As discussed in detail in [Moroz, 2003], there are four
different parameter regimes of interest: Λ ≷ 1 and
σ ≷ 1 + ν [Moroz, 2003]. The conditions on Λ corre-
spond to the existence or otherwise of a codimension-
two double-zero bifurcation, associated with each of
the equilibrium solutions (see [Moroz, 2003] for de-
tails). The conditions on σ relate to the termination
or otherwise of the limit cycle, which emerges from
the multiple bifurcation point, on the subcritical Hopf
bifurcation point, associated with the Lorenz limit of
βˆ = 0. [Moroz, 2004] chose to work with the classic
Lorenz parameter values, which correspond to σ > 1.
In this contribution, we consider the alternative choice
of σ < 1, with Λ ≷ 1. This necessitates a different
selection for r, σ and b, and the ones we chose were
those in [Moroz, 2003] of ν = b = 1, σ = 1.9 and
either Λ = 0.5 or Λ = 1.2. These are cases (iii) and
(iv) of [Moroz, 2003].
2.4 Λ < 1 and 0 < σ < 1 + ν
When Λ = 0.5, ν = 1 and σ = 1.9, a linear stability
analysis shows the curves of Figure 3 of [Moroz, 2003].
The are curves of steady state solutions and curves of
Hopf bifurcation, emanating from the single double-
zero bifurcation point. Because of the violation of the
Hopf bifurcation criterion in the βˆ = 0 limit, the curve
of subcritical limit cycles, asymptotes for the βˆ = 0
axis, but does not touch it.
Numerical integrations reveal a small region of pa-
rameter space to be occupied by chaotic solutions, with
most of the rest dominated by simple periodic be-
haviour [Moroz, 2003].
2.5 Λ > 1 and 0 < σ < 1 + ν
When Λ = 1.2, ν = 1 and σ = 1.9, neither the
double-zero nor the subcritical Hopf bifurcation ex-
ist. A linear stability analysis yields one curve of
steady state bifurcations and a curve of (subcritical)
Hopf bifurcations (associated with the nontrivial equi-
libria) which no longer intersect (see Figure 5 of [Mo-
roz, 2003]).
Numerical integrations there show two regions of
chaotic behaviour, separated by a relatively small re-
gion of periodic behaviour, which increases asR and βˆ
increase. In this and in the other case considered here,
chaos is found near the left hand Hopf stability bound-
ary, marking the transition to steady state solutions.
3 Numerical integrations
In the numerical investigations, we integrated the
EMR dynamo equations for 60, 000 seconds using time
steps of either 0.01, 0.005 or 0.001 seconds, discarding
the rst 100 seconds as representing transients. The se-
lection of a close return on the Poincar·e sectionW = 0,
which corresponds to takingW = (1+ βˆΛσ ) in equation
(2), was accomplished using the criterion
||Wi −Wj || <  (11)
where Wm = (Ym, Zm,Wm, Um) and m = i, j de-
note the ith and jth intersections. We chose  = 0.05
to dene the ball of close returns.
We were able to guarantee that the trajectories land
precisely upon the Poincar·e section, regardless of the
size of the integration step, by adopting the approach
described in [Henon, 1982]. This is achieved by re-
writing the EMR equations withW as the independent
variable instead of t, and integrating the resulting equa-
tions for 1 step fromWa toW = 0. Here ta is the time
just beforeW changes sign.
Histograms of periods of the upos were constructed
based upon close returns to within  of the xed points
of the mth iterates of the return map for m = 1, ..., 10.
In our numerical integrations we used the
same initial conditions of (Y, Z,W , U) =
(−2.159025,−6.5509769,−8.7136071,−4.81492).
3.1 Λ < 1 and 0 < σ < 1 + ν
Figure 3 of [Moroz, 2003] shows the linear and non-
linear stability regions for Λ = 0.5 and σ = 1.9 be-
tween the ranges of 0 ≤ R ≤ 20 and 0 ≤ β ≤ 10. For
our current investigations, we chose R = 20, which,
for β increasing, yields the onset of chaotic behaviour
at β ≈ 0.39 as a sudden jump up from the zero equi-
librium state of no dynamo action. For β decreasing,
the loss occurs at β ≈ 0.31, again as a sudden jump,
demonstrating the presence of hysteresis effects. Such
behaviour is consistent with the subcritical nature of
the Hopf bifurcation curve which no longer terminates
on β = 0. The loss of chaotic behaviour is observed at
0 50 100 150 200 250 300
−8
−6
−4
−2
0
2
4
6
8
10
t
Y
Figure 1. Time series of Y for β = 0.4
4.5 5 5.5 6 6.5 7
4.5
5
5.5
6
6.5
7
Y
n
Y n
+
1
Figure 2. The first return map in (Yn, Yn+1)-space β = 0.4
β ≈ 2.34, via an intermittent destabilisation of a peri-
odic orbit.
[Moroz, 2004] found that for very small perturbations
about the Lorenz limit of β = 0 and for the Lorenz
parameter choices for σ and ν, the Λ < 1 case ap-
pears to be far more sensitive to the number of close
returns than does the case of Λ > 1. We have also
found the same to hold for the current study. When
β = 0.4, we found only 65 close returns with 60 mil-
lion time steps, the lowest order upo having a period
of 3.96sec. Figure 1 shows the time series of Y (t) for
β = 0.4, while Figure 2 shows the rst return map in
the (Yn, Yn+1)-plane. As β increases into the chaotic
regime, the number of close returns increases, while
the period of the dominant upo decreases. For exam-
ple when β = 2.0, we counted 2800 close returns with
the only observed lowest order upo having a period of
3.33sec. Figure 1 shows a section of the time series for
Y and Figure 2 the lowest order upo in the (Y, Z)-plane
for β = 2. The time series shows evidence of an inter-
mittent destabilisation of an underlying periodic orbit,
corroborated in the corresponding phase portrait (not
shown here). As the chaotic regime gives way to stable
0 50 100 150 200 250 300
−8
−6
−4
−2
0
2
4
6
8
10
t
Y
Figure 3. Time series of Y for β = 2.0
−8 −6 −4 −2 0 2 4 6 8
−15
−10
−5
0
5
10
15
Y
Z
Figure 4. Lowest order upo in (Y, Z)-space for β = 2.0
periodic limit cycle behaviour, we found, at β = 2.35,
the limit cycle’s period to be 3.259sec.
3.2 Λ > 1 and 0 < σ < 1 + ν
Figure 4 of [Moroz, 2003] shows the linear and non-
linear stability curves for σ = 1.9 with Λ = 1.2
between the ranges 0 ≤ β ≤ 100 and 0 ≤ R ≤
50. Numerical integrations showed the existence of
two regions of chaotic behaviour, contained within the
domain bounded by the subcritical Hopf bifurcation
curve. The larger region occurred for small values of
β, loss of stability being a discontinuous jump down to
steady (non-trivial) equilibrium states. The smaller re-
gion occurred in an expanding band, for larger values of
β, near the righthand boundary of the subcritical Hopf
curve. Figure 13 of [Moroz, 2003] gave a typical ex-
ample of a bifurcation transition curve for R = 50. For
the current study, we chose R = 20, as in the previous
subsection. In this case, the initial onset of the larger
chaotic band occurs at β ≈ 0.93, terminating in a sta-
ble periodic limit cycle when β = 4.6. This periodic
limit cycle maintains its stability, growing in amplitude
as β increases, until the onset of the second band of
−8 −6 −4 −2 0 2 4 6
−15
−10
−5
0
5
10
15
Y
Z
Figure 5. An L21 upo for β = 2.0
−8 −6 −4 −2 0 2 4 6 8
−15
−10
−5
0
5
10
15
Y
Z
Figure 6. An L22 upo for β = 2.0
chaotic solutions at β ≈ 23.8, with its subsequent loss
at β ≈ 24.1.
Following [Koga, 1986; Moroz, 2004], we adopt the
terminology Lnm, which refers to the trajectory rotat-
ing n times around the left focus and m times the right
focus, where the foci are the nontrivial equilibrium so-
lutions of equation (8). The basic upo of the previous
subsection would therefore be labelled L11.
When β = 2 we observed 857 close returns. The ba-
sic upo with period ≈ 3.79sec resembles that shown in
Figure 2 above. Figures 3 and 4 show examples of L21
and L22 upos with periods 5.8sec and 7.58sec respec-
tively. When β = 4 there were 6031 close returns with
the L11 upo of period 3.454sec accounting for nearly
all of the Lm1 series. Examples of L22 and L33 are
shown in Figures 5 and 6 respectively. The examples
shown had periods of 8.098sec and 11.74sec. When
β = 24, we observed 7520 close returns with ≈ 1800
of those centred about a basic upo of period 9.1sec.
Figure 6 shows a typical L11 upo. As in the Λ > 1
case in [Moroz, 2004], the Poincar·e return maps show
a cleaner banded structure than in the Λ < 1 case (see
for example Figure 7 for β = 1.0). As β increases,
−8 −6 −4 −2 0 2 4 6 8
−15
−10
−5
0
5
10
15
Y
Z
Figure 7. An L22 upo for β = 4.0
−8 −6 −4 −2 0 2 4 6 8
−15
−10
−5
0
5
10
15
Y
Z
Figure 8. An L33 upo for β = 4.0
more of the bands disappear until very few remain (see
Figure 8) prior to the loss of chaotic dynamo action to
a steady state.
4 Discussion
As noted in [Moroz, 2004], the purpose of study-
ing such low-dimensional dynamos lies in their incor-
poration of some of the key features of much larger-
scale naturally-occurring MHD dynamos. They exhibit
steady dynamo action, irregular reversals as well as
sustained periodic/multiply periodic oscillations. How-
ever with so many different low order dynamos of this
family in existence [Moroz, 2001], the issue naturally
arises as to how to distinguish between them. The iden-
tication of the spectrum of unstable periodic orbits, as
well as their behaviour when key parameters vary pro-
vides one means of comparison, especially if similar
approaches are applied to identify upos of the large-
scale dynamo models. The present study is one more
step in that direction.
−3 −2 −1 0 1 2 3
−20
−15
−10
−5
0
5
10
15
20
Y
Z
Figure 9. An L11 upo for β = 24.0
5 5.5 6 6.5 7
5
5.5
6
6.5
7
Y
n
Y n
+
1
Figure 10. First return map for β = 1.0 with Λ = 1.2
References
Henon, M. (1982) On the numerical computation of
Poincar·e maps. Physica D 5 pp. 412-414.
Koga, S. (1986) Phase description method to time aver-
ages in the Lorenz system. Prog. Theor. Phys. 76 pp.
335-355.
Moroz, I. M. (2001). Self-exciting Faraday disk ho-
mopolar dynamos. Int. J. Bifur. Chaos. 11 pp. 2961-
75.
Moroz, I. M. (2003). The Malkus-Robbins dynamo
with a linear series motor. Int. J. Bifur. Chaos. 13 pp.
147-161.
Moroz, I. M. (2004). The extended Malkus-Robbins
dynamo as a perturbed Lorenz system. Nonlinear
Processes. In the press.
Robbins, K. A. (1977) A new approach to subcritical
instability and turbulent transitions in a simple dy-
namo. Math. Proc. Camb. Phil. Soc. 82 pp. 309-325.
2.35 2.4 2.45 2.5 2.55 2.6
2.35
2.4
2.45
2.5
2.55
2.6
Figure 11. First return map for β = 24.0 with Λ = 1.2


A new approach to subcritical instability and turbulent transitions in a simple dynamo.

On the numerical computation of Poincar·e maps.

Phase description method to time averages in the Lorenz system.

Self-exciting Faraday disk homopolar dynamos.

The extended Malkus-Robbins dynamo as a perturbed Lorenz system. Nonlinear Processes. In

The Malkus-Robbins dynamo with a linear series motor.

Unstable periodic orbits of perturbed Lorenz equations

En primer lugar se hace referencia al Método de los Elementos Finitos (MEF), a partir del cual se plantea el  esquema general de un programa en elasticidad en el que, para empezar, se procede a la introducción de datos; despues, para cada elemento se realiza el cálculo de k elevado a e, el cálculo de f elevado a e, la rotación a coordenadas globales y la colocación en matrices globales k y f; los siguientes pasos son: aplicación de las condiciones de contorno, resolución del sistema de ecuaciones, cálculo de tensiones y esfuerzos y finalmente, presentación de resultados. En segundo lugar, en el apartado dedicado a la viscoelasticidad se recuerdan dos alternativas que son de interés en el tratamiento numérico de problemas viscoelásticos. Por un lado se resume la formulación de la ley de comportamiento mediante variables de estado y a continuación se presenta una formulación de tipo potencial. Finalmente se muestra un ejemplo sencillo

Alarcón Álvarez, Enrique

Archivo Digital UPM

http://oa.upm.es/25887/1/Metodos_Numericos.pdf

Unstable periodic orbits of perturbed Lorenz equations

Abstract

Similar works

Full text

Available Versions

Repositorio da Universidade da Coruña

Mathematical Institute Eprints Archive

Archivo Digital UPM