Wave Equation Simulation

I've been recently trying to simulate a wave equation in P5.js.

My aproach is to plot a bunch of points with $$(x, y)$$ coordinates. I've used the wave equation for obtaining the movement of those particles.

This is my numerical integration method:

var pos     = this.masses[i].pos
var nextPos = this.masses[i + 1].pos
var lastPos = this.masses[i - 1].pos

var vel = this.masses[i].vel, acc = this.masses[i].acc // for simplicity

acc.x = 0
acc.y = c * c * ((nextPos.y - 2*pos.y + lastPos.y)) / dxdx

vel.x += acc.x * dt
vel.y += acc.y * dt

pos.x += vel.x * dt
pos.y += vel.y * dt

Edit: variables definition.

dxdx = 1e-5 //small step
dt = deltaTime // value from P5js (calculates time between frames)
c = 10

My acceleration is defined like this: $$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$$ So I convert my equation into a discrete equation like this: $$acc = c^2 \frac{u(x + 2\Delta x, t) - 2 u(x + \Delta x, t) + u(x, t)}{\Delta x^2}$$ As shown here in my code:

acc.y = c * c * ((nextPos.y - 2*pos.y + lastPos.y)) / dxdx

And then simply I integrate by the Euler's method.

But my result is quite strange: As you can see, the wave grows, when the expected result was something like this: What I'm doing wrong? Someone could help me? Is there something I'm missing? Thanks :)

• You did not write the values of delta x, delta y, delta t, and c. That's a piece of key information. Nov 13 '21 at 10:48
• @GiorgioP Now is changed. Thanks! Nov 13 '21 at 10:52
• Ok, but, even if it is set by P5js, it is important to know what is the numerical value of deltaTime. Could you get it from the code producing the wrong picture? Nov 13 '21 at 13:24
• I mean, dt is just a value that depends on my framerate. The dt value is just the inverse of the framerate. So if my computer is always running the simulation at $60$ frames per second, dt would be aproximately $0.01666666$ Nov 13 '21 at 13:50
• Have you analyzed stability of your scheme? Forward Euler's method is not always stable. Nov 13 '21 at 14:00

The numerical integration of hyperbolic partial differential equations (PDE) is apparently a direct extension of the usual discretization algorithms used in the case of ordinary differential equations (ODE) and of the related stability analysis. However, the analogy is not complete. The mathematical theory behind PDE introduces some additional constraints, at least in the case of explicit methods, that are not present in the case of ODE.

This is the case of the so-called Courant–Friedrichs–Lewy (CFL) condition. The physical motivation and general description of the condition go as follows.

Hyperbolic PDE's are characterized by a set of conic-like hypersurface passing through each space point (a couple of lines, in the simple case of one space and one time variables) determining which points in the past can causally influence the given point and which points in the future can be influenced by the same point. It is the familiar light-cone concept of Special Relativity, whose origin is actually in the wave equation controlling the light behavior. Courant, Friedrichs, and Lewy were able to show that a necessary condition for the convergence of an explicit numerical method for hyperbolic PDE is that the physical light-cone must be contained in the numerical light cone.

In the case of your scheme, the condition requires that. $$\frac{c \Delta t}{\Delta x} \leq 1,$$ which is not the case, with the data you have reported.

Therefore, if you play with the simulation parameters ($$\Delta t$$ and $$\Delta x$$) to satisfy the CFL condition, it should be possible to get the correct numerical behavior of the one-dimensional wave.

I am not at all an expert in numerical methods but the stability condition for the algorithm "time stepping" is $$c<∆x/∆t$$.

This condition does not seem to be verified in your case, with $$c=10$$ ?