Lang–Kobayashi方程描述为:
第一部分(Drive laser)是描述的驱动激光器,第二部分(Response laser)描述的是响应激光器。实验结构图如下:
虚线框表示响应激光器中的闭环配置。开环中响应激光器无反馈,产生的时间序列与驱动激光器相同,而闭环配置对于响应激光器添加了反馈,产生的时间序列与驱动激光器不同。
开环配置产生的混沌时间序列如下:
Figure 3开环配置混沌时间序列(蓝色为驱动激光器,红色为响应激光器,黑色为延迟驱动激光)
Figure 4开环配置三种混沌时间序列对比
开环配置产生的混沌时间序列如下:
开环配置中驱动与响应和延时之间的对比:
闭环配置产生的混沌时间序列:
闭环配置中驱动与响应和延时之间的对比:
混沌序列的横坐标为时间序列,点与点之间的间隔为定值。纵坐标为各中电场振幅的平方。混沌与噪声十分相似,又有很大不同:混沌是指非线性动力学系统的复杂演化行为,具有确定性的规律,但对初始条件极其敏感。噪声是指随机性的信号或干扰,以无规律的、随机的成分为特征,会对信号进行干扰和扰动。混沌系统的行为是由特定方程或规则决定的,而噪声是一个随机过程。混沌系统表现出的行为是非周期性的,而噪声通常不具备明确的周期性。
闭环配置数值模拟python代码:
import math
import matplotlib.pyplot as plt
C = 2.99792458e8
M = 6
r3Dri = 0.008
r3Res = 0.008
JRatio = 1.3
LDri = 0.6
LRes = 0.6
LInj = 1.2
detun = -4.0e9
kapInjRatio = 12.5
r2 = 0.556
tauIn = 8.0e-12
lambdaa = 1.537e-6
Gn = 8.4e-13
N0 = 1.4e24
tauP = 1.927e-12
tauS = 2.04e-9
alpha = 3.0
Nth = float
Jth = float
J = float
kapDri = float
kapRes = float
kapInj = float
tauDri = float
tauRes = float
tauInj = float
omega0 = float
phaseShiftDri = float
phaseShiftRes = float
phaseShiftInj = float
DELAY_MAX = 100000
eDelayDri = [0.0]*DELAY_MAX
phiDelayDri = [0.0]*DELAY_MAX
eDelayRes = [0.0]*DELAY_MAX
phiDelayRes = [0.0]*DELAY_MAX
eDelayInj = [0.0]*DELAY_MAX
phiDelayInj = [0.0]*DELAY_MAX
delayDriNum = int
delayResNum = int
delayInjNum = int
delayDriIndex = int
delayResIndex = int
delayInjIndex = int
def calcParameter(h):
global tauRes, tauInj, J, kapInj, kapDri, kapRes, tauDri, omega0, phaseShiftDri, phaseShiftRes, phaseShiftInj, delayResNum, delayInjNum, delayDriNum
Nth = N0 + 1.0 / tauP / Gn
Jth = Nth / tauS
J = JRatio * Jth
kapDri = (1 - r2 * r2) * r3Dri / r2 / tauIn
kapRes = (1 - r2 * r2) * r3Res / r2 / tauIn
kapInj = kapInjRatio * kapDri
tauDri = 2.0 * LDri / C
tauRes = 2.0 * LRes / C
tauInj = LInj / C
omega0 = 2.0 * math.pi*C/lambdaa
phaseShiftDri = math.fmod(omega0 * tauDri, 2.0 * math.pi)
phaseShiftRes = math.fmod((omega0 - 2.0 * math.pi * detun) * tauRes, 2.0 * math.pi)
phaseShiftInj = math.fmod(omega0 * tauInj, 2.0 * math.pi)
delayDriNum = int(tauDri / h)
delayResNum = int(tauRes / h)
delayInjNum = int(tauInj / h)
def initializeDelay(a):
global delayDriIndex,delayInjIndex,delayResIndex,DELAY_MAX
delayDriIndex = 0
delayResIndex = 0
delayInjIndex = 0
for item in range(DELAY_MAX):
eDelayDri[item] = a[0]
phiDelayDri[item] = a[1]
for item in range(DELAY_MAX):
eDelayRes[item] = a[3]
phiDelayRes[item] = a[4]
for item in range(DELAY_MAX):
eDelayInj[item] = 0.0
phiDelayInj[item] = 0.0
def laser(x,b,theta):
global kapRes,kapInj,kapDri,delayResIndex,delayInjIndex,delayDriIndex,J
b[0] = 1.0 / 2.0 * (Gn * (x[2] - N0) - 1.0 / tauP) * x[0] + kapDri * eDelayDri[delayDriIndex] * math.cos(theta[0])
b[1] = alpha / 2.0 * (Gn * (x[2] - N0) - 1.0 / tauP) - kapDri * eDelayDri[delayDriIndex]/ x[0] * math.sin(theta[0])
b[2] = J - x[2] / tauS - Gn * (x[2] - N0) * x[0] * x[0]
b[3] = 1.0 / 2.0 * (Gn * (x[5] - N0) - 1.0 / tauP) * x[3] + kapRes * eDelayRes[delayResIndex] * math.cos(theta[1]) + kapInj * eDelayInj[delayInjIndex] * math.cos(theta[2])
b[4] = alpha / 2.0 * (Gn * (x[5] - N0) - 1.0 / tauP) - kapRes * eDelayRes[delayResIndex] / x[3] * math.sin(theta[1]) - kapInj * eDelayInj[delayInjIndex] / x[3] * math.sin(theta[2])
b[5] = J - x[5] / tauS - Gn * (x[5] - N0) * x[3] * x[3]
def rungeKutta(a, h, t):
global delayDriIndex,delayInjIndex,delayResIndex,phaseShiftRes,phaseShiftInj,phaseShiftDri
x=[0.0]*M
b=[[0.0]*M]*4
theta=[0.0]*3
theta[0] = math.fmod(phaseShiftDri + a[1] - phiDelayDri[delayDriIndex], 2.0 * math.pi)
theta[1] = math.fmod(phaseShiftRes + a[4] - phiDelayRes[delayResIndex], 2.0 * math.pi)
theta[2] = math.fmod(phaseShiftInj + a[4] - phiDelayInj[delayInjIndex] - 2.0 * math.pi * detun * t, 2.0 * math.pi)
for item in range(4):
for jtem in range(M):
if item == 0:
x[jtem] = a[jtem]
if item == 1:
x[jtem] = a[jtem] + h * b[0][jtem] / 2.0
if item == 2:
x[jtem] = a[jtem] + h * b[1][jtem] / 2.0
if item == 3:
x[jtem] = a[jtem] + h * b[2][jtem]
laser(x, b[item], theta)
for item in range(M):
a[item] += h * (b[0][item] + 2.0 * b[1][item] + 2.0 * b[2][item] + b[3][item]) / 6.0
eDelayDri[delayDriIndex] = a[0]
eDelayRes[delayResIndex] = a[3]
eDelayInj[delayInjIndex] = a[0]
phiDelayDri[delayDriIndex] = a[1]
phiDelayRes[delayResIndex] = a[4]
phiDelayInj[delayInjIndex] = a[1]
delayDriIndex = (delayDriIndex + 1) % delayDriNum
delayResIndex = (delayResIndex + 1) % delayResNum
delayInjIndex = (delayInjIndex + 1) % delayInjNum
if __name__ == "__main__":
a = [0.0]*M
t = float
h = 5.0e-12
transient = 5000.0e-9
tMax = 50.0e-9
trans = int(transient / h)
n = int(tMax / h)
div = 10
a[0] = 1.3e10
a[1] = 0.0
a[2] = 1.90e24
a[3] = 1.4e10
a[4] = 0.0
a[5] = 1.85e24
initializeDelay(a)
calcParameter(h)
for item in range(trans):
t = h*item
rungeKutta(a, h, t)
pass
timevalue = []
DI = []
RI = []
DDI = []
for item in range(n):
t = h*(trans + item)
if item % div == 0:
print(h*item*1e9, end=' ')
timevalue.append(h*item*1e9)
print(a[0] * a[0] * 1e-20, end=' ')
DI.append(a[0] * a[0] * 1e-20)
print(a[3] * a[3] * 1e-20, end=' ')
RI.append(a[3] * a[3] * 1e-20)
print(eDelayDri[delayDriIndex] * eDelayDri[delayDriIndex] * 1e-20)
DDI.append(eDelayDri[delayDriIndex] * eDelayDri[delayDriIndex] * 1e-20)
rungeKutta(a, h, t)
plt.subplot(3, 1, 1)
plt.plot(timevalue, DI, color='blue')
plt.title("Drive Intensity")
plt.subplot(3, 1, 2)
plt.plot(timevalue, RI, color='red')
plt.title("Response Intensity")
plt.subplot(3,1,3)
plt.plot(timevalue,DDI,color='black')
plt.title("Delayed Drive Intensity")
plt.show()
plt.plot(timevalue, DI, color='blue', label='Drive Intensity')
plt.plot(timevalue, RI, color='red', label='Response Intensity')
plt.plot(timevalue, DDI, color='black', label='Delayed Drive Intensity')
plt.legend(loc='upper left', bbox_to_anchor=(0.01, 0.87))
plt.ylim(1, 18)
plt.show()
开环配置数值模拟python代码:
import math
import matplotlib.pyplot as plt
C=2.99792458e8#光速
M=6#方程的数量
r3Dri=0.008#驱动器外镜反射率
JRatio=1.3#归一化注入电流
LDri=0.6#驱动器外腔长度
LInj=1.2#驱动器和相应激光器之间的距离
detun=0.0e9#光频失谐频率
kapInjRatio=1.0#注入强度比驱动反馈强度
r2=0.556#内镜反射率
tauIn=8.0e-12#内腔往返时间
lambdaa=1.537e-6#光的波长
Gn=8.4e-13#增益系数
N0=1.4e24#透明载流子密度
tauP=1.927e-12#光子寿命
tauS=2.04e-9#载体寿命
alpha=3.0#参数
DELAY_MAX=100000#延迟的最大数组大小
J=float#注入电流
kapDri=float
kapInj=float
tauDri=float
tauInj=float
delayDriNum=int
delayInjNum=int
delayDriIndex=int
delayInjIndex=int
omega0=float
phaseShiftDri=float
phaseShiftInj=float
eDelayDri=[0]*DELAY_MAX
phiDelayDri=[0]*DELAY_MAX
eDelayInj=[0]*DELAY_MAX
phiDelayInj=[0]*DELAY_MAX
def calcParameter(h):
global J,kapInj,phaseShiftInj,phaseShiftDri,delayInjNum,delayDriNum,tauS,JRatio,kapInj,kapDri,tauInj,tauDri,omega0,kapInjRatio,lambdaa,C
Nth = N0 + 1.0 / tauP / Gn#阈值载流子密度
Jth = Nth / tauS#阈值注入电流
J = JRatio * Jth#注入电流
kapDri = (1 - r2 * r2) * r3Dri / r2 / tauIn#驱动反馈强度
kapInj = kapInjRatio * kapDri#从驱动到响应的注射强度
tauDri = 2.0 * LDri / C#驱动的外腔往返时间
tauInj = LInj / C#从驱动到响应的耦合时间
omega0 = 2.0 * math.pi * C /lambdaa#光角频率
phaseShiftDri = math.fmod(omega0 * tauDri, 2.0 * math.pi)#驱动的初始相移
phaseShiftInj = math.fmod(omega0 * tauInj, 2.0 * math.pi)#耦合的初始相移
delayDriNum = int(tauDri/h)
delayInjNum = int(tauInj/h)
def initializeDelay(a):
global delayDriIndex,delayInjIndex,DELAY_MAX
delayDriIndex = delayInjIndex = 0
for item in range(DELAY_MAX):
eDelayDri[item] = a[0]
phiDelayDri[item] = a[1]
for item in range(DELAY_MAX):
eDelayInj[item] = 0.0
phiDelayInj[item] = 0.0
pass
pass
def laser(x,b,theta):
global delayDriIndex,delayInjIndex,Gn,N0,tauP,kapInj,kapDri,alpha,J,tauS
b[0] = 1.0 / 2.0 * (Gn * (x[2] - N0) - 1.0 / tauP) * x[0] + kapDri * eDelayDri[delayDriIndex] * math.cos(theta[0])
b[1] = alpha / 2.0 * (Gn * (x[2] - N0) - 1.0 / tauP) - kapDri * eDelayDri[delayDriIndex] / x[0] * math.sin(theta[0])
b[2] = J - x[2] / tauS - Gn * (x[2] - N0) * x[0] * x[0];
b[3] = 1.0 / 2.0 * (Gn * (x[5] - N0) - 1.0 / tauP) * x[3] + kapInj * eDelayInj[delayInjIndex] * math.cos(theta[1])
b[4] = alpha / 2.0 * (Gn * (x[5] - N0) - 1.0 / tauP) - kapInj * eDelayInj[delayInjIndex] / x[3] * math.sin(theta[1])
b[5] = J - x[5] / tauS - Gn * (x[5] - N0) * x[3] * x[3];
def rungeKutta(a,h,t):
global delayDriIndex,delayInjIndex,M,phaseShiftDri,phaseShiftInj,detun,delayDriNum,delayInjNum
x=[0]*M
b=[[0]*M]*4
theta=[0]*2
theta[0] = math.fmod(phaseShiftDri + a[1] - phiDelayDri[delayDriIndex], 2.0 * math.pi)
theta[1] = math.fmod(phaseShiftInj + a[4] - phiDelayInj[delayInjIndex] - 2.0 * math.pi * detun * t, 2.0 * math.pi)
for item in range(4):
for jtem in range(M):
if item == 0:
x[jtem] = a[jtem]
if item == 1:
x[jtem] = a[jtem] + h * b[0][jtem] / 2.0
if item == 2:
x[jtem] = a[jtem] + h * b[1][jtem] / 2.0
if item == 3:
x[jtem] = a[jtem] + h * b[2][jtem]
laser(x, b[item], theta)
for item in range(M):
a[item] += h * (b[0][item] + 2.0 * b[1][item] + 2.0 * b[2][item] + b[3][item]) / 6.0
#更新延迟数组
eDelayDri[delayDriIndex] = a[0]
eDelayInj[delayInjIndex] = a[0]
phiDelayDri[delayDriIndex] = a[1]
phiDelayInj[delayInjIndex] = a[1]
delayDriIndex = (delayDriIndex + 1) % delayDriNum
delayInjIndex = (delayInjIndex + 1) % delayInjNum
if __name__ == "__main__":
a = [0] * M
h = 5.0e-12#计算步长
transient = 5000.0e-9#瞬态时间
tMax = 50.0e-9#时间步长
trans = int(transient / h)
n = int(tMax / h)
div = 10#绘图间隔
a[0] = 1.3e10#驱动的电场振幅
a[1] = 0.0#用于驱动的电场相位
a[2] = 1.90e24#驱动载流子密度
a[3] = 1.4e10#响应的电场振幅
a[4] = 0.0#响应的电场相位
a[5] = 1.85e24#响应载流子密度
time_value=[]#时间序列
DI=[]#驱动强度
RI=[]#反应强度
DDI=[]#延迟驱动强度
initializeDelay(a)
calcParameter(h)#过渡过程计算的数学模型
for item in range(trans):
t = h * item
rungeKutta(a, h, t)
for item in range(n):
t = h * (trans + item)
if item % div == 0:
print(h * item * 1e9, end=' ')
time_value.append(h * item * 1e9)
print(a[0] * a[0] * 1e-20, end=' ')
DI.append(a[0] * a[0] * 1e-20)
print(a[3] * a[3] * 1e-20, end=' ')
RI.append(a[3] * a[3] * 1e-20)
print(eDelayDri[delayDriIndex] * eDelayDri[delayDriIndex] * 1e-20)
DDI.append(eDelayDri[delayDriIndex] * eDelayDri[delayDriIndex] * 1e-20)
rungeKutta(a, h, t)
plt.subplot(3,1,1)
plt.plot(time_value,DI,color='blue')
plt.title("Drive Intensity")
plt.subplot(3,1,2)
plt.plot(time_value,RI,color='red')
plt.title("Response Intensity")
plt.subplot(3,1,3)
plt.plot(time_value,DDI,color='black')
plt.title("Delayed Drive Intensity")
plt.show()
plt.plot(time_value,DI,color='blue',label='Drive Intensity')
plt.plot(time_value,RI,color='red',label='Response Intensity')
plt.plot(time_value,DDI,color='black',label='Delayed Drive Intensity')
plt.legend(loc='upper left', bbox_to_anchor=(0.01, 0.87))
plt.ylim(1,19)
plt.show()
运行的代码,依照光通信大书中附件的c语言代码改写成的python代码,运用的积分思想是数学中的四阶龙格库塔的方法,对于上次的模型的代码修改还在进行中。欢迎佬们对上次二阶龙格库塔解混沌方程的方法进行指导!