An attempt to simulate a one-dimensional wave

Question

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

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

(the derivative is time-dependent), so I tried to calculate the dynamics using the following algorithm. Here osc[i] is the ith oscillator as an object with attributes loc (absolute location), vel (speed), acc (acceleration) and bloc (equilibrium position), 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 6th oscillator. You can see that not only does it not generate waves at all, but it also generates movement with growing full energy (even if I added humidification!). What am I doing wrong?

Answer 1

All provided answers provide (great) helpful information. You should do the following:

• use observation ja72 to perform updates - you need to be consistent and calculate the acceleration on node based on positions from the same batch of iterations (i.e. don't mix $ x⁽k⁾ $ and $ x ⁽k + 1⁾ $ in the same expression acceleration, since these are amounts belonging to different steps of the iteration)
• forget about your starting positions, which were discussed in Volume 10 - you only need to consider relative positions. The behavior of this wave equation is similar to a Laplacian smoothing filter applied to a broken line - the point is weakened using two directly connected neighbors, being pulled to the middle of the segment defined by these neighbors.
• finally, if you want your energy to be saved (that is, it does not stop on the surface of the water from vibration), it is worth using a symplectic integrator. A symplectic / semi-implicit Euler will do. You do not need higher-order integrators for this kind of simulation, but if you have time, consider using Verlet or Ruth-Forest, because they are both symplectic and of order 2 or 3 in terms of accuracy. Runge-Kutta will, like the implicit Euler, gain energy (much slower) from the system, so you get internal attenuation. An explicit Euler will ultimately doom your death - it's a bad choice - always (esp for rigid or non-fading systems). There are many other integrators out there, so keep experimenting.

PS , you did not implement the true implicit Euler, as this requires simultaneous exposure to all particles. To do this, you must use the projected conjugate gradient method, derive the Jacobian, and resolve the sparse linear systems for your particle string. Explicit or semi-implicit methods will allow you to "follow the leader" that you expect from the animation.

Now, if you want more, I have 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 shown in the stacked graphs below:

Answer 2

The increasing amplitude in time is probably an artifact of the simple Euler integration you use. You can try to use a shorter time interval or try switching to the half-line Euler , also known as the symplectic Euler, which has better energy-saving properties.

As for the behavior of the odd propagation (where the wave apparently propagates very slowly), you may have too little spring (k) constant relative to the mass of the particle.

Also, the equation you got looks a bit wrong, since it should include k / m, not k * m. That is, the equation should be

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

(cf this wikipedia page ). However, this only affects the overall propagation speed.

Answer 3

You incorrectly formulated the original equation in your code in two important ways :

First, note that the equation expresses only relative distortions; that is, in the equation, if x[i]==x[i-1]==x[i+1] then x"[i]=0 no matter how far x[i] from zero. In your code you measure the absolute distance from the equilibrium for each particle, that is, the equation captures a series of oscillators only at the boundary, like a single spring held at the ends, but you simulate a whole set of small springs, each of which is fixed at some "equilibrium" point.

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

This second problem is probably why you are gaining energy, but this is only a side effect of incorrect modeling, not necessarily implying that your simulation is causing errors. There is currently no evidence that you need to go from simple Euler modeling (as Nathan suggested). First get the right equation, and then come up with a modeling method if you need to.

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

Comment: I rarely want to comment on the answers of others to the question I answered, but both Nathan and ja72 focus on things that simply are not the main problem. First get the right simulation equations, and then maybe worry about more convenient approaches than Euler, as well as how to update your conditions, etc., if you ever need to. For a simple, linear first-order equation, especially with a slight attenuation, the direct Euler method works great if the time step is small enough, so work it out first.

Answer 4

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

You need to first calculate all the accelerations for all nodes, and then in a separate cycle to update the position and speed.

It’s best to create a function that returns accelerations based on time, positions, and speeds.

 // 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 may need to use a higher order integrator circuit. Start with the implicit 2nd order Runge-Kutta .

Answer 5

There is a simple and easy way to do this described here. its only one line of code and one loop and very realistic.

http://users.softlab.ece.ntua.gr/~ttsiod/wavePhysics.html#

however, I cannot establish the boundary conditions for the wave to bounce off the end points without inverting.

the best I can do is set the last to the previous value, that is, the value vel = 0 at the free end, but then the wave just ends and does not leave the raw ones back as I wish.

Does anyone understand this?