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

Question

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) { 

      //adjust base position
      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. :)

Answer 1

I don't think the Nyquist theorem is applicable for the timesteps of direct numerical simulation. I think it would only be useful if you already had the full trajectory of your mass and would do a fourier analysis at a given sampling rate. Since you don't have the actual trajectory yet you might have to do a lot more number crunching.

you should do a back-of-the-napkin estimation of the needed stepsize. Lets say your driving force is at 44.1 kHz. That puts one period of oscialltion in the ballpark of 2e-4 seconds. In order to resolve the forces on your mass correctly I would estimate that you need between 10 and 100 samples per period. That puts your stepsize to be about:

2e-6s to 2e-7s

I also think that switching to RK4 will not allow you to increase that significantly.

Answer 2

the very simple answer is to use an ode solver. There are decades of research on how to efficiently solve ODEs and there are very high quality open source solvers written by groups of people over decades who have made efficiently solving ODEs their job.