There are two types of analytic equiations: explicit
and implicit.
Explicit equations are of the form y = f(x)
whereas implicit equations are of the form f(x,y)
= 0. The former is good for generating points while the
latter is good for testing to see if points satisfy the equation.
The drawback of the explicit form is that it is dependent on the
choice of coordinate axes and is ambiguous if there is more than
one y for a given x (such as y = sqrt(x)).
An alternative to the analytic form of equations is called the parametric form, such as x = f(t), y = g(t). This is good for generating points by varying the parameter t.
Equations can be classified according to the terms contained in them. Equations that only contain variables raised to a power are polynomial equations. If the highest power is one, then the equation is linear. If the highest power is two, then the equation is quadratic. If the highest power is three then it is cubic. If the equation contains sins, cosines, log or a variety of other functions, then it is called transcendental.
Continuity refers to how well behaved the curve is in a mathematical sense. If, for a value arbitrarily close to a value x0 the function is arbitrarily close to f(x0), then it has positional, or zeroth order continuity at that point. If the slope of the curve (or the first derivative of the function) is continuous, then the function has first order continuity. This is extended to all of the functions derivatives.
If a curve is pieced together from individual curve segments then one can speak of piecewise continuity if each individual segment is continuous and the values at the juntions of the curve segments are continuous.
If some data points are used to construct a curve the curve can either pass through the points, in which case the curve is interpolating, or the points can just be used to control the general shape of the curve and the curve doesn't actually pass through the points, in which case the curve is approximating.
Simple linear interpolation is given by
P(t) = (1-u)*P0 + u*P1
Notice that the interpolants, '1-u' and 'u', sum to one. This property insures that the interplating curve (in this case a straight line) falls within the convex hull of the geometric entities being interpolated (in this case the convex hull is the straight line itself). This is called the convex hull property.
Using more general notation, the equation above can be rewritten as:
P(u) = F0(u)*P0 + F1(u)*P1
When written in this form, F0 and F1 are referred to as blending functions. The form is referred to as the geometric formulation because the geometric information, in this case P0 and P1, are explicit in the equation.
The linear interpolation equation can be rewritten as:
P(u) = (P1-P0)*u + P0
Again, using more general notation, the equation is of the form: P(u) = a1*u+ a0 in which the terms are collected into coefficients of the variable to some power. This is the general polynomial form referred to as the algebraic formulation.
Alternatively, both of these forms can be put in a matrix representation . The geometric formulation becomes:
T T
P(t) = [ F0(u) F1(u) ] [P0 P1] = F B
T
(B denotes the transpose of matrix B.)
The algebraic form becomes:
T T
P(t) = [u 1] [a0 a1] = U A
In the fully expanded form, we have:
T T T T
P(u) = [u 1] [ -1 1 ] [P0 P1] = U M B = F B = U A
[ 1 0 ]
Most all of the forms for equations that we will deal with can be written in the format above. Depending on the actual curve type, the U (variable matrix), M (coefficient matrix) and P (geometric information matrix) matrices will have different information in them.
It should be noted that, in general, there is not a linear relationship between changes in the parameter u and the distance traveled along a curve (its arclength). It does happen to be true in the example above concerning a straight line and the parameter u. However, as Mortenson points out, there are other equations which trace out a straight line in space which are fairly convoluted in their relationship between changes in the parameter and distance traveled. The non-linear relationship is almost always true in any parameterized curve unless special care is taken to ensure linearity. (See Reparameterization by Arc Length in the Chapter on Techniques to Aid Motion Specification)
Hermite formulation produces a cubic polynomial through the points. The user needs to supply beginning/end tangents as well as end points.
T
P(u) = U M B
where
3 2
U = [u u u 1]
M = matrix of coefficients found by solving the cubic equation with constraints
P = [P0 P1 P0' P1']
where P0, P0' are the beginning point and tangent vector, respectively, and P1, P1' are the ending point and tangent vector, respectively.
Continuity between beginning and ending tangent vectors of connected segments is ensured by merely using the ending tangent vector of one segment as the beginning tangent vector of the next one. For a single cubic Hermitian curve, we have:
If trying to put a Hermite curve through a large number of points, requiring the user to specify all of the needed tangent vectors can be a burden. There are several techniques to get around this. One such technique is to enforce second degree continuity. This requirement provides enough constraints so that the user doesn't have to provide interior tangent vectors; they can be calculated automatically. See Rogers and Adams' book (p.119-123 in the first edition), or see Mortenson's book (p. 98-112) for an alternative formulation.
The Catmull-Romm curve can be viewed as a Hermite curve in which the tangents at the interior control points are automatically generated according to a relatively simple geometric procedure (as opposed to the more involved numeric techniques referred to above). For each interior point, p[i], the tangent at that point, p'[i], is computed as:
p'[i] = (1/3)*(p[i+1]-p[i-1])
Catmull-Romm spline is a specific type of the cardinal splines. In cardinal splines, the factor of 1/3 is user specified. In addition, the beginning and end tangent vectors can also be automatically determined by various simple geometric means. For example,
p'[0] = (1/3)(p[1]-p[0]) or p'[0] = (1/2)*(2*p[1]-p[2]-p[0])
and these definitions for the tangents can be put in the B matrix of the equation for Hermite interpolation given above.
A cubic curve can also be created in the middle section of four control points in the following way. Take the first three points, P1, P2, P3, and fit a parabola, P(u), through them using the following contraints: P(0) = P1, P(1/2) = P2, P(1) = P3. Take the last three pints, P2, P3, P4, and fit a parabola, R(u), throught them using similar constraints: R(0) = P2, R(1/2) = P3, R(1) = P4. Between points P2 and P3 the two parabolas overlap. Reparametrize this region into the range [0,1] and linearly blend the two parabolic segments:
For interpolating a list of points, interior segments are calculated using the equation above. End conditions can be effected by making parabolic arcs at the very beginning and very end.
This form assumes that all points are equally spaced in parametric space. Often it is the case that even spacing is not present. In such cases, relative cord length can be used to estimate parametric values.
The are several ways to look at Bezier curves. One way to look at cubic Bezier curves is to say that it is similar to the Hermite formulation, except instead of using tangent vectors, the Bezier formulation uses auxiliary control points to define tangent vectors. A cubic curve is defined by four points: P0, P1, P2, P3. P0 and P3 are the beginning and end points, respectively, of the curve. P1 and P2 are the interior control points which define the beginning and ending tangents.
Continuity between adjacent Bezier segments can be controlled by colinearity of the control points on either side of the shared beginning/end point of the two curve segments where they join.
In addition, the Bezier curve formulation allows one to define a curve of arbitrary order. If three interior control points are used, then the resulting curve will be quartic; if four interior control points are used, then the resulting curve will be quintic. See Mortenson for a discussion. For a single cubic Beziercurve:
B-splines are a more general formulation which includes Bezier as a special case. The formulation decouples the number of control points from the degree of the resulting polynomial with the knot vector. Values in the knot vector also allow specification of control points to be interpolated or just approximated by the curve. For a single cubic B-spline curve:
NURBS, which we won't get into here, are even more flexible than basic B-splines. Nonrational Bsplines and Bezier curves are special cases of NURBS. NURBS allow for exact representation of circular arcs whereas Bezier and nonrational Bsplines don't.
Considered by some to be more intuitive control of the curves. P(t) = TMP
Physically-based motion is really just a limited simulation of physical reality. This can be as simple or as complex as the implementor desires. Below are some of the equations of importance in simple physics simulation.
The fundamental equation relating position, distance, and speed is
distance = speed * time
This can be used to control the positioning of an object for any particular frame of the animation because this is tied directly to time
time = frame-number * time-per-frame
The average velocity of a body is the distance moved divided by the time it took to move, and is given by
v = (distance between s(t1) and s(t2))/(t2-t1) (B.1)
where s(t) is a function that gives the position of the object at time t. Notice that the units of velocity is distance per time, e.g., ft/sec.
The instantaneous velocity is determined by taking equation (B.1) and moving t2 closer and closer to t1. In the limiting case this becomes the derivative of the function s(t) with respect to time.
Similarly, the average acceleration of an object is the change in velocity divided by the time it took to effect the change, and is given by
a = (v(t2)-v(t1))/(t2-t1)
where v(t) is a function that gives the velocity of the object at time t. Notice that the units of acceleration is velocity per time or distance per time per time, e.g., ft/sec2.
In the same way, instantaneous acceleration is the derivative of v(t) with respect to t so that we now have
a(t) = v'(t) = s''(t)
For example, in the case of motion due to gravity, we have:
s(t) = (1/2)gt**2 v(t) = gt a(t) = g
where g is the acceleration due to gravity, a constant that has been measured to be 32 ft/sec2.
Circular motion is an important type of motion in physics and comes up when discussing cases from planets to robotics. A simple way to describe it is to specify position based on polar coordinates. The position of a particle orbiting the origin a distance r from it can be described as:
p = (r*cos(theta))i + (r*sin(theta))j (B.2)
where i and j are orthonormal unit vectors (at right angles to each other and unit length). p is the positional vector of the particle. In a constant radius curcular orbit, the distance r is constant, but it is important to keep in mind that p, being a vector, is not.
Uniform circular motion means that theta changes at a constant rate which we call omega:
d(theta)/dt = omega
or theta = omega*time
If theta increases with time, the quantity omega is positive as is called the angular speed of the object. If theta is measured in radians and time in seconds, then omega is measured in radians per second.
The angular speed, omega, is related to the time T that it takes to complete one revolution. T is called the period and is related by the following equation:
2PI = omega*T
since omega moves through 2PI radians in one revolution.
Given (B.2), we can take the derivative of it with respect to time to get the acceleration. After substituting (B.3) and calling it p(t), using w to represent omega and taking the derivative, we get:
v(t) = dp/dt = (-r*w*sin(w*t))i + (r*w*cos(w*t))j (B.4)
Notice that the velocity vector, v(t), is perpendicular to the position vector p(t) (this can be demonstrated by taking the dot product of the two vectors v(t) and p(t) to get zero. By computing the length of v(t) we arrive at:
v = r*w
showing that the velocity is independent of t and is, therefore, constant.
Notice, however, that a constant circular motion still gives rise to an acceleration. We take the derivative of (B.4) to get:
a(t) = dv/dt = -w2*((r*cos(w*t))i + (r*sin(w*t))j) = w**2*r
using v = r*w, we have
a(t) = v**2/r
which is called the centripetal acceleration which is directed radially inward and has constant magnitude.
For any particle in a rigid mass undergoing a rotation, that particle is undergoing the same rotation about its own center and, in addition, because it may be displaced from the center of rotation, it is undergoing an instantaneous positional translation.
First Law: If no force is acting on an object, then it will not accelerate. It will maintain a constant velocity.
This is the principle of inertia.
Second Law: The change of motion of an object is proportional to the forces applied to it.
F = m*a
The change in motion (acceleration) of an object due to a force is proportional to the object's mass. That is, the change in motion involved not only the body's velocity but also its mass. Stated another way this law becomes
F = d(m*v)/dt
where the quantity mv is the object's momentum. In this form the law states 'Force is equal to the rate of change of momentum.
It is important to note that the force F used here is considered to be the sum of all external forces acting on an object and that these equations really represent three sets of equations, one for each coordinate:
Fx = m*ax Fy = m*ay Fz = m*az
Third Law: To every action there is always opposed an equal and opposite reaction.
F12 = -F21 where F12 is the force Body 1 exerts on Body 2 and vice versa
Push on anything, and it pushes back with an equal force
An inertial frame is a frame of reference (e.g., coordinate system) in which the principle of inertia holds. Any frame that is not accelerating is an inertial frame. An observer in an inertail frame deduces the same laws of motion and has no way of determining whether he is at rest or moving in an absolute sense (but then, what is moving in an absolute sense?)
The tendency of a force to produce rotation is called torque. In 2D, we have:
t = F*d where t = torque F = force perpendicular to arm d = length of arm
In the 2D case, the convention is to call counterclockwise torques positive and clockwise torques negative. Considering the full 3D situation:
t = r x F where r = is vector to point of application of force F= force vector t = torque vector
It is important to note that torque refers to a specific axis of rotation. The same force can exert different torques about different axes of rotation. The torque vector for a rigid body is position independent. That is, for any particle in a rigid mass, the particle is undergoing that same torque
In the absense of acceleration, there is equilibrium. Study of bodies in equilibrium is statics. In such as case the net force and torque acting on a body vanishes:
SUM of F = 0 SUM of t = 0
There may be forces present, but the vector sum is equal to 0.
Law of universal gravitation is
F = G*m1*m2/r**2 where G = 6.67x10-11 m3/kg sec2
The resultant force on a mass is the vector sum of the individual forces.
When two spherically symmetric objects are not touching, each acts as if all its mass were concentrated at its center.
The force of gravity from a mass distributed symmetrically over a spherical shell is zero anywhere inside the shell.
Centripetal force is any force which is directed inward, toward the center of an object's motion (cetripetal means 'center seeking'). In the case of a body in orbit, gravity is the cetripetal force that holds the body in the orbit. Consider a body, such as the moon, with mass Mm and a distance r away from the earth. The unit vector from the earth to the moon is R. The earth has a mass Me. The acceleration induced by a circular orbit always points toward the center (centripetal) and has magnitude v2/r:
a = (-v**2/r)R
The force inducing the acceleration is gravity and the force of gravity is:
F = (-GMmMe/(r*r))R
and it causes the moon to have acceleration
a = F/Mm
or a = (-GMe/(r*r))R
using the original equation for acceleration of circular motion and canceling terms we can solve for the velocity of the moon which compensates for the effect of gravity and establishes an orbit to be:
v = sqrt(GMe/r) Similar analyses can be carried out when considering a ferris wheel, an airplane banking to make a turn or the motion of a conical pendulum.
While gravity is one of the four fundamental forces (gravity, electromagnetic, strong, and weak) which act at a distance. A second important category is contact forces.The tension in a wire or rope is an example of a contact force, as is the compression in a rigid rod. These forces arise from the complex interactions of electric forces which tend to keep atoms a certain distance apart. The empirical law which governs such forces, and which is familiar when dealing with springs, is Hook's Law:
F = -kx where x = change in lenght from the equilibrium value of the spring k = spring constant F = restoring force of sping to return to rest position
The constant k is a measure of the stiffness of the spring; the larger k is, the more sensitive the spring is to motion away from the rest position.
Another example of contact force is the result of repulsion between any two objects pressed against each other, known as the normal force. It arises from the repulsion of the atoms of the two objects, is always perpendicular to the surfaces and its magnitude is proportional two how hard pressed the two objects are. Other important examples of contact forces are friction and viscosity.
Frictional forces from a surface do not exceed an amount proportional to the normal force exerted by the surface on the object:
fs <= �sn where �s="coefficient" of static friction n="normal" force�s varies with the two materials that are in contact. The frictional forces acting between surfaces at rest with respect to each other are called forces of static friction. For a force applied to an object sitting on another object that is parallel to the surfaces in contact, there will be a specific force at which the block starts to slip. As F increases from zero up until that threshold force, the lateral force is conteracted by an equal force of friction in the opposite direction. Once the object begins to move Kinetic friction acts on the object and approximately obeys the empirical law:
fk = �kN where �k = coefficient of kinetic friction N = normal force
Kinetic friction is typically less than static friction. The force of kinetic friction is always opposite to the velocity of the object.
The resistive force is viscosity. It is another contact force and extremely difficult to model accurately. When an object moves at low velocity through a fluid, the viscosity is approximately proportional to the velocity:
Fvis = -Knv where k = proportionality constant (depends on size and shape of object) n = coefficient of viscosity (depends on fluid)
The coefficient of viscosity, n, decreases with increasing temperature for liquids, and increases with temperature for gases.
Stokes' law for a sphere of radius R gives:
K = 6�R
An object dropping through a liquid attains a constant speed, called the limiting or terminal velocity, at which gravity, acting downward, and the viscous force, acting upward, balance each other (no acceleration):
mg = 6�Rnv v = mv/6�Rn
In a viscous medium, heavier bodies fall faster than lighter bodies. For sphereical objects falling a low velocity in a viscous medium, not necessarily at terminal velocity:
mdv/dt = mg - 6�Rnv
Consider a frame of reference that rotates in uniform circular motion with respect to an inertial frame, a post, because it is held at a constant distance by a rope. Relative to the inertial frame, each point in the uniformly rotating frame has centripetal acceleration:
a = -(v**2/r)R where r is the distance the point is from the axis of rotation R is the unit vector from the inertial frame to the rotating frame v is the speed of the point
And the tension in the rope supplies the force necessary to produce the centripetal acceleration. Relative to the rotating frame (not an inertial frame), the frame itself does not move and therefore to counteract the force supplied by the rope is the centrifugal force:
Fc = -ma = mw**2rR
For a constant force of magnitude F moving an object a distance h parallel to the force, the work W done by the force:
W = Fh
If a mass m is lifted up so that it doesn't accelerate, then the lifting force is equal to the weight (mg) of the object. Since mg is constant, the work done to raise the object up to a height h is
W = mgh
Energy that a body has by virtue of its location is called potential energy. The work in this case is converted into potential energy.
Energy of motion is called kinetic energy. For a falling body that started at a height h,
v**2 = 2gh (1/2)mv**2 = mgh K = (1/2)mv**2
Potential plus kinetic energy of a closed system is conserved: Ua + Ka = Ub + Kb
Power is defined to be the rate at which work is done with respect to time. P = dW/dt
In a closed system, momentum (mass times velocity) is conserved: d(m1v1 + m2v2 + ... + mnvn)/dt = 0 m1v1 + m2v2 + ... + mnvn = const. The center of mass r of a system of particles moves under external forces as if all the mass were concentrated at r and all forces acted at r.
[Need collision response material]
Stability of an object subject to a linear restoring force:
F = -kx
Associated with this force is the potential energy:
U = (1/2)kx2
These systems have the property that if they are distrubed from equilibrium, the restoring force that acts on them tends to move them back into equilibrium. When disturbed from equilibrium, they tend to overshoot that point when they return, on account of inertia. Then the restoring force acts in the opposite direction, trying to return the system to equilibrium. The result is that the system winds up oscillating back and forth like a mass on the end of a spring.
We have the following:
F = ma F = -kx a = d2x/dt2 m d2x/dt2 = -kx d2x/dt2 = -kx/m
is a differential equation satisfied by the displacement x(t).
We can model the damping force after Stokes' law, in which resistance is assumed to be proportional to velocity. This is usually valid for oscillations of sufficiently small amplitude. The damping force opposes the motion and can be written as:
Fd = -�*dx/dy
where � is aconstand called the damping coefficient. Adding the damping force to the spring force, Newton's second law gives us:
m*d2x/dt2 = -k*x - �*dx/dt
After dividing through by m and rearranging:
d2x/dt2 + �*dx/dt + w20*x = 0
where � = �/m and w20 = k/m.
If there is a sping force but no damping, the general solution can be written as C*cos(w0*t+�0). If there is damping but no spring force, the general solution turns out to be x = C*e-�t+D with C and D constant. If both the spring force and damping force are present, the solution is more compicated.
L = r x p where L = angular momentum r = vector from reference point p = mv, momentum
The time rate of change of angular momentum equals the torque:
= dL/dt
Angular momentum, just like linear momentum, is conserved in a closed system:
�pi = const �Li = const
The center of mass r of a system of particles moves under external forces as if all the mass were concentrated at r and all forces acted at r.
[Need some extensive material on this - it arises in robotics as well as rigid body simulation]
This is a collection of techniques and pointers to other sources that should help the reader understand and implement numerical techniques that are useful for computer animation algorithms. This is intended just to give the reader a starting point.
There is a hierarchy of complexity of numerical problems. In increasing order of complexity:
Linear Algebraic Equations Non-linear Algebraic Equations Ordinary Differential Equations Partial Differential Equations
We will treat each of these in turn to give an overview of various approaches to handle problems in these domains. Usually the techniques will, by approximations, reduce a problem from one domain to a simpler domain until there is a procedure for determining a solution in the simpler domain.
Linear algebraic equations can be expressed in matrix form, for example, the familiar
Ax = b
where A and b are a constant matrix and vector respectively and x is a variable vector which is solved for.
Non-linear algebraic equations involve no derivatives but are more complex than simple linear equations. For example, polynomials fall into this class. Some non-linear equations can be solved directly, for example, a quadratic polynomial can be solved using the quadratic equation. Others can be converted into linear algebraic equations by using simplifying assumptions and then solved using the linear techniques.
Ordinary differential equations (ODEs) involve the derivative of only one variable. Some ODEs can be solved directly. ODEs can be converted into non-linear agebraic equations using simplifying assumptions and solved from there.
Partial differential equations (PDEs) involve functions of more than one variable and derivatives of more than one variable. There are techniques which convert PDEs to ODEs.
A common useful equation that comes up is Taylor's series expansion:
f(x+�) = f(x) + f'(x)� + (1/2)f''(x)�2 + ...
This equation comes in handy when the current situation is known or can be approximated, and future conditions must be calculated. An example is the case where the motion of an object is controlled by physical laws and f(x) is the position of that object, f'(x) is the velocity and f''(x) is the acceleration (e.g., under gravity). The position of the object at the next instant of time (i.e., x+�) is to be computed, acceleration is given, for example, by gravity and velocity can be approximated by adding the acceleration to the last velocity.
Systems of linear equations arise, for example, when considering a mesh of mass points connected by springs for flexible body animation. Each pair of mass points considered to be connected by a spring gives rise to a linear equation of the form F = -kx.
Linear systems of equations can be written
Ax = b where A is a two-dimensional matrix of coefficients x is a vector of unknowns b is a vector of constants
The most obvious approach to this problem may appear to be to calculate the inverse of A and multiply both sides by it:
x = A-1b
The inverse of A, A-1, can be computed by Gauss-Jordan elimination with either partial or full pivoting. It is more efficient, however, to do Gaussian elimination until an upper triangular matrix results and follow with backsubstitution. Both of these solutions require that the right hand side of the original equations be known in advance. A better approach is to decompose A into an upper triangular matrix and a lower triangular matrix:
A = LU where L is lower triangular U is upper triangular
so that
Ax = (LU)x = L(Ux) = b
Note that the LU decomposition of A can be done without knowing the b vector. The reason for doing this is that solving a system of equations with a triangular matrix is accomplished by simple substitiution. So we first solve
Ly = b
called forward substitution, and then
Ux = y
called backsubstitution.
The only thing left to do is specify how the LU decomposition is to be done. In order to keep the calculation stable, the decomposition should be done with pivoting so that a permutation of A is actually computed.
Singular value decomposition is a powerful set of techniques for solving most linear least squares problems and for dealing with systems of equations are either singular or else numerically very close to being singular. The singular value decomposition of a matrix which is MxN, where the number of rows M is greater than or equal to the number of columns N, is an MxN column-orthogonal matrix U, an NxN diagonal matrix W with positive or zero elements, and the transpose of an NxN orthogonal matrix V. Zeroing out small values in W essentially well-conditions the resulting set of linear equations and very often produces a better solution that both the direct-method solution and the SVD solution where the small w's are left nonzero.
Non-linear functions arise in a variety of places, most notably when dealing with higher order surfaces.
In a few cases, as mentioned earlier in this Appendix, non-linear algebraic equations have analytic solutions, such as the quadratic equation. However, most of the time analytic solutions are not known for equations of this type. Simlar to the linear case, a simple way to solve non-linear equation is by functional iteration. If you have an equation in the form of
y = f(y)
you can use an initial guess and use it to generate the next estimate to the answer:
yk+1 = g(yk)
However, functional iteration has poor convergence properties.
A better technique for solving non-linear functions is Newton's method which tries to find zero of non-linear function by truncating the Taylor's series after two terms. !!!
There are other non-linear equation solvers, such as the variable metric method. !!!
Such hill-climbing techniques are applicable in function of two variables because, from a given point on a curve there are only two directions you can proceed. With higher dimensional functions, the problem is much more complicated because a surface is defined and there are an infinite number of directions to proceed from any one point.
Problems involving ordinary differential equations (ODEs) can always be reduced to the study of sets of first order differental equations. For example, the secnd order equation
d2y/dx2 + q(x)dy/dx = r(x)
can be rewritten as two first order equations
dy/dx = z(x) dz/dx = r(x) - q(x)z(x)
where z is a new variable. The usual choice for new variables is to let them be derivatives of eadch other and of the original variable. Occassionally, it is useful to incorporate some other factors into their definition fro the purpose of mitigating singular behavior that could result in overflows or increased roundoff errors.
The generic problem in ordinary differential equations is thus reduced to the study of a set of N coupled first order differential equations.
dyi(x)/dx = fi(x,y1,...,yN) i = 1,...,N
The nature of the problem's boundary conditions are crucial in determining how to attack the problem. Boundary conditions are algebraic conditions on the values of the functions yi. In general they can be satisfied at discrete specified points, but do not hold between those points, i.e., are not preserved automatically by the differential equations. Boundary conditions divide into two broad categories:
The underlying idea of any routine for solving the initial value problem is always to rewrite the dy's and dx's as finite steps �y and �x and multiply the equations by �x. In the limit of making the stepsize very small, a good approximation to the underlying differential equation is achieved. Literal implementation of this procedure results in Euler's method whcih is not recommended for any practical use.
Three major types of practical numerical methods for solving initial value prolems for ODEs are:
Partial differential equations are at the heart of many computer analyses or simulations of continuous physical systems, such as fluids, electromagnetic fields, the human body, tec. Partial differential equations are classified into three categories, hyperbolic, parabolic and elliptic, on the basis of their characteristics, or curves of information propagation. The prototypcial example of a hyperbolic equations is the one-dimensional wave equation:
�2u/�t2 = v2(�2u/�x2)
where v = constant is the velocity of wave propagation. The prototypical parabolic equation is the diffusion equation:
�u/�t = -�(D(�u/�x)/�x
where D is the diffusion coefficient. The prototypical elliptic equation is the Poisson equation:
�2u/�x2 + �2u/�y2 = p(x,y)
where the source term p is given. If the source term is equal to zero, the equation is Laplace's equation.
The first two are initial value problems which describe how u(x,t) propagates itself forward in time. The last one is a boundary value problem where the task is to find a single static function u(x,y) which satisfies the equation within some (x,y) region of interest and which has some desired behavior on the boundary of the region of interest.
Methods for solving PDEs are finite differencing, finite element, Monte Carlo, spectral and variational methods.
A = Axi + Ayj + Azk where i, j, k are unit vectors along the principle axes, x, y, and z, respectively.
|A| = �(Ax2 + Ay2 + Az2)
A.B = |A||B|cos� where � is the angle from A to B, 0 <= � <="180" a.b="Ax*Bx" + ay*by + az*bz a.a="|A|2" a.b="0" if and only if a and b are perpendicular (i.e., cos(�)="cos(90)" or a="0" or b="0" a.i="Ax" a.j="Ay" a.k="Az"Cross Product | i j k | A x B = | Ax Ay Az | | Bx By Bz | A x B = (AyBz - AzBy)i + (AzBx - AxBz)j + (AxBy - AyBx)k A x B = -(B x A) the direction of A x B is perpendicular to both and is determined by the right hand rule (if A and B are in right-hand space) |A x B| = |A||B|sin� where � is the angle from A to B, 0 <= � <="180" |a x b|="0" if and only if a and b are parallel (i.e., sin(�)="sin(0)" or a="0" or b="0"
linear objective function unconstrained - no solution constrained linear constraints - linear programming: simplex method; integer programming non-linear constraints - quadratic programming non-linear object function unconstrained unconstrained search techniques 1D: binary search; Fibonocci search; Golden search first derivative methods nD: gradient search; steepest descent second derivative methods Jacobian & Hessian constrained linear constraints non-linear constraints
Go back to Table of Contents
Go to Last Chapter: Chapter 5. Algorithms and Models to Fully Specify Motion
Go to References