1. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Dec 2012
Posts
5
Rep Power
0

Hi guys,

So I've been stuck on this problem for some time now and I've decided to ask for your help. I am attempting to plot a damped oscillation behaviour.

I have two functions which interact with each other. One being linear and the other being sinusoidal. The linear function is an amplitude controller of the sinusoidal, if you will. It is a way to to put the sinusoidal function to zero.

The sinusoidal function is an ordinary differential equation, and plots as it should on its own. odeint is used from scipy for this operation.

The linear function starts with an original value x and begins decreasing in value starting at time 100 and finally reaches 0 at time 200. This also plots as it should on its own.

The combined plot looks as it should until time 100, which is where it goes a bit crazy. What it should plot is a nice damped oscillation behaviour. Maintaining its period, wave length etc...

What it actually plots is a nice and consistent oscillation going absolutely nuts at time 100, and reaches 0 at time 200. While the amplitude decreases presumably as it should, I lose consistent wave behaviour, which is a nightmare for physicists.

So, my question is the following:

Is there a better way of making a function interact with another function that what I have done?

My first thought was that the linear function is being integrated because of the way I wrote the code...but that shouldn't make the plot go crazy. Then I thought it was an anti-aliasing problem, but changing the time step makes no difference.

Here's the actual code. Thank you so much for your time and I really appreciate your help.

Version 1

Code:
```from scipy.integrate import odeint
import pylab
import math

###############

th =  -1.47389
phi = -0.33449319
b = -200.0
tmax = 300.0
dt = .1
hy = 50.0

###############

t1 = 100.0
t2 = 200.0
tt = t2-t1
hx = -2*hy
h0 = math.sqrt((hy**2)+(hx**2))

###############

def h(t):
if t< t1:
return h0
elif t< t2:
h = (-h0/tt)*(t-t1) + h0
return h
else:
return 0.0

def f (dy, t):
th = dy [0]
phi = dy [1]
calc = [h(t)*math.sin (phi), b*math.cos (th)+ math.cos(phi)*h(t)*math.cos(th)/math.sin(th)]
return calc

dy = pylab.array([th, phi])
t = pylab.arange(0.0, tmax, dt)
y = odeint(f, dy, t)
yth = pylab.cos(y[:,0])
pylab.xlabel('time (t)')
pylab.ylabel('y = cos (theta)')
pylab.title('Theta vs Time')
pylab.grid(True)
pylab.plot (t, yth)
pylab.show ()```
Version 2 (separates x and y components)

Code:
```from scipy.integrate import odeint
import pylab
import math

###############

th =  -1.47389
phi = -0.33449319
b = -200.0
tmax = 300.0
dt = .1
hy = 50.0

###############

t1 = 100.0
t2 = 200.0
tt = t2-t1
hx = -2*hy
h = math.sqrt((hy**2)+(hx**2))

###############

def hyy(t):
if t< t1:
return hy
elif t< t2:
calc = (-hy/tt)*(t-t1) + hy
return calc
else:
return 0.0

def hxx(t):
if t< t1:
return hx
elif t< t2:
calc2 = (hx/tt)*(t-t1) - hx
return calc2
else:
return 0.0

def f (dy, t):
th = dy [0]
phi = dy [1]
calc3 = [math.sin (phi)*hxx(t) - math.cos(phi)*hyy(t), b*math.cos (th)+ (math.sin(phi)*hyy(t) + math.cos(phi)*hxx(t))*math.cos(th)/math.sin(th)]
return calc3

dy = pylab.array([th, phi])
t = pylab.arange(0.0, tmax, dt)
y = odeint(f, dy, t)
yth = pylab.cos(y[:,0])
pylab.xlabel('time (t)')
pylab.ylabel('y = cos (theta)')
pylab.title('Theta vs Time')
pylab.grid(True)
pylab.plot (t, yth)
pylab.show ()```
Linear function on its own

Code:
```import pylab
import math

t1 = 100.0
t2 = 200.0
tt = t2-t1
tmax = 300
dt = 0.1
nmax =int(tmax/dt)

def h(t):
h0 = 115.0
if t < t1:
return h0
elif t < t2:
h = ((-h0/tt)*(t-100)) + h0
return h
else:
return 0.0

harray=pylab.zeros(nmax)
t = pylab.arange(0.0,tmax,dt)
for j in range (0,nmax):
harray[j]=h(t[j])

pylab.plot(t,harray,'-')
pylab.show()```
2. I used gsl initial value ode. Graph of your version 1 time,calc[0] attached. It is not the same as what you expect.

Simple harmonic oscillator (spring) has
x''(t) = -k x(t)

where x(t) is displacement from unstressed position.
Your equation looks significantly more complicated.
3. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Dec 2012
Posts
5
Rep Power
0

My equation is complicated because it's describing a hamiltonian chaos in a single layer magnetic system. While I also have no idea what that actually means, it turns out that a photon oscillates when it interacts with a magnetic field.

But thanks a lot for that plot, that's actually one of the results I should be getting. With that I can prove that a photon has an elliptical polarization after it hits my magnet system.

Thank you so much!
4. I've already overwritten the code that solved your function. Recreation would be quick, but you'd need to install the gsl (gnu scientific library) and python 3, and I will support it (tonight) only on linux. So if you have all that, and want it, I'll tell explain.
5. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Dec 2012
Posts
5
Rep Power
0
That sounds awesome! But unfortunately my project leader insists that I use 2.7.3 because he hasn't gotten around using 3. Alas life has its kinks.

While I appreciate your help and am astonished at how fast you've solved my problem, any effective solution must be found in the old python.

But seriously, thank you. I'm currently rewriting the equation to make it a bit simpler. Hopefully this may solve the problem. I'll post up results once I finish.
6. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Dec 2012
Posts
5
Rep Power
0
I've messed up my plot even more, so I'm back to version 1 of the code. If you are willing to explain your solution in 2.7.3, I am at your mercy. If not, then thanks anyway.
7. I am not willing to explain in 2.7.3. I'm using an interface to the gsl I wrote for python 3.0. It core dumps in python 2. I tried it by mistake earlier tonight.

See attachment to post 4 at this link.

Are matplotlib, numpy, scipy available for python3? I use gnuplot for graphics.
Last edited by b49P23TIvg; December 19th, 2012 at 12:24 AM.
8. No Profile Picture
Registered User
Devshed Newbie (0 - 499 posts)

Join Date
Dec 2012
Posts
5
Rep Power
0
That's unfortunate, any clues you could give me on why it's not plotting as it should?
9. yth = pylab.cos(y[:,0])
I didn't take the cosine of anything to plot.
That won't help, it's different though.

Try some of the other methods found in

>>> help(scipy.integrate.ode)