odeint для переменного во времени возбуждения

Я хочу решить связанное дифференциальное уравнение. введите здесь описание изображения

Здесь A — возбуждение, а y(t) — реакция, где y=[y1,y2]. Я написал следующий код, который работает с постоянным значением одиночного возбуждения. Но это не удается, когда я ставлю переменное во времени возбуждение. Другими словами, мой код работает для A=0.5e8 Не работает для A=1e8*np.sin(np.log10(t)). Как я могу это исправить?

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint

k=1.381e-23
T=300
V0=1e13


def dpdt(y,t,Wb,V0,A,kT):
    y1=y[0]
    y2=y[1]
    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
    
y1_0=0.5
Wb=15936297.62 
t=np.logspace(-5,5,100)

A=1e8*np.sin(np.log10(t))  #Doesn't work
A=0.5e8  #Works

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

plt.semilogx(t,p[:,0])
plt.grid()
plt.show()

python odeint
person Alam    schedule 13.11.2020    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Тип A=1e8*np.sin(np.log10(t)) — это numpy.ndarray, потому что t есть. Вы можете убедиться в этом сами, если оцените type(1e8*np.sin(np.log10(t)))

Но проблема в том, что функция dpdt ожидает, что A будет скаляром (float).

Вы можете вытащить определение A внутри dpdt, чтобы исправить это, переопределив его как:

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

Конечно, затем мы должны изменить вызов odeint, чтобы удалить параметр A следующим образом:

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

Другой подход состоит в том, чтобы написать dpdt так, чтобы он ожидал аргумент A, который на самом деле является функцией t. Например.

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

Теперь вы определяете функцию A и вызываете odeint следующим образом:

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))
person L.Grozinger    schedule 13.11.2020