Unstable solution for simple harmonic motion using dormand-prince and RK4 method

Question

I am trying to solve the differential equation for a simple harmonic motion in javascript that is given by

$$ mx'' = -kx $$

We know the solution is a sinusoidal. The problem I face is that the amplitude of the sinusoidals in my solution are increasing exponentially.

So the first method I tried was using the numeric.js, It uses the doramnd-prince method to solve the equation. The relevant code the library uses to solve this ODE is https://github.com/sloisel/numeric/blob/master/src/numeric.js#L2885-L2999

Next I tried this GMA1D library, I think it uses RK4 method. There relevant function that the library use to solve the ODE is here https://github.com/llarsen71/GMA1D/blob/gh-pages/javascripts/ODE.js#L89-L123

Here is my implementation

simulate(){
    var state = [this.x, this.dx]
    var time = this.T;

    //var F = 2*Math.sin((-2) - state[0])
    var F = -5*state[0];

    console.log(`Position x : ${this.x}, Velocity v : ${this.dx} at time t: ${this.T}`);
    // calculate new position one one is used the other is commented out

    // Method 1 - Dormand prince
    [this.x, this.dx ]= numeric.dopri(0,0.02,state,function(t,x){return [x[1], F]}).at(0.02);
    // Method 1 - END

    /// Method 2 RK4
    var V = function(t,pt){
      var x = pt[0], y =pt[1];
      return [y, F]
    }
    var ode_ = new ODE(V);
    var soln = ode_.solve(1,0.02,this.T,state);
    this.x = soln.pts[1][0]
    this.dx = soln.pts[1][1]
    // Method 2 - END

    this.T += 0.02;

  }

For both these libraries my solution becomes unstable as I iterate, I have tried to change the step size but that does not help. I have the following questions

1. I dont't have much experience with MATLAB ode45, but I would likely get a correct solution from MATLAB. What method does it use to solve the DE, if its dormand-prince why does numeric.js solution is unstable, are there different variationso of dormand-prince. If yes, which one should I use?

2. Can someone provide code for solving the above DE in MATLAB and plotting the result?

3. Can I expect a correct solution from RK4?

4. Which method should I use to solve Differential Equations?

5. While trying to figure out the problem I read about stiff and non-stiff equations, what are these types of equations and how that analysis is useful?

6. Point me to some good resources where I can learn more about all the above questions so next time I have a better understanding. :)

7. Does anyone know any good Javascript libraries to solve differential equations?

Gracias

Answer 1

This is really a javascript question.

The problem is that you return a global variable in the vector field instead of the function of the state passed as argument. Define the force F as proper function of a state vector argument,

var F = function(x) { return -5*x[0]; }

and use that in the ODE system functions,

[this.x, this.dx ]= numeric.dopri(0,0.02,state,function(t,x){return [x[1], F(x)]}).at(0.02);

and

var V = function(t,pt){
  var x = pt[0], y =pt[1];
  return [y, F(pt)]
}

Then the methods should have the corrects slopes at their interior stages and thus the correct error order which should then give a similar error order to the amplitude.