Adam Warner
Adam Warner

Reputation: 1354

Why Does My Code Update

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

Answers (4)

David Hammen
David Hammen

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

coffeepls
coffeepls

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

Daniel
Daniel

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

Blckknght
Blckknght

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

Related Questions