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

    Join Date
    Dec 2012
    Posts
    5
    Rep Power
    0

    Please Help! A Better way of having a function within a function?


    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. #2
  3. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,711
    Rep Power
    480
    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.
    Attached Images
    • File Type: png p.png (3.9 KB, 30 views)
    [code]Code tags[/code] are essential for python code and Makefiles!
  4. #3
  5. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Dec 2012
    Posts
    5
    Rep Power
    0
    Hey, thanks for your reply.

    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!
  6. #4
  7. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,711
    Rep Power
    480
    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.
    [code]Code tags[/code] are essential for python code and Makefiles!
  8. #5
  9. 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.
  10. #6
  11. 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.
  12. #7
  13. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,711
    Rep Power
    480
    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.
    [code]Code tags[/code] are essential for python code and Makefiles!
  14. #8
  15. 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?
  16. #9
  17. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,711
    Rep Power
    480
    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)
    [code]Code tags[/code] are essential for python code and Makefiles!

IMN logo majestic logo threadwatch logo seochat tools logo