python – odeint for time varying excitation

python – odeint for time varying excitation

The type of A=1e8*np.sin(np.log10(t)) is a numpy.ndarray, because t is. You can see this yourself if you evaluate type(1e8*np.sin(np.log10(t)))

But the problem is that the function dpdt expects A to be a scalar (float).

You can pull the definition of A inside of dpdt to fix this, be redefining it as:

def dpdt(y,t,Wb,V0,kT):
    y1=y[0]
    y2=y[1]
    A = 1e8*np.sin(np.log10(t))
    W1=(Wb+0.02*A)*1e-26      
    W2=(Wb-0.02*A)*1e-26      

    w_12=V0*np.exp(-W1/kT)   
    w_21=V0*np.exp(-W2/kT)
    
    
    dP1_dt=- w_12*y1 + w_21*y2
    dP2_dt=+ w_12*y1 - w_21*y2
    dp_dt=[dP1_dt,dP2_dt]     
    return dp_dt

Of course, we then must modify the call to odeint to remove the A parameter as follows:

p=odeint(dpdt,[y1_0,1-y1_0],t,args=(Wb,V0,k*T))

Another approach is to write dpdt so that it expects a argument A, which is actually a function of t. E.g.

def dpdt(y,t,Wb,V0, A,kT):
    y1=y[0]
    y2=y[1]
    W1=(Wb+0.02*A(t))*1e-26      
    W2=(Wb-0.02*A(t))*1e-26      

    w_12=V0*np.exp(-W1/kT)   
    w_21=V0*np.exp(-W2/kT)
    
    
    dP1_dt=- w_12*y1 + w_21*y2
    dP2_dt=+ w_12*y1 - w_21*y2
    dp_dt=[dP1_dt,dP2_dt]     
    return dp_dt

Now you define a function A, and call odeint as follows:

def A(t):
    return 1e8*np.sin(np.log10(t))

p=odeint(dpdt,[y1_0,1-y1_0],t,args=(Wb,V0,A, k*T))

python – odeint for time varying excitation

Leave a Reply

Your email address will not be published.