# How to improve and stabilize this code simulating a damped mass-spring oscillator? Runge-Kutta?

I wrote the following function which is simulating a damped mass-spring oscillator. It is being driven by the audio sample input at 44.1 kHz sampling to create the same effect as a resonant bandpass filter.

The "base" is being directly positioned by the audio input samples (input), then I calculate the spring length, apply the forces to the mass, and solve the mass position (output).

  double processNextSample(double sampleInput) {

basePosition = sampleInput;

//solve spring length
currentSpringLength = massPosition - basePosition;
deltaSpringLength = currentSpringLength - springLengthAtRest;

//solve spring force
double springForceOnMass = -springK * deltaSpringLength - massVelocity * dampCoeff;

//solve mass position
massVelocity += springForceOnMass * deltaTime; //mass = 1 kg
massPosition += massVelocity * deltaTime;

//return output
return (massPosition - springLengthAtRest);

}

double springLengthAtRest = 1000; //unnecessary but just to help me visualize it
double currentSpringLength = 1000;
double deltaSpringLength = 0;
double basePosition = 0;
double massPosition = springLengthAtRest;
double massVelocity = 0;


In theory, at 44.1 kHz sample rate, one can use resonant bandpasses up to 20 kHz (or up to 1/2 Nyquist) with no stability problem, so I presume I should be able to resonate a damped mass-spring oscillator up to the same frequency if done correctly.

However, I can't get anywhere near the high frequencies I'm trying to without losing stability. It works fine at low frequencies.

I don't have any particular math or coding knowledge though of where the problem is or how to fix it. I believe one part of the problem is that I am told I am using the explicit Euler method and this may be unstable.

I have read about the Runge Kutta RK4 but I don't understand how to use it. Wouldn't that also only help me replace the solve mass position portion, and wouldn't that desynchronize it somehow from the rest?

If you are kind enough to answer, please speak to me like someone who has no math past first year university stats, no physics past a physics 100 course, and self taught for coding. Any sample code is appreciated to illustrate. Thanks. :)

• How long does a single run take? If not too long try decreasing daltaTime. Stability of these explicit time schemes are often related to the time step size. Jan 24 at 12:15
• I should have said I don't want to increase the sample rate. Sorry. Jan 24 at 17:27
• Numerical algorithms like this are very very sensitive to step size. If that can't change, then you could maybe rescale the problem so that you are operating in some characteristic time scale that is easier to simulate Jan 24 at 18:18

• +1. Nyquist theorem doesn't say much about whether $x_{n+1} = x_{n} + hf(x_n)$ converges, which is the real issue here. Jan 25 at 15:43