Reputation: 1354
So most of the time I am asking questions on stack exchange I usually have the wrong answer, but this time my code is producing the correct graphs I just want to know why. My question is why does theta update properly despite its dependence on omega coming after theta. By all means run my code if you do not believe me. Just as a warning I am not a computer scientist, I am just a physics student attempting to solve problems using computational methods, but I am interested why it works. Here is my code:
# This program is designed to show the difference between
# undamped, damped, and critically damped oscillators behavior over
# time.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
m = float(raw_input('Enter a mass for the pendelum '))
g = 9.8 # gravity
l = float(raw_input('Enter a number for length '))
theta = float(raw_input('Enter a number for postion (radians) '))
theta1 = float(raw_input('Enter a number for postion (radians) '))
theta2 = float(raw_input('Enter a number for postion (radians) '))
omega = 0
omega1 = 0
omega2 = 0
ArrayTheta = [theta]
ArrayTheta1 = [theta1]
ArrayTheta2 = [theta2]
q1 = float(raw_input('Enter a number for the damping constant '))
q2 = float(raw_input('Enter a number for the damping constant '))
q3 = float(raw_input('Enter a number for the damping constant '))
step = .001
tx = np.arange(0,5,step)
for i in np.arange(step,5,step):
theta = theta + omega*step
ArrayTheta.append(theta)
omega = omega + -(g/l)*np.sin(theta)*step+(-q1*omega*step)
theta1 = theta1 + omega1*step
ArrayTheta1.append(theta1)
omega1 = omega1 - (g/l)*np.sin(theta1)*step+(-q2*omega1*step)
theta2 = theta2 + omega2*step
ArrayTheta2.append(theta2)
omega2 = omega2 - (g/l)*np.sin(theta2)*step+(-q3*omega2*step)
# this does not really make sense to me that theta2 is able to update despite
# omega2 being after it. Clearly I should have put omega2 before theta2 yet
# it still works.
plt.plot(tx,ArrayTheta, color ='blue')
plt.plot(tx,ArrayTheta1, color ='red')
plt.plot(tx,ArrayTheta2, color ='green')
plt.ylabel('Position in Radians')
plt.xlabel('Time')
blue_patch = mpatches.Patch(color='blue', label='Damped q1')
red_patch = mpatches.Patch(color='red', label='Damped q2')
green_patch = mpatches.Patch(color='green', label='Damped q3')
plt.legend(handles=[blue_patch,red_patch,green_patch])
plt.show()
Upvotes: 0
Views: 63
Reputation: 33126
It looks like your thetas are angles, the omegas are their time derivatives, and that the time derivatives of the omega variables depend on the values of theta. You have a second order initial value problem.
Your code as-is works, but you are using naive Euler technique. This "works" because any n-dimensional second degree ODE can be re-expressed as a 2n-dimensional first degree ODE. This however discards a lot of the geometry of the problem.
You would be using the symplectic Euler technique if you change the code so that the derivatives are updated before the angles. This tends to preserve some of the geometry of the problem.
Upvotes: 0
Reputation: 373
The reason it is working is because you have defined the omega variables
omega = 0
omega1 = 0
omega2 = 0
So when you are at iteration 0 of your for loop, all the omega
values are equal to 0 before you calculate theta. Then that makes your theta
values as
theta = theta
theta2 = theta2
theta3 = theta3
You then update the value of omega
. In iteration 1, the omega
variables are equal to what you defined them as in iteration 0. Basically, when you update omega
in step n
, it is used to update theta
in step n+1
.
Upvotes: 0
Reputation: 42768
You set omega
, omega1
and omega2
to zero in the beginning and you update them in each for loop. You simply don't see, that the first two data points of each theta are the same.
Upvotes: 2
Reputation: 104752
All of your omega
variables are initialized to 0
before the loop starts, so on the first iteration they will have that value during the update of their respective theta
variables, so the first theta
value appended to the results will be the same as the initial value. The update to the omega
value changes the loop's behavior on the later iterations.
It might make more sense if you changed your code so that the ArrayTheta
lists start empty and your loop started at 0
. That way you'd only get the one copy of the initial theta
values in the lists, rather than the two your current code puts at the start:
omega = 0
omega1 = 0
omega2 = 0
ArrayTheta = [] # start these as empty lists
ArrayTheta1 = []
ArrayTheta2 = []
# unchanged stuff omitted
for i in np.arange(0,5,step): # start this loop at zero so you have the right number of steps
theta = theta + omega*step # on the first iteration, this doesn't change theta
ArrayTheta.append(theta) # so the first value in the list will be the original theta
omega = omega + -(g/l)*np.sin(theta)*step+(-q1*omega*step) # omega changes if theta != 0
Upvotes: 1