This article is one in a series of articles about particle systems in Processing. To view the first of the series, click here.

In the last article, we discussed forces and left off with questions about our technique for determining the particle’s position at each frame. It turns out that it was wildly inaccurate. In this article, I am first going to try to explain the mathematics behind what we are doing and use that knowledge to expose the source of our problem. Then I am going to try to present a solution.

### Where Does the Problem Come From?

Recall the equation I presented in the last article which gives us distance given the state of a body and a time period:

This equation gives us an **exact** distance for the given acceleration, initial velocity, and time period and we used it to measure the error of our technique. So the question remains, why not use this equation to determine the particle position at each frame? Well, the problem is unfortunately too complicated for me to explain easily but I am going to try. From my understanding, that equation only works when acceleration is constant over time. If you take a hard look at it, this property will become apparent. The equation gives you a distance for a given **step** of time but is not actually a continuous function with respect to time. In other words, you cannot use it to solve acceleration for any given point in time. In a particle system, there may be **one to many** forces acting on a given particle at any time. On top of that, those forces are changing the acceleration of the particle in seemingly erratic ways. In fact, it is said that the function which determines the acceleration with respect to time is **well-defined** but is actually **unknown**. Therefore, we cannot get an exact value at any given point, we can only approximate it. Given that our calculations are **approximations**, we have to assume there will always be an error. All we can do is find techniques to make the error manageable.

### Solving Systems Governed by Unknown Functions

So, how exactly were we solving this unknown function before? Let’s step away from our particle system for a very brief moment. You may remember from high school maths that you can solve the integral of a function by integrating that function over an interval when given boundary conditions. You may also remember that in certain situations, such as when f(t, y) is unknown or cannot be solved analytically, you can solve these problems by taking an **initial value** and solving in steps. This is called an **initial value problem** and follows this form:

Some tricky manipulating of this formula results in a **series** which usually looks something like this:

So, I am pretty bad at all this maths o_0 but hopefully you are starting to see something here. In this series, we are determining the next value [y_n+1] using the current value [y_n] and adding the differential function. Just like we did with our particles, we found velocity by adding the last velocity plus the difference in velocity[acceleration]. Then we did the same for position. Acceleration is a **second order** derivative so when using this method which solves first order diffEQs, we must break it down into two first order derivatives:

where **v** is **velocity** and **p** is **position**.

Now let’s try and make sense out of this mess.

### Euler’s Method for Solving First Order Differential Equatons

You may not have known it, but the method we were using before in our code is known as Euler’s Method and is the defacto, inherent approach to this problem. The only thing I left out is the **time step**. The actual process should look something like this:

velocity.add(PVector.mult(acceleration, h)); location.add(PVector.mult(velocity, h));

where **h** is the time step in the simulation. As an astute reader, you may realize that the h was actually accounted for before. It was just assumed to be ’1′ ;) The common way to make the Euler method more accurate is to reduce the time-step [h]. I will explain why in a second, let’s just quickly alter the code to see what happens. Go check it out in action then come back and examine the source:

ArrayList particles; Particle a; float h = 0.2; //new time step < 1 void setup() { size(600, 600); stroke(0); strokeWeight(3); fill(150); smooth(); particles = new ArrayList(); a = new Particle(); a.location = new PVector(width/2f, height/2f); a.fixed = true; particles.add(a); } void draw() { background(255); for (int i = particles.size()-1; i >= 0; i--) { Particle p = (Particle)particles.get(i); // experiment w/ different friction coefficients applyDissipativeForce(p, 0.01); applyAttractiveForce(a, p, 5000f, 50f); if(!p.exist()) { particles.remove(i); } } } void mouseDragged() { particles.add(new Particle()); } void applyDissipativeForce(Particle p, float friction) { PVector f = PVector.mult(p.velocity, -friction); p.applyForce(f); } void applyAttractiveForce(Particle a, Particle b, float strength, float minDistance) { PVector dir = PVector.sub(a.location, b.location); float d = dir.mag(); if (d < minDistance) d = minDistance; dir.normalize(); float force = (strength * a.mass * b.mass) / (d * d); dir.mult(force); if (!b.fixed) b.applyForce(dir); if (!a.fixed) { dir.mult(-1f); a.applyForce(dir); } } public class Particle { PVector location; PVector velocity; PVector acceleration; float mass; boolean fixed; public Particle() { location = new PVector(mouseX, mouseY); //get velocity from direction and speed of mouse movement velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY); acceleration = new PVector(0, 0, 0); mass = 1; fixed = false; } public boolean exist() { if (fixed) { fill(255,0,0); } else { velocity.add(PVector.mult(acceleration, h)); location.add(PVector.mult(velocity, h)); acceleration.mult(0); fill(150); } ellipse(location.x, location.y, 10, 10); return true; } void applyForce(PVector force) { acceleration.add(PVector.div(force, mass)); } }

All I did was alter the integration part, added the **h** float at 0.2, and increased the strength to **5000** to speed things up a little. What you should have noticed is that it has gotten a lot slower and it should be obvious why: Our framerate is the same but we are stepping ahead in the simulation in smaller steps of time.

What you probably did not notice is the huge increase in accuracy. To prove this, let’s re-implement our python script experiment from the last article and compare it with the exact results again. This time we will have a time step-size of 0.1:

t = 0 h = 0.1 vel = 0 pos = 0 acc = 10 while (t <= 10): print "t = %s" % t print "pos = %s" % pos print "vel = %s" % vel print "-------------------" vel += acc * h pos += vel * h t += h

When we run this script we get [middle cut out for brevity]:

t = 0

pos = 0

vel = 0

——————-

t = 0.1

pos = 0.1

vel = 1.0

——————-

t = 0.2

pos = 0.3

vel = 2.0

——————-

t = 0.3

pos = 0.6

vel = 3.0

——————-

t = 0.4

pos = 1.0

vel = 4.0

——————-

t = 0.5

pos = 1.5

vel = 5.0

——————-

.

.

.

——————-

t = 9.6

pos = 465.6

vel = 96.0

——————-

t = 9.7

pos = 475.3

vel = 97.0

——————-

t = 9.8

pos = 485.1

vel = 98.0

——————-

t = 9.9

pos = 495.0

vel = 99.0

——————-

t = 10.0

pos = 505.0

vel = 100.0

**Position = 505**. You may remember that using these initial values with our equation:

yielded an exact end position of **500**. Now we are only **off by 5** whereas with a **time-step of 1** we were **off by 50**! This should show you that the error for this method is close to [i think?]:

There is a mathematical way to determine the real error but we aren’t going to go into it because we will be abandoning this method soon and it is too complicated for the little value it provides us. If you really want to know, check out the Wikipedia article on Euler’s method.

### Why does a Smaller h Yield More Accurate results?

The question arises, why is it more accurate? The answer may help you better understand what is happening. Let’s look at a graphic representation of the Euler method in action:

Pretend that the blue line is our unknown but well-defined function in which the acceleration is the dependent variable and that A0, A1, A2, A3, and A4 are the acceleration values we calculate. Essentially, at each step, we know the derivative [the tangent to the curve] and multiply it by the time step, then add it to the current **A** value to guess the next value. At that point, we assume we are at the correct spot and do this again and again and again. You can now see why the error accumulates. We are basically walking a dark path using a slow strobe light. Now imagine that the **A** values were taken in smaller intervals [the strobe light is blinking faster] and you can see that we would be closer to the blue line [exact solutions]. But if you took a bigger step size, you would be farther away from the exact values. The following pictures illustrate the point perfectly. The red line is the exact solution and the blue is the approximation:

Here is the source of these photos.

So as h gets smaller, our approximated curve inches closer to the exact curve.

### Picking a New Method

Where do we go from here? Well, we have determined that the Euler method is just not really feasible. Sure, it may have worked out for our little examples but we can achieve more stability and accuracy with some other methods. Plus we can’t pump up the framerate enough to get a realistic speed. It just doesn’t look believable. In the next article, I am going to show you how to implement the Runge-Kutta methods, particularly the Fourth Order Runge-Kutta method. I am also going to try to explain the differences between implicit and explicit integration.

### References

Gaffer on Games : Integration Basics

## 1 comment so far ↓

Laxman Bajagainon 03.23.12 at 8:11 pmDear Sir/Madam !

I have visited several pages of your your website. These are very useful and fruitful.

Thanks.

## Leave a Comment