If we split the time interval (0,T) into N segments of length l =T/N
and iterate
xnew = xold + lvold ; vnew = vold + la.
or more formally
xn+1 = xn + lvn ; vn+1 = vn + la
N times what will we get? This is an example of a recurrance relation.
Recurrance relations describe a sequence of terms by specifying a given
term as an expression involving previous terms in the sequence. The Fibonacci numbers, for example,
1,1,2,3,5,8,13,... are uniquely specified by the recurrance relation
xn+1 = xn + xn-1 together with "initial condition" x1 = x0 = 1.
The advantage of recurrance relations for the programmer is
that the position of an object at time t=nl, though theoretically given by
a quadratic in nl, can be derived without multiplication from xn-1 (its
position last game cycle) and the quantity lv which (if l is sensibly
chosen, often l = 1) is easily calculated.
The movement iteration above in fact gives
xn = x0 + nlu + n(n-1)l2a/2 ;
vn = v0 + nla
since vn is readily solved and then xn follows. We see from this that xn is not the correct position of the body
at time t=nl.
The correct value is x0 + nlu + n2l2a/2 but this differs by
only nl2a/2 which is usually negligible.
If the position is required to be exactly correct then there are two ways to achieve this:
(i) Iterate as above but start with v0 = u + la/2.
(ii) Use the following slightly modified iteration:
vtemp = vold + la/2
xnew = xold + lvtemp
vnew = vtemp + la/2
These iterations extend readily to two or three dimensional motion,
the position, velocity, and acceleartion terms becoming vectors.
Solving Recurrance Relations
To solve a recurrance relation of the form
g(xn+k,...,xn+1,xn) = f(n)
for some linear function g and arbitary function f we find the general
solution to g(xn+k,...,xn+1,xn) = 0 and then add to it any specific
solution to g(xn+k,...,xn+1,xn) = f(n), which is often a function of
similar form to f(n).
The general solution to axn+1 + bxn = 0 (a¹0) is xn= A (-b/a)n.
The form of the general solution to
axn+2 + bxn+1 +cxn = 0 (a¹0)
depends on the nature of the roots of the auxiliary equation am2+bm+c=0.
If the roots are repeated (b2=4ac) then the general solution is
xn = (A+Bn)mn where m is the repeated root -b/2a.
If the roots are real and distinct (b2>4ac) then the general solution
is xn=Am1n + Bm2n where m1 and m2 are the two roots (-b ± Ö (b2-4ac))/2a.
If the roots are imaginary (b2<4ac) then the general solution is
xn = (c/a)n/2(Acos(nq)+Bsin(nq)) or equivalently xn = (c/a)n/2Acos(nq+e) where
tanq = - Ö(4ac-b2)/b.
Example - Damped Motion
It is often desirable to damp the motion of an object each game cycle
so that objects come to rest if not explicitly accelerated. We will
consider motion in N dimensions and write r for the position vector of the
point under consideration at time t, v for its velocity vector (the first
derivative of r with respect to time) at time t, and a for its constant
applied acceleration vector.
We might wish to simulate motion satisfying dv/dt = -sv+a
where s is the damping term.
This differential equation has solution
r = r0 + (1-e-st)(v0-a/s)/s + ta/s
v = (v0-a/s)e-st + a/s
for constant a and we see that in a unit time interval v is decreased by a
damping factor of e-s. As t®¥ these approximate
r = ta/s + (r0 + (v0-a/s))
v = a/s
We would typically simulate this with the iteration
rn+1 = rn + lvn
vn+1 = dvn + la
where l is time step simulated by one iteration cycle and
d = e-sl is the
damping factor for one iteration period.
We will consider here the more general iteration
rn+1 = rn + l1vn
vn+1 = dvn + l2a.
(by taking l1=l, l2=dl, this covers the alternative iteration
rn+1 = rn + lvn  ;  
vn+1 = d(vn + la) ).
To find the explicit forms for this iteration we first find vn as the
solution of vn+1 - dvn = l2a.
The general solution is thus Adn where A is an
arbitary vector. The form of f(n) is a constant vector so we look for a
constant specific solution.
That is, we try to solve x-dx = l2a and
obtain x=l2a/(1-d).
Our full solution is therefore
vn = Adn + l2a/(1-d) and
it remains only to satisfy the initial conditions, ie. to find a value for A that gives
the desired value for v0 .
Once vn is solved we use it to find rn using results A1.1 and A1.9.
The results are
l1(1-dn) | l1l2 | æ | (1-dn) | ö | ||||||||
rn | = | r0 | + | ¾¾¾¾ | v0 | + | ¾¾ | ç | n- | ¾¾ | ÷ | a |
(1-d) | (1-d) | è | (1-d) | ø |
(1-dn) | ||||||
vn | = | dnv0 | + | l2 | ¾¾ | a |
(1-d) |
l1l2 | æ | l2 | l1l2 | ö | ||||||||||
rn | = | n | ¾¾ | a | + | ç | r0 | + | ¾¾ | v0 | - | ¾¾ | a | ÷ |
(1-d) | è | (1-d) | (1-d) | ø |
l2 | |||
vn | = | ¾¾ | a |
(1-d) |
n | x3n | x2n | x1n | x0n |
0 | d | al2+ bl+c | 6al+2b | 6a |
1 | al3+ bl2+ cl+d | 7al2+3bl +c | 12al +2b | |
2 | 8al3+4bl2+2cl+d | 19al2+5bl+c | ||
3 | 27al3+9bl2+3cl+d |
Each stage of the iterative process requires M additions of the form
x=y+lz. If l is a "friendly" value such as 2-p this is okay
but if l is simply 1/N for some general N we would rather not incorporate it into the
iteration if at all possible. Suppose we wish to evaluate the
Xcoordinate of a bezier curve cubic x(t)=at3+bt2+ct+d
for t=n/N " n Î ZN.
To do this we iterate
x0n' = x00 " n ³ 0
x1n' = x1(n-1)' + x00' " n ³ 1
x2n' = x2(n-1)' + x1(n-1)' " n ³ 1
x3n' = x3(n-1)' + x2(n-1)' " n ³ 1
which will work if
x3n' = x3n ;
x2n' = lx2n ;
x1n' = l2x1n ;
x0n' = l3x0n " n ³ 1
so we start from the initial values
x00' = l3x00,
x10' = l2x10,
x20' = lx20,
x30' = x30 where l =1/N .
The price of loosing l from the iteration is that the initial count values become very small for small l but any inaccuracies in them will ripple up into increasing inaccuracies in x3n'.
This method extends readily to polynomials of higher degree and can
be used to fit polynomials to given data.
Succesive Evaluations Via The Midpoint Algorithm
Suppose we wish to iterate a function u=U(x) for x=0,1,2,....
and that the gradient u'=du/dx is known to be DU - d(x) over a
given subrange of x, where DU is a constant integer and d(x) is known to lie between 0 and 1.
DU is sometimes refered to as the "integral gear" of the iteration.
we know that Int[u(x+1)] will be either Int[u(x)]+DU-1 or
Int[u(x)]+DU so we can construct iteratively from x and Int[u(x)] .
Given two candidate points (x+1,u+DU-1) and (x+1,u+DU).
we decide which of these is closest to the true curve
{ (x,u) : u=U(x) } by rewriting u=U(x) in the form f(x,u)=0
and considering the sign of a "decision function"
g(x,u) = mf(x+1,u+DU-1+0.5)
where m is a positive integer of convenience often chosen to make g(x,u) integral.
g corresponds to f evaluated half way between the two candidate points.
If g ³ 0 we select (x+1,u+DU) and if g < 0 we select (x+1,u+DU-1).
The idea is of course to iterate g rather than evaluate it anew for each x,u pair.
Example - Division-Free 2D Line Draw
u(x) | =  | (DY/DX)x + c | where 0<DY<DX |
f(x,u) | = | uDX - xDY - cDX | |
g(x,u) | = | 2f(x+1,u+DU-0.5) | where DU=1 |
  | = | 2uDX - 2xDY + DX -2DY -2cDX |
u(x) | =  | (Ax+B)/(Cx+D) |
f(x,u) | = | Cxu + Du - Ax -B |
g(x,u) | = | 2f(x+1,u+DU-0.5) |
  | = | 2(C(x+1)(u+DU-0.5)+D(u+DU-0.5)-Ax-B) |
dg/du | = | 2C(x+1)+2D |
If g(x,u) <0 then we move to x+1,u+DU-1. The next g is then given by
g(x+1,u+DU-1) = g(x,u) + k1(x,u)
where k1(u,x) = 2Cx(DU-1) + 2Cu + C(6DU-5) + 2D(DU-1) - 2A
is maintained by iteration.
[ Proof :
g(x+1,u+DU-1) = 2f(x+2,u+2DU-1.5)
= 2(C(x+2)(u+2DU-1.5)+D(u+2DU-1.5)-A(x+1)-B)
= 2(C(x+1)(u+2DU-1.5)+D(u+2DU-1.5)-Ax-B) + 2(C(u+2DU-1.5)-A)
= 2(C(x+1)(u+DU-0.5)+D(u+DU-0.5)-Ax-B) + 2(C(x+1)(DU-1) +D(DU-1)) + 2(C(u+2DU-1.5)-A)
= g(x,u) + 2(C(x+1)(DU-1) + D(DU-1) + C(u+2DU-1.5)-A)
= g(x,u) + 2Cx(DU-1) + C(6DU-5) + 2D(DU-1) + 2Cu -2A
= g(x,u) + k1(x,u)
.]
If g(x,u) >0 then we instead move to x+1,u+DU. The next g is then given by
g(x+1,u+DU) = g(x,u) + k2(x,u)
where k2(x,u) = 2CxDU + 2Cu + C(6DU-1) + 2DDU - 2A
is maintained by iteration.
[ Proof :
g(x+1,u+DU) = 2f(x+2,u+2DU-0.5)
= 2(C(x+2)(u+2DU-0.5)+D(u+2DU-0.5)-A(x+1)-B)
= 2(C(x+1)(u+2DU-0.5)+D(u+2DU-0.5)-Ax-B) + 2(C(u+2DU-0.5)-A)
= 2(C(x+1)(u+DU-0.5)+D(u+DU-0.5)-Ax-B) + 2(C(x+1)DU +DDU) + 2(C(u+2DU-0.5)-A)
= g(x,u) + 2(C(x+1)DU + DDU + C(u+2DU-0.5)-A)
= g(x,u) + 2(CxDU + 3CDU + DDU + Cu - C/2 -A)
= g(x,u) + 2CxDU + 2Cu C(6DU-1) + 2DDU - 2A
= g(x,u) + k2(x,u)
.]
Note that k1(x,u) = k2(x,u)-2Cx-4C-2D
At a given stage in the iteration we consider g(x,u).
If g(x,u) < 0 we move to x+1,u'=u+DU-1 setting
|
If g(x,u) ³ 0 we move to x+1,u'=u+DU setting
|
It only remains to detect when and how the "gear" DU needs to be changed.
u'(x) = du/dx = (A(Cx+D)-(Ax+B)C)/(Cx+D)2 = (AD-BC)/(Cx+D)2.
u"(x) = d2u/dx2 = -2C(AD-BC)/(Cx+D)3
Thus we see that u' is either always >0, always <0, or always =0 and that
it increases while (Cx+D)<0.
Our two candidate points (x+1,u+DU-1) and (x+1,u+DU) are only the correct ones
if g takes opposite signs when evaluated at them (or is zero at one of them).
Thus the gear needs increasing if both g+k1 and g+k2 are < 0, and decreasing if both
are > 0.
Gear changes will be upwards while Cx+D<0.
Mears asserts
that if 4CDU - 2C > 0 then all subsequent gear changes will be upwards. This is incorrect.
Inverse Functions Via Iterative Approach
In Square Roots the iteration xn+1 = (xn + y/xn)/2 was shown to approach
Öy extremely rapidly.
This is one example of the general iteration
xn+1 = xn + (y-f(xn))/f '(xn)
which tends to x¥=f -1(y).
Infact, if dn=xn-x¥
then dn+1=f"(x)dn2/2f'(x) + O(dn3), and so the convergence is often extremely
rapid.
The following table shows iterations converging to integral roots,
reciprocals, and arctangents.
f(x) | xn+1 | x ¥ | dn+1 | ||||||||||||
xm | ((m-1)xn + y/xnm-1)/m mÖy | (m-1)dn2/2x¥
| p/x | 2xn - yxn2/p | p/y | -ydn2/p
| tan(x) | xn + cos2xn(y-tanxn) | tan-1(y) | -sinx¥dn2
| tan(x) | ((2-sin2xn)xn + (1+cos2xn)y)/2 | tan-1(y) | -sinx¥dn2
| |
f(x) | xn+1 | dn | dn+1 |
1/(1-x) | 1 + xxn | -xn+1 | xdn |
sinx | -x2xn/2n(2n-1) | (-1)nx2n+1/(2n+1)! | -x2dn/2n(2n-1) |
Writing Xn for åi=1n xi we have Xn+1-Xn = Xn-Xn-1 + abXn + ay0, ie. Xn+1(2+ab)Xn+Xn-1 = ay0 which is a second order recurrance relation having specific solution Xn = -y0/b and an auxilliary equation m2-(2+ab)m+1=0 having roots ((2+ab)± Ö (ab(4+ab))/2.
If ab(4+ab)<0 then we have
Xn = Acos(nq)+Bsin(nq)-y0/b where
q= tan-1( Ö (ab(4+ab))/(2+ab)).
We then have
xn = Xn-Xn-1 = A(cos(nq)-cos((n-1)q))+B(sin(nq)sin((n-1)q))
= -2Asin(q/2)sin(nq)+2Bcos(q/2)cos(nq)
= Ccos(nq) + Dsin(nq).
From this general solution for xn we can obtain the specific solution
in terms of x0 and y0.
We can derive yn similarly.
xn = x0cos(nq) + [(x0+ay0-x0cosq)/sinq]sin(nq)
yn = y0cos(nq) + [(bx0+(1-ab-cos+q)y0)/sinq]sin(nq).
If a=-b and a is small then
q= tan-1( Ö (a2(4-a2))/(2-a2)) » tan-1(a)
so that
sinq » a and cosq » 1-a2/2 then these become
xn = x0cos(nq) + y0sin(nq) + x0asin(nq)/2
yn = y0cos(nq) - x0sin(nq) - y0asin(nq)/2.
Typically we might have x0=R, y0=0 giving
xn = Rcos(nq) + Rasin(nq)/2
yn = -Rsin(nq).
This is our desired circle but with a small "error term"
(x0asin(nq)/2,-y0asin(nq)/2).
If the units used are screen pixels, then for this error
term be less than half a pixel and so have no effect we need
a <Min{1/x0,1/y0}.
Thus taking a < 1/R ensures that xn and yn never change by more than one.
To map out the full circle we need 2p/q » 2p/a iterations.
If the reader is wondering why the general result is asymmetric in x
and y it is because the iteration
xn+1 = xn + ayn
yn+1 = yn + bxn+1
is itself asymmetric. The symmetric iteration
xn+1 = xn + ayn
yn+1 = yn + bxn
is far less useful, giving a spiral of increasing radius (if ab<0):
xn = (1-ab)n/2(x0cos(nq) + ((x0(1-cosq)+ay0)/sinq)sin(nq))
yn = (1-ab)n/2(y0cos(nq) + ((y0(1-cosq)+bx0)/sinq)sin(nq))
where q = tan-1(Ö(ab)).
Ellipse Drawing Using Circle Iterations
It is possible to generate a general (ie. non axi-aligned) ellipse
without using any multiplications other than those implicit in two circle
iterations of the type described above.
An axi-aligned ellipse having centre (0,0) has parametric form:
X = acosq
Y = bsinq
and we see immediately that by running two parallel circle iterations,
one from (a,0) and one from (b,0) we can pick the X term from the Xcoord
of the first circle and the Y term from the Ycoord of the second.
This same ellipse rotated by an angle a about the origin has
parametric form:
X = a cosacosq b sinasinq
Y = a sinacosq + b cosasinq
which we can rewrite using A?.? as
X = Acos(q+a) + Bcos(q-a)
Y = Asin(q+a) - Bsin(q-a)
where A=(a+b)/2 ; B=(a-b)/2. So if we iterate one circle of radius A and
another of radius B such that the second circle is 2a behind the first,
then adding the X-coords of the two circles and subtracting the Y-coord
of the second from that of the first will give, respectively, the X and
Ycoords of the desired ellipes. One such pair of circles is obtained by
starting the first from (A,0) and the second from (Bcos(2a),-Bsin(2a)); which
corresponds to starting with q=-a.
If we wish to start from the point on the ellipse having maximal
Ycoordinate then we need to begin with q=tan-1((bcota)/a). This corresponds
to starting the first circle at
((a+b)(ab)sin2a/4C , (a+b)(bcos2a + asin2a)/2C)
and the second at
((a+b)(ab)sin2a/4C , (a-b)(bcos2a-asin2a)/2C)
where C= Ö (b2cos2a + a2sin2a).
Note that q is not the angle subtended with the major axis. That angle is given by
tan-1( b * tanq /a) .
Next : Bezier Curves