1
\$\begingroup\$

I have written very simple Python code to solve the simple harmonic oscillator using Euler method, but I am not sure if the program is correct or not. I would be very grateful if anyone can look at my code and suggest further improvements since I am very new to this programming thing.

import matplotlib.pyplot as plt
v=0.0   #initial velocity
h=0.01  #time step
x=5.0   #initial position

t=0.0
ta,xa=[],[]

while t<1.0:
    ta.append(t)
    xa.append(x)

    v=v-(10.0/1.0)*x*h    #k=10.0, m=1.0
    x=x+v*h
    t=t+h
plt.figure()
plt.plot(ta,xa,'--')
plt.xlabel('$t(s)$')
plt.ylabel('$x(m)$')
plt.show()
\$\endgroup\$
  • \$\begingroup\$ In what sense is this the "Euler" method? It just looks like a discrete-timestep simulation to me. \$\endgroup\$ – 200_success May 29 '17 at 2:09
  • \$\begingroup\$ For dy/dx=f(x,y), we can find the value of y at x1=x0+h as y(x1)=y(x0)+hf(x0,y0)$. Isn't this what the Euler's method precisely about ? \$\endgroup\$ – Roshan Shrestha May 29 '17 at 2:22
  • \$\begingroup\$ I don't see what role does g play. Besides, h is really dt. That said, it is far from harmonic oscillator. \$\endgroup\$ – vnp May 29 '17 at 4:37
  • \$\begingroup\$ Yes, it seems I did a mistake by calling variable g which is not needed for this one since I was doing free fall motion earlier. So, haven't I used the Euler method correctly in here. \$\endgroup\$ – Roshan Shrestha May 29 '17 at 5:50
2
\$\begingroup\$

The code is quite monolithic: calculation and presentation are separate concerns, and while I appreciate the context given by showing the presentation, it would be better to make clear that they're separate concerns by separating them with whitespace and commenting briefly.

Also on structure: separate constants (h) from variables (v, x, t).


    v=v-(10.0/1.0)*x*h    #k=10.0, m=1.0
    x=x+v*h
    t=t+h

Two things: firstly, what are k and m? It seems that they should be constants with comments explaining their physical significance.

Secondly, this isn't Euler's method. Euler's method is \$\vec{x}_{n+1} = \vec{x}_n + hf(t_n, \vec{x}_n)\$ where \$f = \frac{d\vec{x}}{dt}\$. Here \$\vec{x} = (v, x)\$ and \$f(t, v, x) = (\frac{dv}{dt}, \frac{dx}{dt}) = (-\frac{k}{m}x, v)\$. In other words, to be Euler's method you should update x according to the previous value of v, not the updated value of v. Since you're using Python, you can take advantage of simultaneous assignment:

    v,x=v-(k/m)*x*h,x+v*h
    t=t+h

(As it happens your buggy implementation works better than Euler's method, but if it was intended to implement Euler's method then it's still technically buggy).

\$\endgroup\$
  • \$\begingroup\$ Thanks @Taylor. I would definitely take that into consideration, and I should have updated x first. \$\endgroup\$ – Roshan Shrestha May 29 '17 at 10:34
  • \$\begingroup\$ @RoshanShrestha, no, that would be equally wrong. You should update them simultaneously, which depending on the language may mean using temporary variables. \$\endgroup\$ – Peter Taylor May 29 '17 at 10:52
  • \$\begingroup\$ But, I have used temporary variables for x,v and t outside the while loop. So, does it mean that if I order my code for Euler method x=x+vh v=v-(10.0/1.0)*xh will still be wrong ? \$\endgroup\$ – Roshan Shrestha May 29 '17 at 10:56
  • \$\begingroup\$ @RoshanShrestha, I wouldn't call those temporary variables. Also, I've realised that the code is Python (I had initially misread and thought it was Matlab, which I don't know), so that means it supports simultaneous assignment. I've edited sample code into the answer. \$\endgroup\$ – Peter Taylor May 29 '17 at 11:22

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.