If we split the time interval (0,T) into N segments of length l =T/N
and iterate
x_{new} = x_{old} + lv_{old} ; v_{new} = v_{old} + la.
or more formally
x_{n+1} = x_{n} + lv_{n} ; v_{n+1} = v_{n} + 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
x_{n+1 }= x_{n} + x_{n1} together with "initial condition" x_{1} = x_{0} = 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 x_{n1} (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
x_{n} = x_{0} + nlu + n(n1)l^{2}a/2 ;
v_{n} = v_{0} + nla
since v_{n} is readily solved and then x_{n} follows. We see from this that x_{n} is not the correct position of the body
at time t=nl.
The correct value is x_{0} + nlu + n^{2}l^{2}a/2 but this differs by
only nl^{2}a/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 v_{0} = u + la/2.
(ii) Use the following slightly modified iteration:
v_{temp} = v_{old} + la/2
x_{new} = x_{old }+ lv_{temp}
v_{new} = v_{temp} + 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(x_{n}_{+k},...,x_{n+1},x_{n}) = f(n)
for some linear function g and arbitary function f we find the general
solution to g(x_{n}_{+k},...,x_{n+1},x_{n}) = 0 and then add to it any specific
solution to g(x_{n}_{+k},...,x_{n+1},x_{n}) = f(n), which is often a function of
similar form to f(n).
The general solution to ax_{n+1} + bx_{n} = 0 (a¹0) is x_{n}= A (b/a)^{n}.
The form of the general solution to
ax_{n}_{+2 }+ bx_{n+1} +cx_{n} = 0 (a¹0)
depends on the nature of the roots of the auxiliary equation am^{2}+bm+c=0.
If the roots are repeated (b^{2}=4ac) then the general solution is
x_{n} = (A+Bn)m^{n} where m is the repeated root b/2a.
If the roots are real and distinct (b^{2}>4ac) then the general solution
is x_{n}=Am_{1}^{n} + Bm_{2}^{n} where m_{1} and m_{2} are the two roots (b ± Ö (b^{2}4ac))/2a.
If the roots are imaginary (b^{2}<4ac) then the general solution is
x_{n} = (c/a)^{n/2}(Acos(nq)+Bsin(nq)) or equivalently x_{n} = (c/a)^{n/2}Acos(nq+e) where
tanq =  Ö(4acb^{2})/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 = r_{0} + (1e^{st})(v_{0}a/s)/s + ta/s
v = (v_{0}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 + (r_{0} + (v_{0}a/s))
v = a/s
We would typically simulate this with the iteration
r_{n+1} = r_{n} + lv_{n}
v_{n+1} = dv_{n} + 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
r_{n+1} = r_{n} + l_{1}v_{n}
v_{n+1} = dv_{n} + l_{2}a.
(by taking l_{1}=l, l_{2}=dl, this covers the alternative iteration
r_{n+1} = r_{n} + lv_{n} ;
v_{n+1} = d(v_{n} + la) ).
To find the explicit forms for this iteration we first find v_{n} as the
solution of v_{n+1}  dv_{n} = l_{2}a.
The general solution is thus Ad^{n} 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 xdx = l^{2}a and
obtain x=l^{2}a/(1d).
Our full solution is therefore
v_{n} = Ad^{n} + l^{2}a/(1d) and
it remains only to satisfy the initial conditions, ie. to find a value for A that gives
the desired value for v_{0} .
Once v_{n} is solved we use it to find r_{n} using results A1.1 and A1.9.
The results are
l_{1}(1d^{n})  l_{1}l_{2}  æ  (1d^{n})  ö  
r_{n}  =  r_{0}  +  ¾¾¾¾  v_{0}  +  ¾¾  ç  n  ¾¾  ÷  a 
(1d)  (1d)  è  (1d)  ø 
(1d^{n})  
v_{n}  =  d^{n}v_{0}  +  l_{2}  ¾¾  a 
(1d) 
l_{1}l_{2}  æ  l_{2}  l_{1}l_{2}  ö  
r_{n}  =  n  ¾¾  a  +  ç  r_{0}  +  ¾¾  v_{0}    ¾¾  a  ÷ 
(1d)  è  (1d)  (1d)  ø 
l_{2}  
v_{n}  =  ¾¾  a 
(1d) 
n  x_{3n}  x_{2n}  x_{1n}  x_{0n} 
0  d  al^{2}+ bl+c  6al+2b  6a 
1  al^{3}+ bl^{2}+ cl+d  7al^{2}+3bl +c  12al +2b  
2  8al^{3}+4bl^{2}+2cl+d  19al^{2}+5bl+c  
3  27al^{3}+9bl^{2}+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)=at^{3}+bt^{2}+ct+d
for t=n/N " n Î Z_{N}.
To do this we iterate
x_{0n}' = x_{00} " n ³ 0
x_{1n}' = x_{1(n1)}' + x_{00}' " n ³ 1
x_{2n}' = x_{2(n1)}' + x_{1(n1)}' " n ³ 1
x_{3n}' = x_{3(n1)}' + x_{2(n1)}' " n ³ 1
which will work if
x_{3n}' = x_{3n} ;
x_{2n}' = lx_{2n} ;
x_{1n}' = l^{2}x_{1n} ;
x_{0n}' = l^{3}x_{0n} " n ³ 1
so we start from the initial values
x_{00}' = l^{3}x_{00},
x_{10}' = l^{2}x_{10},
x_{20}' = lx_{20},
x_{30}' = x_{30} 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 x_{3n}'.
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 D_{U}  d(x) over a
given subrange of x, where D_{U} is a constant integer and d(x) is known to lie between 0 and 1.
D_{U} is sometimes refered to as the "integral gear" of the iteration.
we know that Int[u(x+1)] will be either Int[u(x)]+D_{U}1 or
Int[u(x)]+D_{U} so we can construct iteratively from x and Int[u(x)] .
Given two candidate points (x+1,u+D_{U}1) and (x+1,u+D_{U}).
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+D_{U}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+D_{U}) and if g < 0 we select (x+1,u+D_{U}1).
The idea is of course to iterate g rather than evaluate it anew for each x,u pair.
Example  DivisionFree 2D Line Draw
u(x)  =  (D_{Y}/D_{X})x + c  where 0<D_{Y}<D_{X} 
f(x,u)  =  uD_{X}  xD_{Y}  cD_{X}  
g(x,u)  =  2f(x+1,u+D_{U}0.5)  where D_{U}=1 
=  2uD_{X}  2xD_{Y} + D_{X} 2D_{Y} 2cD_{X} 
u(x)  =  (Ax+B)/(Cx+D) 
f(x,u)  =  Cxu + Du  Ax B 
g(x,u)  =  2f(x+1,u+D_{U}0.5) 
=  2(C(x+1)(u+D_{U}0.5)+D(u+D_{U}0.5)AxB)  
dg/du  =  2C(x+1)+2D 
If g(x,u) <0 then we move to x+1,u+D_{U}1. The next g is then given by
g(x+1,u+D_{U}1) = g(x,u) + k_{1}(x,u)
where k_{1}(u,x) = 2Cx(D_{U}1) + 2Cu + C(6D_{U}5) + 2D(D_{U}1)  2A
is maintained by iteration.
[ Proof :
g(x+1,u+D_{U}1) = 2f(x+2,u+2D_{U}1.5)
= 2(C(x+2)(u+2D_{U}1.5)+D(u+2D_{U}1.5)A(x+1)B)
= 2(C(x+1)(u+2D_{U}1.5)+D(u+2D_{U}1.5)AxB) + 2(C(u+2D_{U}1.5)A)
= 2(C(x+1)(u+D_{U}0.5)+D(u+D_{U}0.5)AxB) + 2(C(x+1)(D_{U}1) +D(D_{U}1)) + 2(C(u+2D_{U}1.5)A)
= g(x,u) + 2(C(x+1)(D_{U}1) + D(D_{U}1) + C(u+2D_{U}1.5)A)
= g(x,u) + 2Cx(D_{U}1) + C(6D_{U}5) + 2D(D_{U}1) + 2Cu 2A
= g(x,u) + k_{1}(x,u)
.]
If g(x,u) >0 then we instead move to x+1,u+D_{U}. The next g is then given by
g(x+1,u+D_{U}) = g(x,u) + k_{2}(x,u)
where k_{2}(x,u) = 2CxD_{U} + 2Cu + C(6D_{U}1) + 2DD_{U}  2A
is maintained by iteration.
[ Proof :
g(x+1,u+D_{U}) = 2f(x+2,u+2D_{U}0.5)
= 2(C(x+2)(u+2D_{U}0.5)+D(u+2D_{U}0.5)A(x+1)B)
= 2(C(x+1)(u+2D_{U}0.5)+D(u+2D_{U}0.5)AxB) + 2(C(u+2D_{U}0.5)A)
= 2(C(x+1)(u+D_{U}0.5)+D(u+D_{U}0.5)AxB) + 2(C(x+1)D_{U} +DD_{U}) + 2(C(u+2D_{U}0.5)A)
= g(x,u) + 2(C(x+1)D_{U} + DD_{U} + C(u+2D_{U}0.5)A)
= g(x,u) + 2(CxD_{U} + 3CD_{U} + DD_{U} + Cu  C/2 A)
= g(x,u) + 2CxD_{U} + 2Cu C(6D_{U}1) + 2DD_{U}  2A
= g(x,u) + k_{2}(x,u)
.]
Note that k_{1}(x,u) = k_{2}(x,u)2Cx4C2D
At a given stage in the iteration we consider g(x,u).
If g(x,u) < 0 we move to x+1,u'=u+D_{U}1 setting

If g(x,u) ³ 0 we move to x+1,u'=u+D_{U} setting

It only remains to detect when and how the "gear" D_{U} needs to be changed.
u'(x) = du/dx = (A(Cx+D)(Ax+B)C)/(Cx+D)^{2} = (ADBC)/(Cx+D)^{2}.
u"(x) = d^{2}u/dx^{2} = 2C(ADBC)/(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+D_{U}1) and (x+1,u+D_{U}) 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+k_{1} and g+k_{2} are < 0, and decreasing if both
are > 0.
Gear changes will be upwards while Cx+D<0.
Mears asserts
that if 4CD_{U}  2C > 0 then all subsequent gear changes will be upwards. This is incorrect.
Inverse Functions Via Iterative Approach
In Square Roots the iteration x_{n}_{+1} = (x_{n} + y/x_{n})/2 was shown to approach
Öy extremely rapidly.
This is one example of the general iteration
x_{n}_{+1} = x_{n} + (yf(x_{n}))/f '(x_{n})
which tends to x_{¥}=f^{ 1}(y).
Infact, if d_{n}=x_{n}x_{¥}
then d_{n+1}=f"(x)d_{n}^{2}/2f'(x) + _{O}(d_{n}^{3}), and so the convergence is often extremely
rapid.
The following table shows iterations converging to integral roots,
reciprocals, and arctangents.
f(x)  x_{n}_{+1}  x_{ ¥ }  d_{n+1} 
x^{m}  ((m1)x_{n }+ y/x_{n}^{m1})/m  ^{m}Öy  (m1)d_{n}^{2}/2x_{¥} 
p/x  2x_{n}  yx_{n}^{2}/p^{}  p/y  yd_{n}^{2}/p 
tan(x)  x_{n} + cos^{2}x_{n}(ytanx_{n})  tan^{1}(y)  sinx_{¥}d_{n}^{2} 
tan(x)  ((2sin2x_{n})x_{n} + (1+cos2x_{n})y)/2  tan^{1}(y)  sinx_{¥}d_{n}^{2} 
f(x)  x_{n+1}  d_{n}  d_{n+1} 
1/(1x)  1 + xx_{n}  x^{n+1}  xd_{n}^{ } 
sinx  x^{2}x_{n}/2n(2n1)  (1)^{n}x^{2n+1}/(2n+1)!  x^{2}d_{n}/2n(2n1) 
Writing X_{n} for å_{i=1}^{n} x_{i} we have X_{n+1}X_{n} = X_{n}X_{n1} + abX_{n} + ay_{0}, ie. X_{n+1}(2+ab)X_{n}+X_{n1} = ay_{0} which is a second order recurrance relation having specific solution X_{n} = y_{0}/b and an auxilliary equation m^{2}(2+ab)m+1=0 having roots ((2+ab)± Ö (ab(4+ab))/2.
If ab(4+ab)<0 then we have
X_{n} = Acos(nq)+Bsin(nq)y_{0}/b where
q= tan^{1}( Ö (ab(4+ab))/(2+ab)).
We then have
x_{n} = X_{n}X_{n1} = A(cos(nq)cos((n1)q))+B(sin(nq)sin((n1)q))
= 2Asin(q/2)sin(nq)+2Bcos(q/2)cos(nq)
= Ccos(nq) + Dsin(nq).
From this general solution for x_{n} we can obtain the specific solution
in terms of x_{0} and y_{0}.
We can derive y_{n} similarly.
x_{n} = x_{0}cos(nq) + [(x_{0}+ay_{0}x_{0}cosq)/sinq]sin(nq)
y_{n} = y_{0}cos(nq) + [(bx_{0}+(1abcos+q)y_{0})/sinq]sin(nq).
If a=b and a is small then
q= tan^{1}( Ö (a^{2}(4a^{2}))/(2a^{2})) » tan^{1}(a)
so that
sinq » a and cosq » 1a^{2}/2 then these become
x_{n} = x_{0}cos(nq) + y_{0}sin(nq) + x_{0}asin(nq)/2
y_{n} = y_{0}cos(nq)  x_{0}sin(nq)  y_{0}asin(nq)/2.
Typically we might have x_{0}=R, y_{0}=0 giving
x_{n} = Rcos(nq) + Rasin(nq)/2
y_{n} = Rsin(nq).
This is our desired circle but with a small "error term"
(x_{0}asin(nq)/2,y_{0}asin(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/x_{0},1/y_{0}}.
Thus taking a < 1/R ensures that x_{n} and y_{n} 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
x_{n+1} = x_{n} + ay_{n}
y_{n+1} = y_{n} + bx_{n+1}
is itself asymmetric. The symmetric iteration
x_{n+1} = x_{n} + ay_{n}
y_{n+1} = y_{n} + bx_{n}
is far less useful, giving a spiral of increasing radius (if ab<0):
x_{n} = (1ab)^{n/2}(x_{0}cos(nq) + ((x_{0}(1cosq)+ay_{0})/sinq)sin(nq))
y_{n} = (1ab)^{n/2}(y_{0}cos(nq) + ((y_{0}(1cosq)+bx_{0})/sinq)sin(nq))
where q = tan^{1}(Ö(ab)).
Ellipse Drawing Using Circle Iterations
It is possible to generate a general (ie. non axialigned) ellipse
without using any multiplications other than those implicit in two circle
iterations of the type described above.
An axialigned 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(qa)
Y = Asin(q+a)  Bsin(qa)
where A=(a+b)/2 ; B=(ab)/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 Xcoords of the two circles and subtracting the Ycoord
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)(bcos^{2}a + asin^{2}a)/2C)
and the second at
((a+b)(ab)sin2a/4C , (ab)(bcos^{2}aasin^{2}a)/2C)
where C= Ö (b^{2}cos^{2}a + a^{2}sin^{2}a).
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