Category : C Source Code
Archive   : DYNAM_S1.ZIP
Filename : YDES.C

/******************* (C) 1986,7,8 by JAMES A. YORKE **************************/
/******************************* YDEs.C ***************************************
This file is where the differential equations are defined. They are
called (generally) by DifEqu() which is called by rungekutta(),
and rungekutta() is in YTOOLS.C. Notice that Y[] is a temporary
position vector which is set by rungekutta to be either y[] or temp[].
DED(k,Y) DUFFING: x'' + C1*x' -C2*x +C3*x*x*x = rho*sin(C4*t)
initDED()
DEPD(k,Y) PARAMETRICALLY DRIVEN DUFFING: x'' + C1*x'
-x + C2*(1+C5*sin(C4*t))*x +C3*(1+C6*sin(C4*t))*x*x*x = 0
initDED()
DEP(k,Y) Forced Damped Pendulum x''+C1*x'+C2*sinx =rho*(C3+cos(C4*t))
initDEP()
DEV(k,Y) V: a VanDerPol differential equation due to UEDA
initDEV() for V: a VanDerPol differential equation due to UEDA

********************** (C) 1986,7,8 by JAMES A. YORKE ************************/
#include "yinclud.h"





/* Y varies in different stages of the rungekutta
and is usually not equal to y */

DED(k,Y) /* DUFFING: x'' + C1*x' -C2*x +C3*x*x*x = rho*sin(C4*t) */
double k[],Y[];
{
int base,i;
double t,x,y,xx,xxx;
t = Y[0];
x = Y[1];
y = Y[2];
xx = x*x;
xxx = x*xx;

k[0] = 1;
k[1] = y ;
k[2] = -C1*y + C2*x -C3*xxx + rho*sin(C4*t) ;
k[3] = 1; /* coordinate 0 is set to 0 after each cycle */
if (num_lyap > 0)
{
for (i = 0; i < num_lyap; i++)
{
base = lyapzero + vec_dim*i;
k[base] = Y[ base + 1];
k[base +1 ] = (C2 - 3*C3*xx) *Y[ base ] -C1*Y[ base + 1];
}
}
}
initDED()
{
int d_mod();
zeroth = 1;
vec_dim = 2; /*the dimension of the
Lyapunov vectors = phase space dim*/
num_lyap = 0; /*the number of Lyapunov numbers
to be computed <= vec_dim */
lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
/* variables: time mod twopi/C4 ,x,y,t */

X_upper = 2.2 ; /* x scale */
X_lower = -2.2;
Y_upper = 2.; /* y scale */
Y_lower = -2.;
X_coord = 1;
Y_coord = 2;
C1 = 1.;
C2 = 10;
C3 = 100;
rho = .85;
C4 = 3.5;
modPointer = d_mod;
DE = DED; /* DE is a pointer to a function */
steps_per_cycle = 80;
/* now PickMap in YPickMap.C will automatically set:
step = (twopi/C4)/steps_per_cycle; */
}
d_mod()/* this is invoked by Duffing and by Parametrically Driven Duffing */
{
double moduloAB();
modFlag = NO;
y[0] = moduloAB(y[0],0.,twopi/C4);
}

/* DEPD PARAMETRICALLY DRIVEN DUFFING:
x'' +C1*x' -x +C2*(1+rho*sin(C4*t))*x^2 +C3*(1+rho*sin(C4*t))*x^3 = 0 */
DEPD(k,Y)
double k[],Y[];
{
int base,i;
double t,x,y,xx,xxx,s,xxCoef,xxxCoef;
t = Y[0];
x = Y[1];
y = Y[2];
xx = x*x;
xxx = x*xx;
s = 1 + rho*sin(C4*t);
xxCoef = C2*s;
xxxCoef = C3*s;
k[0] = 1;
k[1] = y ;
k[2] = -C1*y +x -xxCoef*xx -xxxCoef*xxx;
k[3] = 1; /* coordinate 0 is set to 0 after each cycle */
if (num_lyap > 0)
{
for (i = 0; i < num_lyap; i++)
{
base = lyapzero + vec_dim*i;
k[base] = Y[base + 1];
k[base +1 ] = (1. -2.*xxCoef*x - 3.*xxxCoef*xx) *Y[ base ]
-C1*Y[ base + 1];
}
}
}

initDEPD()/* Parametrically Driven Duffing */
{
int d_mod();
zeroth = 1;
vec_dim = 2; /*the dimension of the
Lyapunov vectors = phase space dim*/
num_lyap = 0; /*the number of Lyapunov numbers
to be computed <= vec_dim */
lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
/* variables: time mod twopi/C4 ,x,y,t */

X_lower = -5.;
X_upper = 3. ; /* x scale */
Y_lower = -3.;
Y_upper = 6.; /* y scale */
X_coord = 1;
Y_coord = 2;
C1 = 1.;
C2 = 1.;
C3 = 1.;
C4 = 1.;
rho = .5;
rho_step = .05;


modPointer = d_mod;
DE = DEPD; /* DE is a pointer to a function */
steps_per_cycle = 50;
/* now PickMap in YPickMap.C will automatically set:
step = (twopi/C4)/steps_per_cycle; */
}

/* y[0] = moduloAB(y[0],0.,1.);
y[1] = moduloAB(y[1],0.,1.);
*/




DEP(k,Y)/* Forced Damped Pendulum x''+C1*x'+C2*sinx =rho*(C3+cos(C4*t))*/
double k[],Y[];
{
int i;
double t,x,y,sinx,C2cosx;
int base;
t = Y[0];
x = Y[1];
y = Y[2];
sinx=sin(x);

k[0] = 1;
k[1] = y;
k[2] = -C1*y -C2*sinx + rho*(C3+cos(C4*t));
k[3] = 1; /* coordinate 0 is set to 0 after each cycle */

if (num_lyap > 0)
{
C2cosx = C2*cos(x);

for (i = 0; i < num_lyap; i++)
{
base = lyapzero + vec_dim*i;
k[base]
= Y[ base + 1];
k[base +1 ]
= -C2cosx *Y[ base ]
-C1*Y[ base + 1];
}
}
}


/* y[0] = moduloAB(y[0],0.,1.);
y[1] = moduloAB(y[1],0.,1.);
*/

v_mod()
{
double moduloAB();
y[0] = moduloAB(y[0],0.,twopi/C4);
modFlag = NO;
}
p_mod()
{
double moduloAB();
y[0] = moduloAB(y[0],0.,twopi/C4);/* time */
modFlag = NO;
if (X_coord != 1)
y[1] = moduloAB(y[1],-pi,pi);
else
y[1] = moduloAB(y[1],X_lower,X_lower + twopi);
}


y_mod_for_des()
{
return;
}


initDEP()
{
int mod_p();
zeroth = 1;
vec_dim = 2; /*the dimension of the
Lyapunov vectors = phase space dim*/
num_lyap = 0; /*the number of Lyapunov numbers
to be computed <= vec_dim */
lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
/* variables: time mod twopi/C4 ,x,y,t */

X_upper = pi ; /* x scale */
X_lower = -pi;
Y_upper = 4.; /* y scale */
Y_lower = -2.;
X_coord = 1;
Y_coord = 2;
C1 = .2;
C2 = 1.0;
C3 = 0.;
rho = 2.5;
C4 = 1.;
steps_per_cycle = 30;
period = 1;/* for Newton method */
its_per_plot = steps_per_cycle;
/* now PickMap in YPickMap.C will automatically set:
step = (twopi/C4)/steps_per_cycle; */
DE = DEP; /* DE is a pointer to a function */
modPointer = p_mod;/* takes angle mod two pi, more or less */
}


DEV(k,Y) /* V: a VanDerPol differential equation due to UEDA */
double k[],Y[];
{
int i,base;
double t,xx,x,y;
t = Y[0];
x = Y[1];
y = Y[2];
xx = x*x;
k[0] = 1;
k[1] = y;
k[2] = C1*y*(1-xx) -C2*x -C3*x*xx + rho*sin(C4*t);
k[3] = 1; /* coordinate 0 is set to 0 after each cycle */
if (num_lyap > 0)
{
for (i = 0; i < num_lyap; i++)
{
base = lyapzero + vec_dim*i;
k[base] = Y[ base + 1];
k[base +1 ] = (-2*x*C1*y -C2 - 3*C3*xx) *Y[ base ]
+C1*(1-xx)*Y[ base + 1];
}
}
}

initDEV()/* for V: a VanDerPol differential equation due to UEDA */
{
int v_mod();

zeroth = 1;
vec_dim = 2; /*the dimension of the
Lyapunov vectors = phase space dim*/
num_lyap = 0; /*the number of Lyapunov numbers
to be computed <= vec_dim */
lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
/* variables: time mod twopi/C4 ,x,y,t */
dim = lyapzero + num_lyap*vec_dim;
/* needed for rungekutta */

X_upper = 2. ; /* x scale */
X_lower = -2.;
Y_upper = 2.; /* y scale */
Y_lower = -2.;
X_coord = 1;
Y_coord = 2;
C1 = .2;
C2 = 0.;
C3 = 1.;
rho = 2.0;
C4 = 1.;
rho = 2.0;
steps_per_cycle = 40;
/* now PickMap in YPickMap.C will automatically set:
step = (twopi/C4)/steps_per_cycle; */
DE = DEV; /* DE is a pointer to a function */
modPointer = v_mod;/* specifies a modulo calculation */
}