I am writing code to solve the simple harmonic oscillator system using the Euler Method. The second order ODE for the system is given as two first order ODEs, x' = v and v' = -k/m x. The question says to solve the pair of equations for x and v as functions of time and plot x vs t. The model being used is a sodium atom so the mass is 3.82x10^-26 kg and k = 12.2 N/m. I am told to use an initial position of 1.0x10^-10 m and initial velocity v = 0 m/s.

I had a previous programme for solving the SHO system for a helical spring which worked perfectly. However the numbers for that system were much larger. In my code I just changed the values for the parameters and left everything else as it was for the helical spring. The graph should look like a sinusoidal curve but it doesn't and I keep getting overflow warnings.

My question is, what else do I need to change for this code to work for these numbers?

Python code for Simple Harmonic Oscillator

Graph

CodePudding user response：

Having tried out your code and made some changes to parameters, I believe that your problem is caused by the size of your time step.

If you consider the acceleration of that sodium atom when acted upon by the spring constant, it will be very large, around 10^16 ms^-2, so the atom will be moving very quickly around its equilibrium. The period of oscillation is of the order of \sqrt{m/k}, which in this case is about 1e-27 seconds.

So, if you use a time step of 0.1 seconds the atom will have oscillated a very large number of times and so you won't see anything looking like a sine wave. The overflow occurs because the distance and velocity are being estimated with a step too large for the Euler approximation to be appropriate.

I changed your code by setting:

t_max = 1e-12
delta_t = 1e-15

and got a nice sine wave out of it.