|
|
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.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|