Compare the Energy of Three Integrators

Goal

Discretize the same physical law in three ways and run them, and make the differences show up in a picture. When this lab is done, you can explain with evidence why a physics engine chooses a particular integrator.

Why it matters

The outward face of a physics engine is collision and friction, but what runs at every step underneath is the integrator. And the integration method does not merely introduce error; it changes the properties of the system. If you run a spring with explicit Euler, the energy grows by exactly (1 + dt²) times at every step and diverges within a few seconds. It is not that the force calculation is wrong; the order is wrong.

Semi-implicit Euler, which only swaps the order of two lines, is stable at the same cost. That is why game physics engines use it, and this fact does not fade once you pull out the numbers yourself and see them in a graph. So this lab saves the results as CSV and overlays them on a single PNG.

Steps

Put the toolbox in /root/integrator.
/root/integrator/out/euler.csv — explicit Euler.
/root/integrator/out/semi.csv — semi-implicit Euler.
/root/integrator/out/verlet.csv — velocity Verlet.
/root/integrator/out/energy.png — the three energy curves on one image.
/root/integrator/out/06-dt.txt — growth rates with different dt.
/root/integrator/out/damped.csv and damped.png — damping as a force.

Notes

Start a server with nohup python3 -m http.server 8080 -d /root/integrator/out & and open http://localhost:8080/energy.png in the Web preview.
The CSV is one header line + 401 lines. Do not leave out step 0 (the initial state).
Common mistake 1: writing explicit Euler as two lines and using the new x in the velocity calculation. At that moment it becomes semi-implicit Euler, and the two CSVs become identical.
Common mistake 2: not computing the initial acceleration before the loop in Verlet. The first step goes wrong.

Put the drawing toolbox in place

Save /root/integrator/gfxlib.py exactly as in the example, and use /root/integrator/check.py to draw a test pattern and make /root/integrator/out/00-check.png. The pattern is a 64x64 black background with a white (255,255,255) diagonal line from (0,0) to (63,63), and over it a red (255,0,0) horizontal line from (0,32) to (63,32).

From this lab on, you do not rebuild the PNG encoder. We hand you the same code you built by hand in the first lab as a tool — because file formats are not what you learn here.

The lab Pod has no volume, so the files you made in the previous lab are not kept. That is why each lab starts by putting the toolbox in place again.

Create Canvas(w, h, bg), draw the two lines with line(x0, y0, x1, y1, rgb), and then save with write_png(path). Draw the horizontal line later, so that the intersection (32,32) becomes red.

In this lab you use it to draw the energy curves. Divergence is hard to see from the numbers alone.

Explicit Euler

Use /root/integrator/euler.py to run a spring (mass 1, stiffness 1, acceleration a = -x) for 400 steps with dt=0.05 starting from x=1.0, v=0.0, and write to /root/integrator/out/euler.csv the header step,x,v,energy and 401 lines. The energy is 0.5*v*v + 0.5*x*x.

Explicit Euler uses only the values at the start of the step. Think of it as updating position and velocity at the same time.

x, v = x + dt * v, v - dt * x    # 오른쪽은 둘 다 옛 값

If you write it on one line, you will not make the mistake. If you split it into two lines, it is easy to use the new x in the v calculation, and then it becomes semi-implicit Euler.

You must also write step 0 (the initial state) as a line to get 401 lines. The final energy should be about 2.7 times the initial one — that is what this step is meant to show.

Semi-implicit Euler

Use /root/integrator/semi.py to run the same conditions by updating the velocity first and then updating the position with that new velocity, and make /root/integrator/out/semi.csv. The format is the same as in the previous step.

Only the order changes.

v = v - dt * x        # 속도 먼저
x = x + dt * v        # 방금 구한 새 속도를 쓴다

You can tell the two methods apart by the value at the first step. In explicit Euler x stays 1.0, but in semi-implicit it becomes 1 - dt*dt = 0.9975. This one digit is the whole difference between the two methods.

The energy should only rise and fall within a small range and not grow. The computational cost is the same as explicit Euler, yet the result is completely different.

Velocity Verlet

Use /root/integrator/verlet.py to run the same conditions with velocity Verlet and make /root/integrator/out/verlet.csv. One step is x += v*dt + 0.5*a*dt*dt, a_new = -x, v += 0.5*(a + a_new)*dt, and you replace a with a_new and go on to the next step.

You use the acceleration twice per step. The key is to update the velocity by averaging the a at the start of the step and the a_new after the position has moved.

Before entering the loop, you must compute the initial acceleration with a = -x. If you leave this out, the first step is wrong.

The first step's x comes out as 1 - 0.5*0.0025 = 0.99875. The three methods give different values from the first step, so you can check here.

The energy should stay almost unchanged — even after 400 steps it stays within a thousandth of the initial value.

Three curves on one image

Use /root/integrator/plot.py to draw the energy of the three CSVs on /root/integrator/out/energy.png (256x256, black background). The point for step i is px = int(i*255/400) and py = 255 - int(E/2.0*255) (clip anything outside 0–255), and the colors are Euler (255,80,80), semi-implicit (80,255,80) and Verlet (80,120,255). The drawing order is also this order.

The y axis puts energy 0 at the very bottom of the screen (255) and energy 2.0 at the very top (0). The initial energy 0.5 is plotted near 255 - int(0.5/2*255) = 192.

You do not need to connect the lines; plotting only points is fine. The 401 points fall into 256 columns, so they look naturally connected.

Read the CSV with open(path).read().splitlines(), skip the first line (the header), and then split on commas.

Only the red curve should shoot upward, and the green and blue should stay low. That picture is the conclusion of this lab.

How much worse does it get when you change dt

Run explicit Euler for the same simulated time T=20 with dt=0.01 (2000 steps), dt=0.05 (400 steps) and dt=0.2 (100 steps), and write how many times the initial value the final energy is to /root/integrator/out/06-dt.txt as three lines, dt001=, dt005= and dt02= (six decimal places).

The number of steps is T / dt. All three cases imitate the same 20 seconds, but the results are entirely different.

For this system there is a closed-form expression. At every step the energy is multiplied by exactly (1 + dt²), so after n steps it is (1 + dt²)^n times. Check whether the values you ran directly match this expression — if they match, it is strong evidence that your implementation is correct.

You only made dt four times larger, yet the growth rate becomes dozens of times larger. That is why a physics engine chooses to run several steps instead of making dt large.

Put damping in as a force

Add damping to semi-implicit Euler as a force (a = -x - 0.3*v, where v is the value before the update), run 400 steps with dt=0.05, and make /root/integrator/out/damped.csv (same format as before) and /root/integrator/out/damped.png (the energy curve, the same coordinate transform as in step 5, color (255,200,80)).

One step is like this.

a = -x - 0.3 * v      # v 는 아직 갱신하기 전 값
v = v + dt * a
x = x + dt * v

The energy should decrease without ever increasing. After 400 steps it drops below one percent of the initial value.

Compare it with the common method of multiplying the velocity by 0.99 at every step. That method is simple, but the amount of damping depends on dt, so objects move differently on machines with different frame rates. If you put it in as a force, it imitates the same physics even when dt differs.