Can't understand why numbers are decreasing?

[jon93]

Hello everyone,
I am very puzzled and need to have my sanity restored.
I am writing a trivial program in which a ball (represented by a pset point) bounces inside a box. (Planning to use it for more complicated things)
There is gravity (g) which pulls the ball down. (increasing y direction values )
There is a potentially infinite repetition of calculation of new y velocity (yvel) (increased by gravity) changing the y position, and x velocity (xvel) being constant. 
When it hits the side it bounces in the other direction.
When it hits the bottom it bounces back up.
It should just go on forever bouncing around.
But for some reason the height of the bounce (y position value) gradually decreases. The y velocity (yvel) gradually decreases.

I cannot see why this happens.

Is it a result of a cumulative effect of rounding numbers when there are many iterations? 



global Xmax, Ymax, Pi, Deg, Rad

global g, iters

global xpos, ypos, previousypos, previousypos, xvel, yvel, itersToDo

global previousxvel, previousyvel

global minxpos, maxxpos, minypos, maxypos, iterstodo

global xrange 'difference between x left of box and x right of box

global xshow , yshow

global slownum

global howoftentocls

global multiplier 'i used this to change the values of xpos etc that are being calculated

'then divide by it to change it into values that can be displayed

'it serves no purpose except to see if the problem with the program

'is associated with the size of the numbers being calculated

'are they for 0 to 2 or 0 to 100 to 100 etc

'as far as I can see the problme has nothing to do with the size of the number



ulx = 10 'upper left x value of box on screen

uly = 10

ww =1200

wh =600

WindowWidth = ww

WindowHeight = wh


UpperLeftX = ulx

UpperLeftY = uly



open "Graphic Title" for graphics_nsb_nf as #gr '<======================= title

#gr "setfocus"

#gr "trapclose [quit]"

#gr "when leftButtonUp lButtonUp"

#gr "when characterInput charIn"

#gr "down"

#gr "fill black"



'#gr "trapclose [quit]"

'wait



call back 255,0,255

call fore 255,0,0

b$ = ""



call StartVals

for iters = 1 to iterstodo

call Sloww 'slow things down so it can be seen

call SavePreviousXY

call MoveBall

call BounceSub

call pset xpos/multiplier, ypos/multiplier

'call circ xpos/distancemultipiernum, ypos/distancemultipiernum, 2



'sometimes need to clear display

if iters mod howoftentocls = 1 then

print #gr, "fill black"

end if

next





'============================== sets drawing

#gr "flush"

wait

[quit]

close #gr

end

'***********************************************



sub StartVals

multiplier = 1

slownum= 100 'slow things down

howoftentocls = 3000 'how often to clear the display

iterstodo = 19500

g = .03 * multiplier';gravity constant

'min and max x and y positions. This determines where the ball bounces.

minxpos = 0

maxxpos = 1150 * multiplier',1

minypos = 0

maxypos = 555 * multiplier'1

xpos = 20 * multiplier'0 .200

ypos = 20 * multiplier'0' .200

yvel = 0

xvel =.09 * multiplier

end sub



sub BounceSub

'bounce off walls and bottom

'bounce off left or right wall

if xpos <= minxpos and previousxpos > minxpos then

xvel = xvel * -1



else

if xpos >= maxxpos and previousxpos < maxxpos then

xvel = xvel * -1

end if

end if

'bounce off bottom

if ypos >= maxypos and previousypos < maxypos then

print "now in bounce sub iters",iters

print "previous xvel yvel", previousxvel, previousyvel

print "currenr xvel yvel", xvel, yvel

yvel = yvel * -1

end if

end sub





sub MoveBall

yvel= yvel + g

xpos= xpos + xvel

ypos= ypos + yvel

end sub



sub SavePreviousXY 'probably not actually needed

previousxvel = xvel

previousyvel = yvel

previousxpos = xpos

previousypos = ypos

end sub



sub Sloww' slow things down

for j = 1 to slownum

dummy=1

next

end sub



sub fore r, g, b

#gr "color ";r;" ";g;" ";b

end sub



sub back r, g, b 'backcolor is used for fills

#gr "backcolor ";r;" ";g;" ";b

end sub



sub pset x, y

#gr "set ";x;" ";y

end sub



sub circ x, y, radius

#gr "place ";x;" ";y;"; circle ";radius

end sub



'Need line: #gr "trapclose quit"

sub quit H$

close #gr '<=== this needs Global H$ = "gr"

end 'Thanks Facundo, close graphic wo error

end sub

[plus]

I would think the gravity would take away from the rebound and eventually the ball would cease to bounce.
yeah right here!

sub MoveBall
yvel= yvel + g
xpos= xpos + xvel
ypos= ypos + yvel
end sub


This code template looks familiar! ;-))

from
global Xmax, Ymax, Pi, Deg, Rad

to
'Need line: #gr "trapclose quit
sub quit H$
close #gr '<=== this needs Global H$ = "gr"
end 'Thanks Facundo, close graphic wo error
end sub

[tsh73]

It has no
previousxpos
in line
global xpos, ypos, previousypos, previousypos, xvel, yvel, itersToDo
so it would not bounce from left side.

As for going down, I think it's because you use most simple method (Euler's)
Next order method would be better (some R-K second order method)
(it would still accumulate errors but they would not be visible)
But I have problems applying them to formulas without visible time step.

For some weird reason, two changes seem to fix "doing down" thing
1) add line
ypos = maxypos
then bouncing from bottom,
Reason that it start from floor level each time, not beneath
- but this one is not enough
2) try different values of multiplier
It should not matter - but some values (0.01, 0.02, 0.04, 27) make upper value stay constant
(and 0.03 goes down for a while then stays constant)
I really don't have idea why.

[jon93]

Problem solved.
tsh73 suggested adding the line
ypos = maxypos
This sort of worked.
Then I realised that when I reversed the yvel when bouncing at the bottom, the velocity I was dealing with was the velocity at the PREVIOUS position, but the ball is lower down at the PRESENT ypos.
So now the ball was lower down (nearer bottom of screen) but only travelling at the velocity it had before crossing the line.
It lacked the acceleration it would have had in that last iteration.
So it could not get high enough at the apex.

So I added
ypos = previousypos
before changing the sign on the yvel.
works fine.
goes on forever without drooping.

Thank you for you help.

[tenochtitlanuk]

Just a quick example of using Runge-Kutta mentioned by Anatoly. ( Based on Wikipediam and the Rosetta Code problem)

ALL methods that step forward by finite steps, following a 'recipe' about values and their changes, are likely to accumulate errors.

The code below steps forward by smaller and smaller steps, and I have stopped it when the steps have become small enough to ALMOST match the known solution. Just keep clicking the 'notice' that pops up.



Available as RK.zip which includes the needed bmp file.


     nomainwin

     WindowWidth  = 750
     WindowHeight = 700

     open "Runge-Kutta and exact solution ( with diff't dt steps) compared" for graphics_nsb as #w

     #w "trapclose quit"

     #w "down ; fill darkblue ; size 1 ; flush"
     #w "up ; goto 0 650 ; down ; color white ; goto 1000 650"

     for t =0 to 10
       #w "up ; goto "; 80 *t; " 650 ;   down ; color white ; goto ";      80 *t; " 0"
     next t

     for y =0 to 700 step 100
       #w "up ; goto 0 "; 650 -y; " ; down ; color white ; goto 750 "; 650 -y
     next y

     loadbmp "scr", "RK.bmp"
     #w "drawbmp scr 20 4"

     y   =1.0

     dt  =2.0

[start]
     #w "up ; goto 0 650 ; down"

     for t = 0 to 10 step dt

       dy1  =dt *( t        ) *sqr( y             )
       dy2  =dt *( t + dt /2) *sqr( y +dt /2 *dy1 )
       dy3  =dt *( t + dt /2) *sqr( y +dt /2 *dy2 )
       dy4  =dt *( t + dt   ) *sqr( y +dt    *dy3 )

       y   = y   +( dy1 +2 *dy2 +2 *dy3 + dy4) /6

       #w "color cyan   ; goto "; int( 80 *t); " "; int( 650 -y)
     next t

     #w "flush"

     notice "Added with step "; using( "#.####", dt)
     dt =dt /1.2
     if dt >=0.1 then y =1: goto [start]

     notice "Adding known solution "

     #w "up ; goto 0 650 ; down ; size 2"

     for t = 0 to 10 step dt

       dy1  =dt *( t        ) *sqr( y             )
       dy2  =dt *( t + dt /2) *sqr( y +dt /2 *dy1 )
       dy3  =dt *( t + dt /2) *sqr( y +dt /2 *dy2 )
       dy4  =dt *( t + dt   ) *sqr( y +dt    *dy3 )

       y   = y   +( dy1 +2 *dy2 +2 *dy3 + dy4) /6

       #w "color    red ; goto "; int( 80 *t); " "; int( 650 -((t^2 + 4)^2) / 16)
     next t

   wait

   sub quit h$
       close #w
       end
   end sub


[tsh73]

Glad you fixed it.

Please use code tags (button with C in rectangle on post editor toolbar) when posting code.
This way code will be in a block with a scrollbar, easier to go by to next message in a thread.
(also extra empty lines after each line of code is prevented by posting in "BBCode" editor tab, not in "Previw" tab.
Tabs are visible below post editor)