Lang–Kobayashi方程实现混沌python实现混沌序列图像

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代码,运用的积分思想是数学中的四阶龙格库塔的方法,对于上次的模型的代码修改还在进行中。欢迎佬们对上次二阶龙格库塔解混沌方程的方法进行指导!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mfbz.cn/a/263616.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

beaglebone black狗板,交叉编译Qt5(eglfs)

1. 下载buildroot-2023.023.7版本 make beaglebone_qt5_defconfig 然后编译,出现错误大多数是因为下载不了包,用bing搜索找到放到对应的dl目录下,最终完成编译。 备注:用系统默认配置,不要参考网上的,网…

天津web前端就业培训班,Web机构选择重点

Web前端培训是目前非常热门的培训领域之一。很多领域都会涉及到web前端开发,比如传统互联网、房地产、金融、游戏、影视传媒等行业都需要web前端技术的支持。越来越多的企业和个人也需要建立自己的网站和移动应用程序,因此市场对web前端工程师的需求是非…

Docker 学习总结(80)—— 轻松驾驭容器,玩转 LazyDocker

前言 LazyDocker 是一个用户友好的命令行工具,简化了 Docker 的管理。它能够通过单一命令执行常见的 Docker 任务,如启动、停止、重启和移除容器。LazyDocker 还能轻松查看日志、清理未使用的容器和镜像,并自定义指标。 简绍 LazyDocker 是一个用户友好的 CLI 工具,可以轻…

Linux基本内容学习

Linux 命令 文件命令 命令释义语法格式lslist,用于显示目录中文件及其属性信息ls [参数名] [文件名]cdchange directory,用于更改当前所处的工作目录,路径可以是绝对路径,也可以是相对路径,若省略不写则会跳转至当前…

使用GitZip下载GitHub指定文件

目录 一、GitZip二、安装GitZip三、链接GitHub四、检验是否安装成功五、总结 一、GitZip GitZip是一个非常实用的浏览器插件,它主要有以下几个优点: 下载指定文件:在我们浏览Github时,如果只想下载某个子目录的内容,…

【爬虫软件】孔夫子二手书采集

项目演示 孔网爬取图书信息 目录结构 [ |-- api-ms-win-core-synch-l1-2-0.dll, |-- api-ms-win-core-sysinfo-l1-1-0.dll, |-- api-ms-win-core-timezone-l1-1-0.dll, |-- api-ms-win-core-util-l1-1-0.dll, |-- api-ms-win-crt-conio-l1-1-0.dll, |-- api…

STL中优先队列(堆)的详解

文章目录 priority_queue的基本介绍堆(heap)堆的概念与结构 priority_queue 的介绍与使用 priority_queue的基本介绍 这个priority_queue翻译成中文就是优先级队列,但其实我们很难去一眼看出他的意思到底是什么,他的逻辑结构实际上类似于数据结构中的堆…

王美莲(瑞美)在博鳌乡村振兴产业发展大会荣获“乡村振兴带头人”殊荣

近日,博鳌乡村振兴产业发展大会2023年会,在博鳌亚洲论坛国际会议中心隆重召开。本届大会由中国乡村发展协会指导,中国民族贸易促进会、博鳌乡村振兴产业发展大会组委会主办,深圳市益米播科技有限公司承办。 右:原文化部副部长、国家博物馆首任馆长潘震宙 左:广州市越秀区大塘…

电容内容介绍

0 Preface/Foreword 电容,Capacitance,i.e. 电容量,指在给定电位差下自由电荷的储存量,符号为C,单位为F(法拉)。 电容,指容纳电荷的能力。任何静电场都是由许多电容组成&#xff0…

HTML---盒子模型

文章目录 前言一、pandas是什么?二、使用步骤 1.引入库2.读入数据总结 一.盒子模型概述 HTML中的盒子模型是一种用于描述和布局元素的概念。每个 HTML 元素都可以被表示为一个矩形的盒子,这个盒子包括四个部分:内容区域、内边距、边框和外边距…

图灵日记之java奇妙历险记--数据类型与变量运算符

目录 数据类型与变量字面常量数据类型变量语法格式整型变量浮点型变量字符型变量希尔型变量类型转换自动类型转换(隐式)强制类型转换(显式) 类型提升不同数据类型的运算小于4字节数据类型的运算 字符串类型 运算符算术运算符关系运算符逻辑运算符逻辑与&&逻辑或||逻辑非…

SpringIOC之MethodBasedEvaluationContext

博主介绍:✌全网粉丝5W,全栈开发工程师,从事多年软件开发,在大厂呆过。持有软件中级、六级等证书。可提供微服务项目搭建与毕业项目实战,博主也曾写过优秀论文,查重率极低,在这方面有丰富的经验…

MyBatis的增删改查(CRUD)

目录 查询 1.单个参数绑定 2.序号参数绑定 3.注解参数绑定(推荐) 4.对象参数绑定(推荐) 5.Map参数绑定 6.模糊查询 7.sql注入 8.聚合函数查询 删除 修改 添加(主键回填) 数据库 程序结构 查询 …

GCC:GNU编译器

GCC(GNU Compiler Collection)是一款广泛使用的开源编译器套件,支持多种编程语言,包括C、C、Objective-C、Fortran、Ada和Go等。在本文中,我们将通过一个简单的C程序来介绍GCC的编译过程,包括预处理、编译、…

ECMAScript基础入门:猫头虎博主的技术分享

🌷🍁 博主猫头虎 带您 Go to New World.✨🍁 🦄 博客首页——猫头虎的博客🎐 🐳《面试题大全专栏》 文章图文并茂🦕生动形象🦖简单易学!欢迎大家来踩踩~🌺 &a…

怎么定义一套完成标准的JAVA枚举类型

一、背景 在java代码中,接口返回有各种各样的状态,比如400 401 200 500 403等常见的http状态码,也有我们自定义的很多业务状态码。如果系统比较复杂,制定一套完整的标准的状态码是非常有必要的,这样比较方面BUG排查。…

六个探索性数据分析(EDA)工具,太实用了!

当进行数据分析时,探索性数据分析(EDA)是一个至关重要的阶段,它能帮助我们从数据中发现模式、趋势和异常现象。而选择合适的EDA工具又能够极大地提高工作效率和分析深度。在本文中,笔者将介绍6个极其实用的探索性数据分析(EDA)工具&#xff0…

【Python小知识 - 6】:QLabel设置图片

文章目录 QLabel设置图片 QLabel设置图片 from PyQt5.QtWidgets import * from PyQt5.QtGui import * import sysapp QApplication(sys.argv)window QWidget()hbox QHBoxLayout(window)# 设置标签图片 lable QLabel() lable.setPixmap(QPixmap(./img/window.png).scaled(1…

Python中实现消息队列:构建高效异步通信系统的完整指南

更多资料获取 📚 个人网站:ipengtao.com 消息队列的基础概念 在开始之前,先了解一些消息队列的基础概念。 1 什么是消息队列? 消息队列是一种通信方式,它允许将消息从一个应用程序传递到另一个应用程序。消息队列提…

2022xctf-final hole

这个题是做到的第一个利用hole和map来制造oob的题目,挺有意思的记录一下 首先根据题目给出的信息可知涉及到此漏洞 https://crbug.com/1263462 poc如下: let theHole %TheHole(); m new Map(); m.set(1, 1); m.set(theHole, 1); m.delete(theHole);…