• Print

Author Topic: Stability equals reliability of QB64 computations  (Read 170 times)

Mrwhy

  • Hero Member
  • *****
  • Posts: 2890
  • My Dad called me Mr Why when I was 5.
    • Email
Stability equals reliability of QB64 computations
« on: March 24, 2013, 09:55:07 AM »
Not all situations are either stable or not.

Here is one that starts looking very stable (for 1 million "hours") but then unexpectedly and suddenly at around 1.6 million "hours" diverges to extremes (but stays on screen).
But then, to our amazement it tries to return to its "stable state".
This it has almost managed by 1.85 million hours, but thereafter suffers a relapse before RETURNING evey so often (around 800 thousand hours!) to and from the region of its stable state.

So results such as these MIGHT be reliable or BECOME reliable again every 0.8 million hours or so ::)
 
Code: [Select]
DEFDBL A-H, J-M, O-Z: DEFINT I, N
n = 12: OPTION BASE 1
DIM x(n), xc(n), f(n), c1(n), c2(n), c3(n), c4(n)
GOSUB graphics
m1 = 1#: m2 = 1.0001#: m3 = .9999#
x3 = -.2: y3 = -.2
'LOCATE 13, 1: INPUT "x1"; x1
t = 0#
' ** initial coordinates **
xc(1) = 1: xc(2) = -.5: xc(3) = -.5 'x coords
xc(4) = 0#: xc(5) = .3215#: xc(6) = -.3215# 'y coords
' ** initial velocities **
xc(7) = 0: xc(8) = -1.135#: xc(9) = -xc(8) 'x vel
xc(10) = -.4865#: xc(11) = -.5 * xc(10): xc(12) = -.5 * xc(10) 'y vel
12 hscale = .000063#
agn:
j = j + 1: LOCATE 1, 1: PRINT j
GOSUB Runge
d12 = (xc(1) - xc(2)) ^ 2 + (xc(4) - xc(5)) ^ 2
d13 = (xc(1) - xc(3)) ^ 2 + (xc(4) - xc(6)) ^ 2
lit = I / d12 + 1 / d13
'IF j < 123600 THEN PSET (j / 123600 - 1, lit / 15 - 1.2)
'PRINT lit-4
PSET (xc(1), xc(4)), 12
PSET (xc(2), xc(5)), 9
PSET (xc(3), xc(6)), 14
t = t + h
'IF j = 100000 THEN CLS: j = 0: GOSUB energy
h = hscale * (r12 * r23 * r31)
GOTO agn
END

rotn:
xo = x: x = x - wd * y: y = y + wd * xo
RETURN

Equations:
d21 = x(2) - x(1): d32 = x(3) - x(2): d13 = x(1) - x(3)
d54 = x(5) - x(4): d65 = x(6) - x(5): d46 = x(4) - x(6)
r12 = (d21 ^ 2# + d54 ^ 2#) ^ .5#
r23 = (d32 ^ 2# + d65 ^ 2#) ^ .5#
r31 = (d13 ^ 2# + d46 ^ 2#) ^ .5#
p12 = r12 ^ 3#: p23 = r23 ^ 3#: p31 = r31 ^ 3#
f(1) = x(7): f(2) = x(8): f(3) = x(9)
f(4) = x(10): f(5) = x(11): f(6) = x(12)
f(7) = m2 * d21 / p12 - m3 * d13 / p31
f(8) = m3 * d32 / p23 - m1 * d21 / p12
f(9) = m1 * d13 / p31 - m2 * d32 / p23
f(10) = m2 * d54 / p12 - m3 * d46 / p31
f(11) = m3 * d65 / p23 - m1 * d54 / p12
f(12) = m1 * d46 / p31 - m2 * d65 / p23
RETURN

Runge:
FOR I = 1 TO n: x(I) = xc(I): NEXT
GOSUB Equations
FOR I = 1 TO n: c1(I) = h * f(I): NEXT
FOR I = 1 TO n: x(I) = xc(I) + c1(I) / 2#: NEXT
GOSUB Equations
FOR I = 1 TO n: c2(I) = h * f(I): NEXT
FOR I = 1 TO n: x(I) = xc(I) + c2(I) / 2#: NEXT
GOSUB Equations
FOR I = 1 TO n: c3(I) = h * f(I): NEXT
FOR I = 1 TO n: x(I) = xc(I) + c3(I): NEXT
GOSUB Equations
FOR I = 1 TO n: c4(I) = h * f(I): NEXT
FOR I = 1 TO n
    xc(I) = xc(I) + (c1(I) + 2# * c2(I) + 2# * c3(I) + c4(I)) / 6#
NEXT
RETURN

graphics:
SCREEN 12
'PAINT (1, 1), 9
xm = 1: ym = 1
'VIEW (180, 17)-(595, 460), 0, 13
WINDOW (-xm, -ym)-(xm, ym)
'LINE (-xm, 0)-(xm, 0), 8: LINE (0, -ym)-(0, ym), 8
'LOCATE 15, 76: PRINT xm
RETURN

energy:
' ** Is energy conserved? **
kin1 = .5# * m1 * (xc(7) ^ 2# + xc(10) ^ 2#)
kin2 = .5# * m2 * (xc(8) ^ 2# + xc(11) ^ 2#)
kin3 = .5# * m3 * (xc(9) ^ 2# + xc(12) ^ 2#)
pot = -(m1 * m2 / r12 + m2 * m3 / r23 + m3 * m1 / r31) '
energy = kin1 + kin2 + kin3 + pot
LOCATE 21, 1: PRINT "Energy"
LOCATE 22, 1: PRINT energy
LOCATE 17, 7: PRINT "             "
LOCATE 17, 1: PRINT "Time ="; CSNG(t)
RETURN


bobtheunplayer

  • Jr. Member
  • **
  • Posts: 70
  • I'd rather be coding.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #1 on: April 05, 2013, 01:52:29 AM »
Quote
Not all situations are either stable or not.
Actually, all situations are stable or not.  ;)

The interestingly calculated Lamniscate requires infinite precision in a finite system; hence, the appearance of instability then stability as values overflow and reset.  Try a fixed approach.  Generate a graph using a granular loop to minimize gaps as sin / cos of theta approach zero, then save only the unique coordinates.  Lastly, loop through the array of pre-calculated values to draw "infinity" for infinity without devating from the curve.

Code: [Select]
' file: polar.bas
' date: 4-Apr-2013
' auth: bsmith
' info: a graphed Lemniscate that chases its tail

'
' *** initialize ***
SCREEN 12
_TITLE "Lamniscate Snake"

TYPE coordType
    x AS INTEGER
    y AS INTEGER
END TYPE

OPTION BASE 0
DIM coords(3000) AS coordType
DIM offset AS coordType
DIM deg2rad AS DOUBLE
DIM radius AS DOUBLE
DIM radiusSquared AS DOUBLE
DIM theta AS DOUBLE
DIM index AS INTEGER
DIM head AS INTEGER
DIM tail AS INTEGER
DIM lead AS INTEGER
DIM leadCount AS INTEGER

magnitude = 300
deg2rad = 3.14159265359 / 180
offset.x = 320
offset.y = 240
index = 1
coords(index - 1).x = -1
coords(index - 1).y = -1

' pre-render a Lemniscate curve
' radiusý = magnitudeý * cos(2é)

' this loop iterates 361,000 times.  testing shows that only about 2000 coordinates are produced
' the granular iterations are important to capture the complete curve as sin é and cos é approach zero
FOR degrees = 0 TO 360 STEP .001
    ' calculate (r,é)
    ' radiusSquared cannot be negative, this curve ignores imaginary numbers
    radius = 0
    theta = degrees * deg2rad
    radiusSquared = (magnitude ^ 2) * COS(2 * theta)
    IF radiusSquared > 0 THEN radius = SQR(radiusSquared)

    ' convert polar to cartesian, add offsets
    coords(index).x = INT(radius * COS(theta)) + offset.x
    coords(index).y = INT(radius * SIN(theta)) + offset.y

    ' manipulate the quadrant for a proper figure eight path
    IF degrees > 90 AND degrees < 270 THEN coords(index).y = ((coords(index).y - offset.y) * -1) + offset.y

    ' save only unique coordinates
    IF (coords(index - 1).x <> coords(index).x) OR (coords(index - 1).y <> coords(index).y) THEN
        index = index + 1
    END IF
NEXT degrees

' setup "snake"
head = 0
tail = 0
lead = 900
leadCount = 0
key$ = ""


'
' *** run program ***
'

' instructions
PRINT "Esc to quit."

' commence tail chasing
DO
    head = intRotate(1, index, head, 1)
    PSET (coords(head).x, coords(head).y)

    IF leadCount > lead THEN
        tail = intRotate(1, index, tail, 1)
        PSET (coords(tail).x, coords(tail).y), 0
    ELSE
        leadCount = leadCount + 1
    END IF

    _DELAY .001
    key$ = INKEY$
LOOP UNTIL (key$ = CHR$(27))
END

'
' *** subs ***
'
FUNCTION intRotate (iFloor AS INTEGER, iCeil AS INTEGER, iValue AS INTEGER, iIncrament AS INTEGER)
' rotates an integer
' "i" prefix used to avoid namespace collision
iValue = iValue + iIncrament
IF iValue > iCeil THEN
    iDiff = iValue - iCeil
    iValue = iFloor + iDiff
END IF
IF iValue < iFloor THEN
    iDiff = iFloor - iValue
    iValue = iCeil - iDiff
END IF
intRotate = iValue
END FUNCTION
« Last Edit: April 05, 2013, 02:04:31 AM by bobtheunplayer »

OlDosLover

  • Hero Member
  • *****
  • Posts: 3859
  • OlDosLover
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #2 on: April 05, 2013, 01:57:14 AM »
Hi all,
    Very nice graphic demo bobtheunplayer. Looks pretty perfect to me. Not that i understand maths very well! A really good job , thanks for sharing!
OlDosLover.

bobtheunplayer

  • Jr. Member
  • **
  • Posts: 70
  • I'd rather be coding.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #3 on: April 05, 2013, 02:10:16 AM »
Thanks OlDosLover.  The math isn't so important.  I was hoping the take-away would be the idea of using a pre-calculated set of values to achieve predictable results.

It does look kind of mesmerizing though.  Been running it on my second monitor for almost an hour for giggles.  ;D

Mrwhy

  • Hero Member
  • *****
  • Posts: 2890
  • My Dad called me Mr Why when I was 5.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #4 on: April 05, 2013, 02:56:59 AM »
Quote from: bobtheunplayer on April 05, 2013, 01:52:29 AM
Quote
Actually, all situations are stable or not. 

Hello Bob
It is great to meet someone who really enjoys thinking and what you say is very interesting.

So you believe thet stable/unstable is a "black or white" situation with no  "shades of grey".
I wonder why you feel that way?

When we attempt, by computer, to follow a path we are bound to make small errors and if each next computation reliies on earliier results  then these errors may accumulate without limit. Or they may give an oscillation about the "true path".
The extent of this oscillation in most "stable" situations is small enough that the computed path does not completely change its shape. So the insight provided is valid.

But in the example I gave, it does while in the example you gave it doesn't

bobtheunplayer

  • Jr. Member
  • **
  • Posts: 70
  • I'd rather be coding.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #5 on: April 05, 2013, 03:17:49 AM »
Quote
So you believe thet stable/unstable is a "black or white" situation with no  "shades of grey".
I wonder why you feel that way?

Not trying to go down the philisophical path.  Just speaking of computer systems only.


The path will never deviate if the computations produced results where none of the precision is lost from computation to computation.  Given the chaos of the calculations provided in the original program, the loss of precision is unavoidable which leads to variations in the path.  It may be possible to contain the path by using limits and truncating values.

Mrwhy

  • Hero Member
  • *****
  • Posts: 2890
  • My Dad called me Mr Why when I was 5.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #6 on: April 05, 2013, 04:57:10 AM »
Yes, you are dead right
It MAY not be possible.
In the case I gave, the calculations are neither exact nor totally useless, but give an INTERMEDIATE result - the lemoniscate curve is NOT closed, not periodic and not-at-all the same curve as the one that for a long time DID seem stable. It eventually and slowly changes beyond recognition and, amazingly, thereafter "tries to return" to its former shape more than once - but each is a mere distant cousin.

The locus given in my program is that of 3 bodies moving under their mutual gravity.
When their masses are the same they can follow a closed curve that stays closed (periodic) no matter how long we continue the calculation.
When one of them is 0.1% too heavy, they cannot and do not.

All I am trying to say is that sometimes computations can give useful results and it is vital we think whether the current computation is one of those!
If you know any way to tell, please let us know ;D

bobtheunplayer

  • Jr. Member
  • **
  • Posts: 70
  • I'd rather be coding.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #7 on: April 05, 2013, 05:56:55 AM »
I'm going to go with that it's not possible, because of the limited resources on a computer, so I would suggest using limiters that cause the enitity to "bounce" the other direction when hit.  The best you can hope for on a computer is "close enough", but given that a computer can offer many digits of precision, "close enough" should fit most realistic situations.

Mrwhy

  • Hero Member
  • *****
  • Posts: 2890
  • My Dad called me Mr Why when I was 5.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #8 on: April 05, 2013, 06:45:02 AM »
The reason computers offer such great "precision" is nobody has any idea how much is required.
If you are going to divide something by x-y, which differ by less than 1% (in a world where 1st class precision of electrical instruments is 1% of full scale deflection), then I guess you get what you deserve ;D
The problem is you get a "result" but it is not possible to verify it experimentally.
Thus one of the principal uses of computers - to simulate reality - fails.

My calculation was not "To get a preordained (="guessed") curve". My aim was to discover how 3 such heavenly bodies would move under Newton's rule-of-thumb (gravity)!


As for computed loci, there seem to be four types
1. Those that go off the screen, no matter how big it is
2. Those that dwindle to a single point
3. Those that repeat periodically (closed loops)
4. Those that never repeat but for limited screen resolution (i.e. after they have visited all the points of a region at least once)

I'd like to know your thoughts on each of those.

Humans like to think of the world in terms of "events"
A thing like a bonfire neither repeats, nor is stable nor unstable. Thus it is an event - a hapening.

bobtheunplayer

  • Jr. Member
  • **
  • Posts: 70
  • I'd rather be coding.
    • Email
Re: Stability equals reliability of QB64 computations
« Reply #9 on: April 05, 2013, 09:02:30 AM »
I don't really have an opion on such matters.  Never took the time to think about them.

  • Print