Numerical Methods For Particle Tracing in Vector Fields: On-Line Visualization Notes
Numerical Methods For Particle Tracing in Vector Fields: On-Line Visualization Notes
Numerical Methods For Particle Tracing in Vector Fields: On-Line Visualization Notes
Kenneth I. Joy
Visualization and Graphics Research Laboratory
Department of Computer Science
University of California, Davis
Overview
If we consider a vector field in space and we have a rectilinear volume whose vertices are sampled from
this vector field, one method of visualizing the field is to place particles in the field and animate their motion
through the field. This time animation is the basis for a number of visualization methods – streamlines,
stream surfaces, and stream ribbons. In this paper we review the numerical techniques necessary for tracing
the path of a particle through the field. We will study these techniques in two-dimensional vector fields, but
they are directly extensible to three-dimensional fields.
!v0,1
!v1,1
p0,1 p1,1
!v
p
!v0,0 v
p0,0 p1,0
u
1 !v1,0
We can calculate the vector at p by using the values u and v and bilinearly interpolating the vectors at
the corners. Here
!v = (1 − v)!v0∗ + v!v1∗
where
or
!v = (1 − v)(1 − u)!v0,0 + (1 − v)u!v1,0 + v(1 − u)!v0,1 + vu!v1,1
Euler’s Method
Given a point p in a cell of our rectilinear grid, an obvious way to trace a particle through a cell is to
step along the vector !v associated with p using a very short step size ∆t. The new point p1 can be written
as
p1 = p + ∆t!v
We now calculate a point p2 by stepping along the vector associated with p1 (calculated by the bilinear
interpolation algorithm of the section above), and continue in this fashion. This is illustrated in the following
figure.
!v p3
p2
p1
This figure suggests the following strategy: (1) Select a ”step size” ∆t; (2) Starting with the point pi , where
p0 = p, we perform the following iteration:
• Calculate !vi by bilinearly interpolating the vectors at the corners of the cell (Note: You first have to
determine the cell in which pi lies.), and
This algorithm is called Euler’s Algorithm.2 It has an error estimate which is bounded by O(∆t2 )3 ,
which is not very good for the estimates that we require in visualization. The algorithm is mostly used as an
initial approximation method for the better algorithms discussed below.
!vi
pP
i+1
pi+1
pi
2
Leonhard Euler was a Swiss mathematician of the 18th century.
3
This estimate comes from Taylor’s theorem which can be found in any elementary numerical analysis book.
In this figure, we begin at the point pi . We use Euler’s method to predict the location of the next point in the
iteration pPi+1 . We then use the calculated vector !vi+1
P at pP as additional information that can be used to
i+1
better predict pi+1 . To obtain pi+1 , we begin at pi and move half the step size (12 ∆t) along !vi , then half the
P . Mathematically, this can be written as
step size along the vector !vi+1
1 1
pi+1 = pi + ∆t!vi + ∆t!vi+1P
2 2
1 ! "
= pi + ∆t !vi + !vi+1
P
2
This numerical method is called a predictor-corrector method (for obvious reasons). It has an error
estimate which is bounded by O(∆t3 ), which is substantially better than Euler’s algorithm. The following
illustration shows the difference between Euler’s method and the improved Euler method on our example
cell.
The white dots are the approximation due to Euler’s method. The black dots are the approximation due to
the improved Euler method. This “improved” version is substantially used in the visualization community,
but a better version, generated by two German mathematicians, Runge and Kutta, is most frequently used.
1
!vi+1
!vi
2
!vi+1
3
!vi+1
p1i+1
p2i+1
p3i+1
pi
3
• !vi+1 – the vector corresponding to the point pi + ∆t!vi+1
2 .
We note that each predictor is used to obtain the subsequent predictor values. These vectors are blended into
the final result as follows:
1 1 1 1
pi+1 = pi + ∆t!vi + ∆t!vi+1
1
+ ∆t!vi+1
2
+ ∆t!vi+1
3
6 3 3 6
pi+1
1
6
∆t
1
6
∆t
1
3
∆t
pi
1 ! "
pi+1 = pi + ∆t !vi + 2!vi+1
1
+ 2!vi+1
2
+ !vi+1
3
6
The differences between Euler’s method, the improved Euler method, and the Runge-Kutta method for our
sample cell are shown in the following illustration.
The white dots are the approximations due to Euler’s method, the gray dots are the approximations due to
the improved Euler method, and the black dots are the approximations due to the fourth-order Runge-Kutta
method.
The Runge-Kutta method presented here has an error bounded by O(∆t4 ), which is significantly better
than the other two algorithms.
Summary
We have presented three numerical algorithms that can be used for particle tracing in vector fields. These
algorithms, the Euler method, the improved Euler method and the Runge-Kutta method, have errors bounded
by O(∆t2 ), O(∆t3 ), and O(∆t4 ), respectively. They are each easy to implement, and the user can pick the
algorithm necessary to achieve the errors desired by the application.
These algorithms lift easily to three-dimensional scalar fields. The only difference is that we must utilize
trilinear interpolation to calculate the vectors from the eight corner points of a cell.