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

    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?


  • 1
    $\begingroup$ What did you use for $m$, $k$ and the initial conditions as I don't want to look through that code - can you add those to the problem statement? $\endgroup$
    – Moo
    Sep 21, 2018 at 17:44
  • $\begingroup$ I used m =1 and k = 5, will they affect if the solution is unstable or not? $\endgroup$
    – ngbtwby
    Sep 21, 2018 at 18:06
  • $\begingroup$ Also the initial conditions are x(0)= -1 and x'(0) = 0 $\endgroup$
    – ngbtwby
    Sep 21, 2018 at 18:12
  • $\begingroup$ The solution to this ODE is not stable: think of it as a rock on the top of a hill, it seems stable but a little variation makes it crumble. So your numeric solution can easily accumulate rounding errors. According to how much time you have to simulate I'd add a little damping. $\endgroup$
    – N74
    Sep 21, 2018 at 18:18
  • 1
    $\begingroup$ All higher order numeric methods (than Euler) evaluate the function at more than one point, that's why you need to pass a functional and not a constant value. Anyway I'm happy it helps. $\endgroup$
    – N74
    Sep 21, 2018 at 19:29

1 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);


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.

  • $\begingroup$ Still this doesn't account for the marginal stability of the problem, as the OP observed in the next to last comment to the question. $\endgroup$
    – N74
    Sep 22, 2018 at 15:17
  • $\begingroup$ The code as it is ignores the Dopri results, overwrites them directly with the RK4 results. $0.02$ is a sensible step size for this problem and RK4, but far from the optimal $0.0005$. At time somewhere between $t=6$ to $t=8$ the step error will be magnified into visibility, giving an increasing phase and amplitude error. $\endgroup$ May 29, 2021 at 12:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.