Trying to simulate a 1-dimensional wave

Question

I'm trying to make a simplified simulation of a 1-dimensional wave with a chain of harmonic oscillators. The differential equation for the position x[i] of the i-th oscillator (assuming equilibrium at x[i]=0) turns out to be - according to textbooks - this one:

m*x[i]''=-k(2*x[i]-x[i-1]-x[i+1])

(the derivative is wrt the time) so i tried to numerically compute the dynamics with the following algorithm. Here osc[i] is the i-th oscillator as an object with attributes loc (absolute location), vel (velocity), acc (acceleration) and bloc (equilibrium location), dt is the time increment:

for every i:
    osc[i].vel+=dt*osc[i].acc;
    osc[i].loc+=dt*osc[i].vel;
    osc[i].vel*=0.99;    
    osc[i].acc=-k*2*(osc[i].loc-osc[i].bloc);

    if(i!=0){
      osc[i].acc+=+k*(osc[i-1].loc-osc[i-1].bloc);
    }
    if(i!=N-1){
      osc[i].acc+=+k*(osc[i+1].loc-osc[i+1].bloc);
    }

You can see the algorithm in action here, just click to give an impulse to the 6-th oscillator. You can see that not only it doesn't generate a wave at all but it also generate a motion with growing total energy (even if I added the dampening!). What am I doing wrong?

Answer 1

All given answers provide a (great) deal of useful information. What you must do is:

• use ja72´s observation for performing the updates - you need to be coherent and compute the acceleration at a node based on the positions from the same iteration batch (i.e. do not mix $x⁽k⁾$ and $x⁽k+1⁾$ in the same expression of the acceleration as these are amounts belonging to different iteration steps)
• forget about your original positions as tom10 said - you need to consider only relative positions - the behaviour of this wave equation is like a laplacian smoothing filter applied to a polygonal line - a point is relaxed using its two directly connected neighbours, being pulled towards the middle of the segment determined by those neighbours
• finally, if you want your energy to be conserved (i.e. not have the water surface stop from vibrating), it's worth using a symplectic integrator. Symplectic/semi-implicit Euler will do. You don't need higher order integrators for this kind of simulation, but if you have the time, consider using Verlet or Ruth-Forest as they're both symplectic and of order 2 or 3 in terms of accuracy. Runge-Kutta will, like implicit Euler, draw energy (much slower) from the system, thus you'll get inherent damping. Explicit Euler will spell doom for you eventually - that's a bad choice - always (esp for stiff or undamped systems). There are tons of other integrators out there, so go ahead and experiment.

P.S. you did not implement true implicit Euler since that one requires simultaneously affecting all particles. For that, you must use a projected conjugate gradient method, derive Jacobians and solve sparse linear systems for your string of particles. Explicit or semi-implicit methods will get you the "follow the leader" behaviour you expect for an animation.

Now, if you want more, I implemented and tested this answer in SciLab (too lazy to program it in C++):

n=50;
for i=1:n
    x(i) = 1;
end

dt = 0.02;

k = 0.05;
x(20) = 1.1;
xold = x;
v = 0 * x;
plot(x, 'r');
plot(x*0,'r');
for iteration=1:10000
    for i = 1:n
        if (i>1 & i < n) then
            acc = k*(0.5*(xold(i-1) + xold(i+1)) - xold(i));
            v(i) = v(i) + dt * acc;
            x(i) = xold(i) + v(i) *dt;
        end
    end
    if (iteration/500  == round(iteration / 500) ) then
        plot(x - iteration/10000,'g');

    end
    xold = x;
end
plot(x,'b');

The wave evolution is seen in the stacked graphs below:

Answer 2

The growing amplitude over time is likely an artifact of the simple Euler integration you're using. You can try using a shorter timestep, or try switching to semi-implicit Euler, aka symplectic Euler, which has better energy conservation properties.

As for the odd propagation behavior (where the wave seems to propagate very slowly), it might be that you have too low of a spring constant (k) relative to the particle mass.

Also, the equation you posted looks slightly wrong as it should involve k/m, not k*m. That is, the equation should be

x[i]''=-k/m(2*x[i]-x[i-1]-x[i+1])

(c.f. this Wikipedia page). However this only affects the overall speed of propagation.

Answer 3

You've expressed the initial equation incorrectly in your code in two important ways:

First, note that the equation only expresses relative distortions; that is, in the equation, if x[i]==x[i-1]==x[i+1] then x"[i]=0 regardless of how far x[i] is from zero. In your code, you're measuring the absolute distance from an equilibrium for each particle. That is, the equation fixes the row of oscillators only at the boundary, like a single spring held at the ends, but you're simulating a whole set of little springs, each fixed at some "equilibrium" point.

Second, your terms like osc[i].acc+=+k*(osc[i-1].loc-osc[i-1].bloc); don't make much sense. Here you're setting osc[i].acc based solely on the absolute position of the particle next to it, and not on the relative position between the two.

This second problem is probably why you're gaining energy, but that's only a side effect of doing an incorrect simulation, not necessarily implying that your simulation is producing errors. Currently there's no evidence that you need to change from the simple Euler simulation (as Nathan suggested). First get the equation correct, and then get fancy with the simulation method if you need to.

Instead, write the equation for relative displacements, ie, with terms like osc.loc[i]-osc.loc[i-1].

comment: I rarely like to comment on others' answers to a question I've answered, but both Nathan and ja72 are focusing on things that just aren't the main issue. First get your simulation equations correct, then, maybe, worry about fancier approaches than Euler, and the order of updating your terms, etc, if you ever need to. For a simple, linear, first order equation, especially with a little damping, the direct Euler method works fine if the time step is small enough, so get it working like this first.

Answer 4

For one your acceleration osc[i].acc here depends on the updated position osc[i-1].loc and the not-updated yet position osc[i+1].loc.

You need to calculate all the accelerations first for all the nodes, and then in a separate loop update the positions and velocities.

It is best to create a function returning the accelerations given time, positions and velocities.

// calculate accelerations
// each spring has length L
// work out deflections of previous and next spring
for(i=1..N)
{
   prev_pos = if(i>1, pos[i-1], 0)
   next_pos = if(i<N, pos[i+1], i*L)
   prev_del = (pos[i]-prev_pos)-L
   next_del = (next_pos-pos[i])-L
   acc[i] = -(k/m)*(prev_del-next_del)
}

// calculate next step, with semi-implicit Euler
// step size is h
for(i=1..N)
{
    vel[i] = vel[i] + h*acc[i]
    pos[i] = pos[i] + h*vel[i]
}

You might need to use a higher order integrator scheme. Start by an 2nd order implicit Runge-Kutta.