|
|
The following program was used to explore the highly nonlinear and potentially
chaotic
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.
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
