The following program was used to explore the highly nonlinear and potentially chaotic
models that are analyzed in chapter 6 of Serpents in the Sand (University of Michigan
Press 1995). The program is written in Think Pascal, and it was run on a Macintosh
Quadra 900. It is capable of producing postscript output that can be sent directly to
a postscript printer, or input into a word processing program capable of accepting
postscript files (such as Word). The program is copyrighted 1995 by Courtney Brown.
However, researchers may copy, paste, and share this program freely for their own
use as long as the program itself is never sold. Simply copy and paste it into your
own program editor. To understand a bit more about how the program works, I
recommend my book, Chaos and Catastrophe Theories (Sage 1995).
program chaos (input, output, Image1, Image2, Image3);
const
minwidth = 15;
maxsize = 1000;
f1name = 'ImageLyap';
{f1name = 'ImageSdev';}
f2name = 'ImageSum';
f3name = 'ImagePower';
{nn=# of equations, nnn=total # of equations;}
matnn = 4;
matnnn = 20;
inc = 0.001;
vv2 = 2048;
nnpp = 9000;
gleft = 680;
gtop = 460;
gright = 1080;
gbottom = 860;
type
HugePtr = ^BigArray;
BigArray = array[1..200, 1..200] of real;
HugePtr2 = ^BigArray2;
BigArray2 = array[1..20000, 1..matnn] of real;
HugePtr3 = ^BigArray3;
BigArray3 = array[1..10000] of real;
GraphPtr = ^GraphBig;
GraphBig = array[1..20000, 1..2] of real;
RealArrayVV2 = array[1..vv2] of real;
RealArrayVV3 = array[1..1025, 1..2] of real;
RealArrayNP = array[1..nnpp] of real;
var
dim, c1, c2, d1, d2, iter, iternum, toss, lyap, kitu, limit, multiplier, alert, partychoice, choice, shocks: integer;
funcdem, funcrep, N, x1, x2, z1, z2, h, thing, min, max, weigh, totvote, time, perspect, CINI: real;
t, cumkt, upper, lower, decay, sense, limitpol, halflife, shift, policy, shock, varsum, varold, change: real;
lx1, ux1, lx2, ux2, xprim, lowp1, lowp2, highp1, highp2, chaoslagg, volume, Et, mold, rDR, timeone: real;
totderiv1, totderiv2, min2, max2, min3, max3, totderiv3, totpower, avefourier, varfourier, sdevfour: real;
ff1, ff2, ff3: text;
gizmo, demdum, repdum, row, col, ticker, ticker2, ticker3, lagg, yearslag, graphtype, traject, sasa: integer;
politics, algorithm, nn, nnn, eqnum, laggedid, count, item, fourier, anza, graphonly, totsignal: integer;
signal1, suppress, p1lag: integer;
filein: text;
y: HugePtr2;
p: array[1..20] of real;
party: array[1..4] of integer;
image1: HugePtr;
image2: HugePtr;
image3: HugePtr;
newx: GraphPtr;
fourdata: RealArrayVV2;
spectral: RealArrayVV3;
xone: array[1..1000] of real;
xlag: HugePtr3;
mx: array[1..matnnn] of real;
x: array[1..matnnn] of real;
znorm: array[1..matnn] of real;
gsc: array[1..matnn] of real;
cum: array[1..matnn] of real;
xnew: array[1..matnnn] of real;
finame: string;
r: Rect;
r2: Rect;
procedure initialize;
begin
r.left := 0;
r.top := 50;
r.right := 640;
r.bottom := 460;
SetTextRect(r);
r2.left := gleft;
r2.top := gtop;
r2.right := gright;
r2.bottom := gbottom;
SetDrawingRect(r2);
new(y);
new(newx);
new(xlag);
alert := 0;
ticker2 := 0;
yearslag := 0;
laggedid := 0;
end;
procedure parms;
begin
p[4] := 0.08847;
p[2] := 0.23665;
p[11] := 0.27780;
p[1] := 0.11699;
p[5] := 1.54509;
p[3] := 0.54922;
p[12] := 0.29619;
p[8] := 0.87930;
p[10] := 0.16780;
p[9] := 0.09635;
p[6] := 0.77792;
p[7] := 0.50896;
p[13] := 0.32081;
end;
function lag (lagvar: real): real;
var
iii: integer;
begin
if (ticker3 = 1) then
begin
for iii := 1 to lagg do
begin
xlag^[iii] := lagvar;
end;
lag := lagvar;
end;
if (ticker3 > 1) then
begin
for iii := (lagg - 1) downto 1 do
begin
xlag^[(lagg - iii)] := xlag^[(lagg - iii + 1)];
end;
xlag^[(lagg)] := lagvar;
lag := xlag^[1];
end;
end;
function power (base, exponent: real): real;
begin
if (base < 0) then
writeln('The power exponential function cannot use a negative base argument!');
if (base > 0) then
power := exp(exponent * ln(base));
if (base = 0) then
power := 0.0;
end;
function eq1: real;
begin
if (kitu = 1) then
eq1 := -(mx[1] / p[2]) - (sin(mx[2])) + (p[1] * cos(mx[3]));
if (kitu = 2) then
eq1 := p[1] * (mx[2] - mx[1]);
if (kitu = 3) then
begin
limitpol := 1;
if (politics = 2) then
policy := p[2];
if ((politics = 1) and (eqnum = 3)) then
begin
if (choice = 1) then
policy := 0.0;
if (choice = 0) then
policy := p[9];
end;
if ((politics = 1) and (eqnum = 4)) then
policy := mx[4];
eq1 := ((p[1] * mx[1]) * (1 - mx[1])) - (policy * mx[2] * mx[1]) - (decay * mx[1]);
end;
if (kitu = 4) then
eq1 := (p[1] * mx[1] * (1 - mold)) - (p[5] * mx[1] * mx[2]) - (p[2] * mx[1] * sqr(sin(mx[3])));
if (kitu = 5) then
eq1 := (((p[1] * mx[1]) + (p[2] * mx[2])) * (1 + (p[3] * CINI) + (p[13] * mx[3]))) * mx[1];
end;
function eq2: real;
begin
if (kitu = 1) then
eq2 := mx[1];
if (kitu = 2) then
eq2 := (p[2] * mx[1]) - (mx[1] * mx[3]) - mx[2];
if (kitu = 3) then
eq2 := (p[3] * mold) * (1 - mx[2]) - (p[4] * mx[3]);
if (kitu = 4) then
eq2 := (p[4] * mx[1] * (1 - mx[2])) - (p[6] * mx[2] * (1 - mx[1]));
if (kitu = 5) then
eq2 := (((p[4] * mx[1]) - (p[5] * mx[1] * mx[3])) * (1 + (p[6] * mx[2]))) * mx[2];
end;
function eq3: real;
begin
if (kitu = 1) then
eq3 := p[3];
if (kitu = 2) then
eq3 := (mx[1] * mx[2]) - (p[3] * mx[3]);
if (kitu = 3) then
eq3 := ((p[5] * mx[2] * mx[1]) * (1 - mx[3])) - (p[6] * mx[3]);
if (kitu = 4) then
eq3 := (2 * pi);
if (kitu = 5) then
eq3 := (p[7] + (p[8] * ln(abs((p[9] + (p[10] * mx[1]) + (p[11] * mx[2])) / mx[3]))) - (p[12] * mx[3])) * mx[3];
end;
function eq4: real;
var
partydum1, partydum2: real;
begin
if (kitu = 3) then
begin
if (choice = 1) then
begin
partydum1 := 1.0;
partydum2 := 0.0;
end;
if (choice = 0) then
begin
partydum1 := 0.0;
partydum2 := 1.0;
end;
eq4 := p[2] * mx[4] * ((p[7] * partydum1) + (p[8] * partydum2) - mx[4]);
end;
end;
procedure special3;
var
iii, tround: integer;
begin
if (timeone > 4.0) then
begin
partychoice := partychoice + 1;
party[1] := 1;
party[2] := 1;
party[3] := 0;
party[4] := 0;
choice := party[partychoice];
if (partychoice = 4) then
partychoice := 0;
timeone := 0.0;
end;
mold := lag(x[1]);
laggedid := 1;
Et := 0.1 * (Random / 32768.0);
{if (Et > 0) then Et := 0;}
{writeln('Et=', Et);}
{writeln('x1=', x[1] : minwidth, ' m1old= ', m1old : minwidth, 'lag= ', lagg);}
tround := round(time);
if ((shocks = 1) and (tround = (toss div 2))) then
x[1] := x[1] + shock;
end;
procedure special4;
begin
mold := lag(x[1]);
laggedid := 1;
end;
procedure special5;
begin
if (ticker3 = 1) then
begin
x[1] := 0.0;
x[2] := 0.86;
x[3] := 0.02;
CINI := 0.23;
end;
if (ticker3 = 11) then
begin
x[1] := 0.0;
CINI := 0.28;
end;
if ((ticker3 > 20) and (x[1] > 0.25)) then
begin
x[1] := 0.0;
CINI := 0.25;
end;
end;
function par1x1: real;
var
it1, it2: real;
begin
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
it1 := eq1;
mx[1] := mx[1] + (2 * inc);
if (laggedid = 1) then
mold := mold + (2 * inc);
it2 := eq1;
par1x1 := (it2 - it1) / (2 * inc);
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
end;
function par1x2: real;
var
it1, it2: real;
begin
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
it1 := eq1;
mx[2] := mx[2] + (2 * inc);
if (laggedid = 2) then
mold := mold + (2 * inc);
it2 := eq1;
par1x2 := (it2 - it1) / (2 * inc);
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
end;
function par1x3: real;
var
it1, it2: real;
begin
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
it1 := eq1;
mx[3] := mx[3] + (2 * inc);
if (laggedid = 3) then
mold := mold + (2 * inc);
it2 := eq1;
par1x3 := (it2 - it1) / (2 * inc);
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
end;
function par1x4: real;
var
it1, it2: real;
begin
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
it1 := eq1;
mx[4] := mx[4] + (2 * inc);
if (laggedid = 4) then
mold := mold + (2 * inc);
it2 := eq1;
par1x4 := (it2 - it1) / (2 * inc);
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
end;
function par2x1: real;
var
it1, it2: real;
begin
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
it1 := eq2;
mx[1] := mx[1] + (2 * inc);
if (laggedid = 1) then
mold := mold + (2 * inc);
it2 := eq2;
par2x1 := (it2 - it1) / (2 * inc);
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
end;
function par2x2: real;
var
it1, it2: real;
begin
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
it1 := eq2;
mx[2] := mx[2] + (2 * inc);
if (laggedid = 2) then
mold := mold + (2 * inc);
it2 := eq2;
par2x2 := (it2 - it1) / (2 * inc);
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
end;
function par2x3: real;
var
it1, it2: real;
begin
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
it1 := eq2;
mx[3] := mx[3] + (2 * inc);
if (laggedid = 3) then
mold := mold + (2 * inc);
it2 := eq2;
par2x3 := (it2 - it1) / (2 * inc);
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
end;
function par2x4: real;
var
it1, it2: real;
begin
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
it1 := eq2;
mx[4] := mx[4] + (2 * inc);
if (laggedid = 4) then
mold := mold + (2 * inc);
it2 := eq2;
par2x4 := (it2 - it1) / (2 * inc);
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
end;
function par3x1: real;
var
it1, it2: real;
begin
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
it1 := eq3;
mx[1] := mx[1] + (2 * inc);
if (laggedid = 1) then
mold := mold + (2 * inc);
it2 := eq3;
par3x1 := (it2 - it1) / (2 * inc);
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
end;
function par3x2: real;
var
it1, it2: real;
begin
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
it1 := eq3;
mx[2] := mx[2] + (2 * inc);
if (laggedid = 2) then
mold := mold + (2 * inc);
it2 := eq3;
par3x2 := (it2 - it1) / (2 * inc);
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
end;
function par3x3: real;
var
it1, it2: real;
begin
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
it1 := eq3;
mx[3] := mx[3] + (2 * inc);
if (laggedid = 3) then
mold := mold + (2 * inc);
it2 := eq3;
par3x3 := (it2 - it1) / (2 * inc);
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
end;
function par3x4: real;
var
it1, it2: real;
begin
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
it1 := eq3;
mx[4] := mx[4] + (2 * inc);
if (laggedid = 4) then
mold := mold + (2 * inc);
it2 := eq3;
par3x4 := (it2 - it1) / (2 * inc);
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
end;
function par4x1: real;
var
it1, it2: real;
begin
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
it1 := eq4;
mx[1] := mx[1] + (2 * inc);
if (laggedid = 1) then
mold := mold + (2 * inc);
it2 := eq4;
par4x1 := (it2 - it1) / (2 * inc);
mx[1] := mx[1] - inc;
if (laggedid = 1) then
mold := mold - inc;
end;
function par4x2: real;
var
it1, it2: real;
begin
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
it1 := eq4;
mx[2] := mx[2] + (2 * inc);
if (laggedid = 2) then
mold := mold + (2 * inc);
it2 := eq4;
par4x2 := (it2 - it1) / (2 * inc);
mx[2] := mx[2] - inc;
if (laggedid = 2) then
mold := mold - inc;
end;
function par4x3: real;
var
it1, it2: real;
begin
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
it1 := eq4;
mx[3] := mx[3] + (2 * inc);
if (laggedid = 3) then
mold := mold + (2 * inc);
it2 := eq4;
par4x3 := (it2 - it1) / (2 * inc);
mx[3] := mx[3] - inc;
if (laggedid = 3) then
mold := mold - inc;
end;
function par4x4: real;
var
it1, it2: real;
begin
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
it1 := eq4;
mx[4] := mx[4] + (2 * inc);
if (laggedid = 4) then
mold := mold + (2 * inc);
it2 := eq4;
par4x4 := (it2 - it1) / (2 * inc);
mx[4] := mx[4] - inc;
if (laggedid = 4) then
mold := mold - inc;
end;
procedure equations (k: integer);
var
i: integer;
begin
if (eqnum = 3) then
begin
if (k = 1) then
xprim := eq1;
if (k = 2) then
xprim := eq2;
if (k = 3) then
xprim := eq3;
if (k > 3) then
begin
if (k < 7) then
begin
i := k - 4;
xprim := (mx[4 + i] * par1x1) + (mx[7 + i] * par1x2) + (mx[10 + i] * par1x3);
end;
end;
if (k > 6) then
begin
if (k < 10) then
begin
i := k - 7;
xprim := (mx[4 + i] * par2x1) + (mx[7 + i] * par2x2) + (mx[10 + i] * par2x3);
end;
end;
if (k > 9) then
begin
i := k - 10;
xprim := (mx[4 + i] * par3x1) + (mx[7 + i] * par3x2) + (mx[10 + i] * par3x3);
end;
end;
if (eqnum = 4) then
begin
if (k = 1) then
xprim := eq1;
if (k = 2) then
xprim := eq2;
if (k = 3) then
xprim := eq3;
if (k = 4) then
xprim := eq4;
if (k > 4) then
begin
if (k < 9) then
begin
i := k - 5;
xprim := (mx[5 + i] * par1x1) + (mx[9 + i] * par1x2) + (mx[13 + i] * par1x3) + (mx[17] * par1x4);
end;
end;
if (k > 8) then
begin
if (k < 13) then
begin
i := k - 9;
xprim := (mx[5 + i] * par2x1) + (mx[9 + i] * par2x2) + (mx[13 + i] * par2x3) + (mx[17] * par2x4);
end;
end;
if (k > 12) then
begin
if (k < 17) then
begin
i := k - 13;
xprim := (mx[5 + i] * par3x1) + (mx[9 + i] * par3x2) + (mx[13 + i] * par3x3) + (mx[17] * par3x4);
end;
end;
if (k > 16) then
begin
i := k - 17;
xprim := (mx[5 + i] * par4x1) + (mx[9 + i] * par4x2) + (mx[13 + i] * par4x3) + (mx[17] * par4x4);
end;
end;
end;
procedure model;
var
k: integer;
rk1, rk2, rk3, rk4: array[1..matnnn] of real;
xx1, xx2, xx3, xx4: array[1..matnnn] of real;
begin
{The fourth order Runge-Kutta;}
for k := 1 to nnn do
begin
mx[k] := x[k];
end;
for k := 1 to nnn do
begin
equations(k);
rk1[k] := h * xprim;
end;
for k := 1 to nnn do
begin
xx1[k] := x[k] + (rk1[k] * 0.5);
end;
for k := 1 to nnn do
begin
mx[k] := xx1[k];
end;
for k := 1 to nnn do
begin
equations(k);
rk2[k] := h * xprim;
end;
for k := 1 to nnn do
begin
xx2[k] := x[k] + (rk2[k] * 0.5);
end;
for k := 1 to nnn do
begin
mx[k] := xx2[k];
end;
for k := 1 to nnn do
begin
equations(k);
rk3[k] := h * xprim;
end;
for k := 1 to nnn do
begin
xx3[k] := x[k] + rk3[k];
end;
for k := 1 to nnn do
begin
mx[k] := xx3[k];
end;
for k := 1 to nnn do
begin
equations(k);
rk4[k] := h * xprim;
end;
for k := 1 to nnn do
begin
xnew[k] := x[k] + ((1 / 6) * (rk1[k] + (2 * rk2[k]) + (2 * rk3[k]) + rk4[k]));
end;
end;
procedure initcond;
var
ii: integer;
begin
{Initial conditions for the nonlinear system;}
if (kitu = 1) then
begin
x[1] := 0.5;
x[2] := 0.0;
x[3] := 0.0;
end;
if (kitu = 2) then
begin
x[1] := 5;
x[2] := 5;
x[3] := 5;
end;
if (kitu = 3) then
begin
x[1] := 0.1;
x[2] := 0.1;
x[3] := 0.1;
if (eqnum = 4) then
x[4] := 0.1;
end;
if (kitu = 4) then
begin
x[1] := 0.1;
x[2] := 0.1;
x[3] := 0.1;
end;
if (kitu = 5) then
begin
x[1] := 0.0;
x[2] := 0.86;
x[3] := 0.02;
end;
{Initial conditions fo rthe orthonormal frame;}
for ii := (nn + 1) to nnn do
begin
x[ii] := 0.0;
end;
for ii := 1 to nn do
begin
x[(nn + 1) * ii] := 1.0;
cum[ii] := 0.0;
end;
end;
procedure initnums;
begin
mold := 0.0;
partychoice := 0;
ticker3 := 0;
time := 0.0;
timeone := 0.0;
count := 0;
volume := 0.0;
Getdatetime(randSeed);
varsum := 0.0;
fourier := 1;
totpower := 0.0;
totsignal := 0;
end;
procedure avevar (data2: RealArrayVV3; uu: integer; var ave, svar: real);
var
L: integer;
s: real;
weightsum: real;
begin
ave := 0.0;
svar := 0.0;
weightsum := 0.0;
for L := 1 to uu do
begin
weightsum := weightsum + data2[L, 2];
ave := ave + (data2[L, 1] * data2[L, 2]);
{writeln('ave=', ave, ' weightsum=', weightsum);}
end;
ave := ave / weightsum;
for L := 1 to uu do
begin
s := data2[L, 1] - ave;
svar := svar + (s * s * data2[L, 2] / weightsum);
end;
end;
procedure four1 (var data: RealArrayVV2; vv, isign: integer);
var
ww, qq, v, mmax, m, q, istep, w: integer;
wtemp, wr, wpr, wpi, wi, theta: double;
tempr, tempi, wrs, wis: real;
begin
v := 2 * vv;
q := 1;
for ww := 1 to vv do
begin
w := 2 * ww - 1;
if (q > w) then
begin
tempr := data[q];
tempi := data[q + 1];
data[q] := data[w];
data[q + 1] := data[w + 1];
data[w] := tempr;
data[w + 1] := tempi;
end;
m := v div 2;
while (m >= 2) and (q > m) do
begin
q := q - m;
m := m div 2;
end;
q := q + m;
end;
mmax := 2;
while (v > mmax) do
begin
istep := 2 * mmax;
theta := 6.28318530717959 / (isign * mmax);
wpr := -2.0 * sqr(sin(0.5 * theta));
wpi := sin(theta);
wr := 1.0;
wi := 0.0;
for ww := 1 to mmax div 2 do
begin
m := 2 * ww - 1;
wrs := wr;
wis := wi;
for qq := 0 to (v - m) div istep do
begin
w := m + qq * istep;
q := w + mmax;
tempr := wrs * data[q] - wis * data[q + 1];
tempi := wrs * data[q + 1] + wis * data[q];
data[q] := data[w] - tempr;
data[q + 1] := data[w + 1] - tempi;
data[w] := data[w] + tempr;
data[w + 1] := data[w + 1] + tempi;
end;
wtemp := wr;
wr := wr * wpr - wi * wpi + wr;
wi := wi * wpr + wtemp * wpi + wi;
end;
mmax := istep;
end;
for ww := 1 to (vv div 2) do
begin
spectral[ww, 2] := sqr(data[((ww * 2) - 1)]) + sqr(data[(ww * 2)]);
spectral[ww, 1] := (ww / (vv * h));
if (suppress = 1) then
spectral[1, 2] := 0.0;
if (ww < 20) then
begin
writeln('freq=', spectral[ww, 1] : 12 : 5, ' periodogram= ', spectral[ww, 2] : 15 : 3);
{writeln('fourier numbers 1=', data[((ww * 2) - 1)] : 8 : 6, ' fourier numbers 2=', data[(ww * 2)] : 8 : 6);}
end;
totpower := totpower + spectral[ww, 2];
end;
totpower := totpower / vv;
writeln('Total spectral power=', totpower : 20 : 6);
avevar(spectral, (vv div 2), avefourier, varfourier);
sdevfour := sqrt(varfourier);
writeln('avefourier=', avefourier : 10 : 3, ' sdevfour=', sdevfour : 10 : 3);
end;
procedure graph;
var
k, j, width, height, transx1, transx2, switch, kwisha, fourmax, transform: integer;
minx1, minx2, maxx1, maxx2, miny1, miny2, maxy1, maxy2, miny3, maxy3, axedraw1, axedraw2: real;
axeminx, axemaxx, axeminy, axemaxy, scale: real;
axedraw: array[1..3, 1..2] of real;
begin
width := gright - gleft;
height := gbottom - gtop;
if (algorithm = 1) then
kwisha := multiplier - 1;
if (algorithm = 2) then
kwisha := multiplier + 5;
EraseRect(0, 0, 400, 400);
if (fourier = 1) then
begin
transform := 1;
if ((multiplier - anza) >= 1024) then
fourmax := 1024;
if (((multiplier - anza) >= 512) and ((multiplier - anza) < 1024)) then
fourmax := 512;
if (((multiplier - anza) >= 256) and ((multiplier - anza) < 512)) then
fourmax := 256;
if (((multiplier - anza) >= 128) and ((multiplier - anza) < 256)) then
fourmax := 128;
for k := 1 to fourmax do
begin
fourdata[((k * 2) - 1)] := y^[((anza - 1) + k), 1];
fourdata[(k * 2)] := 0.0;
end;
four1(fourdata, fourmax, transform);
end;
for k := anza to (multiplier - 1) do
begin
if (k = anza) then
begin
miny1 := y^[anza, 1];
maxy1 := y^[anza, 1];
miny2 := y^[anza, 2];
maxy2 := y^[anza, 2];
miny3 := y^[anza, 3];
maxy3 := y^[anza, 3];
end;
if (k > anza) then
begin
if (y^[k, 1] < miny1) then
miny1 := y^[k, 1];
if (y^[k, 1] > maxy1) then
maxy1 := y^[k, 1];
if (y^[k, 2] < miny2) then
miny2 := y^[k, 2];
if (y^[k, 2] > maxy2) then
maxy2 := y^[k, 2];
if (y^[k, 3] < miny3) then
miny3 := y^[k, 3];
if (y^[k, 3] > maxy3) then
maxy3 := y^[k, 3];
end;
end;
for k := multiplier to (multiplier + 5) do
begin
if (k = multiplier) then
begin
y^[k, 1] := miny1;
y^[k, 2] := miny2;
y^[k, 3] := miny3;
end;
if (k = multiplier + 1) then
begin
y^[k, 1] := maxy1;
y^[k, 2] := miny2;
y^[k, 3] := miny3;
end;
if (k = multiplier + 2) then
begin
y^[k, 1] := miny1;
y^[k, 2] := miny2;
y^[k, 3] := miny3;
end;
if (k = multiplier + 3) then
begin
y^[k, 1] := miny1;
y^[k, 2] := maxy2;
y^[k, 3] := miny3;
end;
if (k = multiplier + 4) then
begin
y^[k, 1] := miny1;
y^[k, 2] := miny2;
y^[k, 3] := miny3;
end;
if (k = multiplier + 5) then
begin
y^[k, 1] := miny1;
y^[k, 2] := miny2;
y^[k, 3] := maxy3;
end;
end;
scale := (maxy1 - miny1) / (maxy3 - miny3);
scale := 1;
for k := anza to kwisha do
begin
for j := 1 to 2 do
begin
if (graphtype = 2) then
begin
if (algorithm = 1) then
begin
newx^[k, j] := (y^[k, j]) * y^[k, 3] * perspect;
end;
if (algorithm = 2) then
begin
if (j = 1) then
newx^[k, j] := (y^[k, j]) - (perspect * scale * y^[k, 3]);
if (j = 2) then
newx^[k, j] := (y^[k, j]) + (perspect * y^[k, 3]);
end;
end;
if (graphtype = 1) then
begin
newx^[k, j] := y^[k, j];
{write(newx^[k, j] : 8 : 6);}
end;
if ((k = anza) and (j = 2)) then
begin
minx1 := newx^[anza, 1];
maxx1 := newx^[anza, 1];
minx2 := newx^[anza, 2];
maxx2 := newx^[anza, 2];
end;
end;
if (k > anza) then
begin
if (newx^[k, 1] < minx1) then
minx1 := newx^[k, 1];
if (newx^[k, 1] > maxx1) then
maxx1 := newx^[k, 1];
if (newx^[k, 2] < minx2) then
minx2 := newx^[k, 2];
if (newx^[k, 2] > maxx2) then
maxx2 := newx^[k, 2];
end;
end;
moveto(30, 20);
WriteDraw('x1=', miny1 : 5 : 3, ' x2=', miny2 : 5 : 3, ' x3=', miny3 : 5 : 3);
moveto(30, 40);
WriteDraw('x1=', maxy1 : 5 : 3, ' x2=', maxy2 : 5 : 3, ' x3=', maxy3 : 5 : 3);
TextSize(12);
PenSize(2, 2);
for k := anza to kwisha do
begin
for j := 1 to 2 do
begin
transx1 := round(height - (25 + ((newx^[k, 1] - minx1) / (maxx1 - minx1)) * (height - 50)));
transx2 := round(25 + ((newx^[k, 2] - minx2) / (maxx2 - minx2)) * (width - 50));
if ((k = anza) and (j = 1)) then
moveto(transx2, transx1);
if ((k = multiplier) and (j = 1)) then
moveto(transx2, transx1);
lineto(transx2, transx1);
{writeln(transx2, transx1);}
end;
if (traject = 1) then
write('x1=', y^[k, 1] : 8 : 6, ' x2=', y^[k, 2] : 8 : 6, ' x3=', y^[k, 3] : 8 : 6, ' count=', y^[k, 4]);
if ((eqnum = 3) and (traject = 1)) then
writeln;
if ((eqnum = 4) and (traject = 1)) then
writeln(' x4=', y^[k, 4] : 8 : 6);
end;
if (lyap = 2) then
writeln('p1=', p[1], ' p2=', p[2]);
writeln(' Lowx1=', miny1 : 8 : 6, ' Highx1=', maxy1 : 8 : 6, ' Lowx2=', miny2 : 8 : 6, ' Highx2=', maxy2 : 8 : 6, ' Lowx3=', miny3 : 8 : 6, ' Highx3=', maxy3 : 8 : 6);
end;
procedure sumchange;
var
creature: real;
k: integer;
begin
if (ticker3 > 1) then
begin
change := abs(x[item] - varold);
varsum := varsum + change;
{writeln('var for sum=', x[item] : 8 : 6, ' varsum=', varsum : 8 : 6);}
end;
varold := x[item];
for k := 1 to nn do
begin
mx[k] := x[k];
end;
if (nn = 3) then
creature := eq1 * eq2 * eq3;
if (nn = 4) then
creature := eq1 * eq2 * eq3 * eq4;
if ((creature <= 0) and (ticker3 = 1)) then
signal1 := 1;
if ((creature > 0) and (ticker3 = 1)) then
signal1 := 2;
if ((creature < 0) and (signal1 = 2)) then
begin
totsignal := totsignal + 1;
signal1 := 1;
end;
if ((creature > 0) and (signal1 = 1)) then
begin
totsignal := totsignal + 1;
signal1 := 2;
end;
{writeln('creature=', creature, ' totsignal=', totsignal, ' signal1=', signal1, ' ticker3=', ticker3);}
end;
procedure engine;
var
k, j, l: integer;
begin
initcond;
initnums;
while (time <= (iternum + toss)) do
begin
{Initializing the integrator;}
timeone := timeone + h;
ticker3 := ticker3 + 1;
if (kitu = 3) then
special3;
if (kitu = 4) then
special4;
if (kitu = 5) then
special5;
model;
multiplier := round((1 / h) * toss);
if (multiplier > 10000) then
multiplier := 10000;
count := count + 1;
sumchange;
if (count < multiplier) then
begin
for j := 1 to nn do
begin
if (graphtype = 2) then
y^[count, j] := x[j];
if (graphtype = 1) then
begin
y^[count, 1] := x[1];
y^[count, 2] := time;
y^[count, 3] := choice;
end;
end;
y^[count, 4] := count;
end;
for k := 1 to nnn do
begin
x[k] := xnew[k];
end;
if (count = multiplier) then
graph;
{Construct a new orthonormal basis by Gram-Scmidt method;}
{Normalize first vector}
znorm[1] := 0.0;
for j := 1 to nn do
begin
znorm[1] := znorm[1] + sqr(x[nn * j + 1]);
end;
znorm[1] := sqrt(znorm[1]);
for j := 1 to nn do
begin
x[nn * j + 1] := x[nn * j + 1] / znorm[1];
end;
{Generate the new orthonormal set of vectors;}
for j := 2 to nn do
begin
{Generate j-1 coefficients}
for k := 1 to (j - 1) do
begin
gsc[k] := 0.0;
for l := 1 to nn do
begin
gsc[k] := gsc[k] + (x[nn * l + j] * x[nn * l + k]);
end;
end;
{Construct a new vector;}
for k := 1 to nn do
begin
for l := 1 to j - 1 do
begin
x[nn * k + j] := x[nn * k + j] - (gsc[l] * x[nn * k + l]);
end;
end;
{Calculate the vector's norm;}
znorm[j] := 0;
for k := 1 to nn do
begin
znorm[j] := znorm[j] + sqr(x[nn * k + j]);
end;
znorm[j] := sqrt(znorm[j]);
{Normalize the new vector;}
for k := 1 to nn do
begin
x[nn * k + j] := x[nn * k + j] / znorm[j];
end;
end;
{Updating running vector magnitudes;}
if (time > toss) then
begin
for k := 1 to nn do
begin
cum[k] := cum[k] + log2(znorm[k]);
end;
{Normalize exponent and print exponents;}
if ((time > (iternum + toss - h)) and (time <= (iternum + toss))) then
begin
write('time=', time : 8 : 3, ' it #=', ticker2);
for k := 1 to nn do
begin
cumkt := cum[k] / time;
volume := volume + cumkt;
if (graphonly = 2) then
write(' ', cumkt : minwidth);
if (k = 1) then
begin
totderiv1 := cumkt;
{totderiv1 := sdevfour;}
totderiv2 := (varsum * totsignal) / time;
totderiv3 := (sdevfour * totpower);
end;
if (k > 1) then
begin
if ((cumkt > totderiv1) and (cumkt < -0.05)) then
totderiv1 := cumkt;
end;
end;
if (volume > 0) then
alert := alert + 1;
if (graphonly = 2) then
writeln(' vol=', volume);
writeln;
writeln('turbulence 2=', totderiv2 : minwidth, ' varsum=', varsum, ' totsignal=', totsignal);
writeln('turbulence 3=', totderiv3 : minwidth);
end;
end;
time := time + h;
end;
end;
procedure answer;
var
parnum, anza1, ii: integer;
begin
write('What model do you want? 1-pendulum 2-Lorenz 3-Ecosphere 4-Maharishi 5-Nazi ');
readln(input, kitu);
{set up number of equations}
nn := 3;
nnn := 12;
eqnum := 3;
if (kitu = 3) then
begin
write('Do you want the four or three equation model? 3-three 4-four ');
readln(input, eqnum);
end;
if (eqnum = 4) then
begin
nn := 4;
nnn := 20;
politics := 1;
end;
write('Do you want one analysis or an image? 1-One analysis 2-image ');
readln(input, lyap);
write('For what variable do you want to sum the change for turbulence? ');
readln(input, item);
write('Do you want to suppress the first term of the power series? 1-yes 2-no ');
readln(input, suppress);
write('Do you want to view the variable values? 1-yes 2-no ');
readln(input, traject);
write('Do you want Lyaponov exponents? 1-graph only 2-Lyaponov exponents ');
readln(input, graphonly);
if (graphonly = 1) then
nnn := 3;
if ((graphonly = 1) and (kitu = 3) and (eqnum = 4)) then
nnn := 4;
write('What type of graph do you want? 1-Overtime 2-Phase ');
readln(input, graphtype);
perspect := 1;
if (graphtype = 2) then
begin
write('What type of algorithm do you want to draw the perspective? 1-multiplicative 2-additive ');
readln(input, algorithm);
write('Do you want a 2D or a 3D phase diagram? 0=2D 1=3D ');
readln(input, perspect);
end;
if (kitu = 3) then
begin
write('How many years of lag do you want for the public pollution reaction? ');
readln(input, yearslag);
write('How many years half-life do you want for the current pollution dissipation? ');
readln(input, halflife);
decay := -ln(0.5) / halflife;
writeln('The decay constant in the model is: ', decay);
write('Do you want shocks? 1-yes 2-no ');
readln(input, shocks);
shock := 0;
if (shocks = 1) then
begin
write('What size shock? ');
readln(input, shock);
end;
if (eqnum = 3) then
begin
write('Do you want political policy differences? 1-yes 2-no ');
readln(input, politics);
end;
end;
if (kitu = 4) then
begin
write('How many time periods of lag do you want for the conflict response? ');
readln(input, yearslag);
end;
write('How many time periods do you want for the graph? ');
readln(input, toss);
write('Do you want to lag the start of the graph? 1-yes 2-no ');
readln(input, sasa);
if (sasa = 1) then
begin
write('At what time period do you want the graph to begin? ');
readln(input, anza1);
end;
iternum := 1;
if (graphonly = 2) then
begin
write('How many more time periods do you want for Lyaponov analysis? ');
readln(input, iternum);
end;
write('Step size? ');
readln(input, h);
if (sasa = 1) then
anza := round(anza1 / h);
if (sasa = 2) then
anza := 1;
if (kitu <> 5) then
begin
if (lyap = 1) then
begin
write('How many parameters? ');
readln(input, parnum);
for ii := 1 to parnum do
begin
write('What is the value for parameter ', ii, ' ?');
readln(input, p[ii]);
end;
end;
end;
if (kitu = 5) then
parms;
if (lyap = 2) then
begin
writeln('You can vary your first two parameters between ranges.');
write('What is your lower value in the range for your first parameter? ');
readln(input, lowp1);
write('What is your higher value in the range for your first parameter? ');
readln(input, highp1);
write('What is your lower value in the range for your second parameter? ');
readln(input, lowp2);
write('What is your higher value in the range for your second parameter? ');
readln(input, highp2);
write('How many other parameters are there in the model? ');
readln(input, parnum);
for ii := 3 to (parnum + 2) do
begin
write('What is the value for parameter ', ii, ' ?');
readln(input, p[ii]);
end;
write('How large do you want the image matrix? ');
readln(input, dim);
end;
lagg := round(yearslag / h);
if (yearslag = 0) then
lagg := 1;
end;
procedure minmax;
begin
ticker := ticker + 1;
if (ticker = 1) then
begin
min := totderiv1;
max := totderiv1;
min2 := totderiv2;
max2 := totderiv2;
min3 := totderiv3;
max3 := totderiv3;
end;
if totderiv1 < min then
min := totderiv1;
if totderiv1 > max then
max := totderiv1;
if totderiv2 < min2 then
min2 := totderiv2;
if totderiv2 > max2 then
max2 := totderiv2;
if totderiv3 < min3 then
min3 := totderiv3;
if totderiv3 > max3 then
max3 := totderiv3;
end;
procedure axes;
var
d1, d2: integer;
begin
for d2 := 1 to 2 do
begin
for d1 := 1 to dim do
begin
write(ff1, ' ', d1);
write(ff2, ' ', d1);
write(ff3, ' ', d1);
end;
writeln(ff1);
writeln(ff2);
writeln(ff3);
end;
end;
procedure matfill;
var
c1, c2: integer;
begin
ticker := 0;
for c2 := dim downto 1 do
begin
shift := lowp2 + ((highp2 - lowp2) * ((c2 - 1) / (dim - 1)));
p[2] := shift;
for c1 := 1 to dim do
begin
ticker2 := ticker2 + 1;
p[1] := lowp1 + ((highp1 - lowp1) * ((c1 - 1) / (dim - 1)));
engine;
minmax;
image1^[(dim - c2 + 1), c1] := totderiv1;
image2^[(dim - c2 + 1), c1] := totderiv2;
image3^[(dim - c2 + 1), c1] := totderiv3;
end;
end;
end;
procedure printit;
var
c1, c2: integer;
begin
writeln(ff1, dim, ' ', dim);
writeln(ff2, dim, ' ', dim);
writeln(ff3, dim, ' ', dim);
writeln('dimension=', dim, ' by ', dim);
writeln(ff1, max : minwidth, ' ', min : minwidth);
writeln(ff2, max2 : minwidth, ' ', min2 : minwidth);
writeln(ff3, max3 : minwidth, ' ', min3 : minwidth);
writeln('max=', max : minwidth, ' min=', min : minwidth);
writeln('max=', max2 : minwidth, ' min=', min2 : minwidth);
writeln('max=', max3 : minwidth, ' min=', min3 : minwidth);
writeln;
axes;
for c2 := dim downto 1 do
begin
for c1 := 1 to dim do
begin
write(ff1, ' ', image1^[(dim - c2 + 1), c1] : minwidth);
write(ff2, ' ', image2^[(dim - c2 + 1), c1] : minwidth);
write(ff3, ' ', image3^[(dim - c2 + 1), c1] : minwidth);
if c1 < 6 then
write(' ', image1^[(dim - c2 + 1), c1] : minwidth);
end;
writeln(ff1);
writeln(ff2);
writeln(ff3);
writeln;
end;
writeln;
for c2 := dim downto 1 do
begin
for c1 := 1 to dim do
begin
if c1 < 6 then
write(' ', image2^[(dim - c2 + 1), c1] : minwidth);
end;
writeln;
end;
writeln;
for c2 := dim downto 1 do
begin
for c1 := 1 to dim do
begin
if c1 < 6 then
write(' ', image3^[(dim - c2 + 1), c1] : minwidth);
end;
writeln;
end;
end;
begin
initialize;
showtext;
showdrawing;
answer;
if (lyap = 1) then
engine;
if (lyap = 2) then
begin
rewrite(ff1, f1name);
rewrite(ff2, f2name);
rewrite(ff3, f3name);
new(image1);
new(image2);
new(image3);
matfill;
printit;
close(ff1);
close(ff2);
close(ff3);
writeln('lpwp1=', lowp1, ' highp1=', highp1, ' lowp2=', lowp2, ' highp2=', highp2, ' p3=', p[3]);
end;
writeln('A nondisipative system was detected', alert, ' times.');
end.