G'MIC Adventure #6: Playing with the Wave Equation

Software G'MIC
1

Hey everyone!

Context.

It’s been a while since I wrote a new episode of #gmic-adventures, so I’ve decided to write a new one!

I’m going to try to show you how to use the G’MIC’s scripting language to implement an evolution of the Wave Equation, and see if we can create cool effects on images with that.

I don’t have a specific plan in mind, so I’ll be improvising a bit, but I’ll start by briefly describing what the wave equation is, how to discretize it, and how to implement it on a regular grid (2D to start with). Then we’ll see where that takes us.

Implementing the Wave Equation.

According to the wikipedia page,

The wave equation is a second-order linear partial differential equation for the description of waves or standing wave fields such as mechanical waves (e.g. water waves, sound waves and seismic waves) or electromagnetic waves (including light waves)

The equation itself is given by:

\frac{\partial^2 u}{\partial t^2} = c^2\;\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)

(I removed the z component in the sum because in our case, we’ll implement it only in a 2D space).

u is the 2D field that is subject to the wave equation. Think of it as an image of pointwise temperatures or as a height field, or whatever physical quantity you prefer :), but as an image that evolves over time. So basically, u depends on three variables: x,y (spatial coordinates) and t (evolution time).
Also in the equation, u is considered to be a continuous enough function, that’s why we can write derivatives (In fact, u must actually be C^², which means twice continuously differentiable, but that’s a math detail we won’t necessarily be concerned with here.)

Now, the question is: How do we “convert” this equation into a piece of code ? That seems a bit tricky at a first glance, but that’s actually not that hard and I’ll show you how to deal with it.

1. Discretization of second derivatives

In the Wave Equation, you mainly have second derivatives, one along the time dimension (left part of the equation, i.e. \frac{\partial^2 u}{\partial t^2}), and the others along the spatial dimensions (right part of the equation, i.e. \frac{\partial^2 u}{\partial x^2} and \frac{\partial^2 u}{\partial y^2}).

So, the first thing is we really need to understand how to estimate second derivatives from an algorithmic point of view. One solution actually comes from the Talor Expansion, which gives an approximation of a continuous function u, around a point x, assuming we also know the values of its N first derivatives at the same point x.

Suppose you want to evaluate the function u at a point that is very close to x (but not at x directly), e.g at x + dx, where dx is a small quantity that is known (and close to 0).
The Taylor expansion of u tells you that a good approximation of u(x + dx) is:

u(x + dx) \approx u(x) + dx\;\frac{\partial u}{\partial x}(x) + \frac{dx^2}{2}\;\frac{\partial^2 u}{\partial x^2}(x)

Now, if you do the same for a point located at x - dx, then:

u(x - dx) \approx u(x) - dx\;\frac{\partial u}{\partial x}(x) + \frac{dx^2}{2}\;\frac{\partial^2 u}{\partial x^2}(x)

And the trick is: summing these two equations gives:

u(x - dx) + u(x - dx) \approx 2\;u(x) + dx^2\;\frac{\partial^2 u}{\partial x^2}(x)

And thus you get one way to estimate the second derivative \frac{\partial^2 u}{\partial x^2} of the function u, if you know the values of u(x), u(x - dx) and u(x + dx):

\frac{\partial^2 u}{\partial x^2}(x) \approx \frac{1}{dx^2}\;\left(u(x - dx) + u(x - dx) - 2\;u(x)\right)

Now you get it. It’s relatively easy to compute the second derivatives of a function u along one of its variable.
In practice, u is defined on a discrete grid, so for estimating the spatial second derivatives \frac{\partial^2 u}{\partial x^2} and \frac{\partial^2 u}{\partial y^2}, we usually choose dx=1 and dy=1 (that’s not so close to 0, but that’s the minimal distance we can choose if we want to avoid spatial value interpolation). And for the variable t, we can choose a dt closer to 0, as we won’t have to actually store all the values of the u function for all possible t.

2. Putting all this together.

At the end of the day, it’s easy to show that a possible discrete approximation of the right part of the Wave Equation is:

c^2\;\left(\frac{\partial^2 u}{\partial x^2}(x,y) + \frac{\partial^2 u}{\partial y^2}(x,y)\right) \approx c^2\;\left( u(x+1,y) + u(x-1,y) + u(x,y+1) + u(x,y-1) - 4\;u(x,y) \right)

And same for the left part (except we keep dt as a parameter):

\frac{\partial^2 u^t}{\partial t^2}(x,y) = u^{t + dt}(x,y) + u^{t - dt}(x,y) - 2\;u^t(x,y)

Putting all this together leads to a discretized Wave Equation, the one we will actually implement:

u^{t+dt}(x,y) = \left(\frac{\partial^2 u^t}{\partial t^2}(x,y) + 2\;u^t(x,y)\right) - u^{t-dt}(x,y)

3. Stop talking, show me the code!

In G’MIC, a first, naive implementation of this equation will be:

foo :
  400,400,1,1x2 => u_t_1,u_t # Initialize 2D functions u() at time t-1 and t
  noise[u_t] 0.01,2 # Put some spikes in u^t
  repeat inf { # Evolution loop
    +f[u_t] " # Compute u() at time t+1
      const boundary = 1;
      const dt = 0.5;
      const c = 1;
      uxx = j(-1) + j(1) - 2*i; # Second derivative of u along x
      uyy = j(0,-1) + j(0,1) - 2*i; # Second derivative of u along y
      dt^2*c^2*(uxx + uyy) + 2*i - i(#$u_t_1)"
    w.
    wait 20

    # Prepare for next iteration.
    rm[u_t_1] =>[u_t,-1] u_t_1,u_t
  }

And the result, while being not impressive, looks like waves evolution anyway :wink: :

BTW, this code is not optimal and does not produce something very interesting right now. But improving all this will be the goal of the next post. Stay tuned!