2
\$\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\$
4
  • \$\begingroup\$ In what sense is this the "Euler" method? It just looks like a discrete-timestep simulation to me. \$\endgroup\$ May 29, 2017 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\$ May 29, 2017 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, 2017 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\$ May 29, 2017 at 5:50

2 Answers 2

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\$
4
  • \$\begingroup\$ Thanks @Taylor. I would definitely take that into consideration, and I should have updated x first. \$\endgroup\$ May 29, 2017 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\$ May 29, 2017 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\$ May 29, 2017 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\$ May 29, 2017 at 11:22
0
\$\begingroup\$

This is well working code, although it's not Euler algorithm - it's Euler-Cromer (if you want to make sure it's showing cos function just make the maximum t in your while more 2 or 3 will do the thing)

Euler method uses current velocity for the next x, so you need to change the places of those

v=v-(10.0/1.0)*x*h with x=x+v*h

\$\endgroup\$
1
  • \$\begingroup\$ Welcome to Code Review! You typed "cos" - is that intended to be "because", or (trigonometric function) cosine? \$\endgroup\$ Dec 13, 2021 at 16:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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