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

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 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 ()

Linear function on its ownCode: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 ()

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()

Tweet This+ 1 thisPost To Linkedin