蔡氏电路仿真实验

模拟电路运行

采用MATLAB进行模拟。本电路为一常微分方程的初值问题。

取定参数

matlab代码实现

首次模拟的图像对比

参数调整后，模拟的图像对比

分析频谱

模拟非线性电路同步

非线性电路同步瞬态过程

取定参数

模拟实验中的电感电容都采用测量值，为了方便，取为定值。

非线性负阻的I-V特性取成奇对称，数值和测量值稍作对称性修正。

代码实现

这里使用Matlab中ode45，四阶龙格库塔法四阶龙格库塔法求解常微分方程初值问题。C和C++的代码

I-V函数 myfun.m

function fun=myfun(x)

if x<0

fun=-myfun(-x);

elseif x<1.664849891

fun=x*(-7.44316E-04)-5.065E-06;

elseif x<11.29059205

fun=x*(-4.003E-04)-5.778E-04;

else

fun=x*3.574E-03-4.545E-02;

end

微分方程组 myode.m

function dy=myode(t,y)

global gg;

dy = zeros(3,1);

G=1/gg;

C1=9.91E-9;

C2=98.2E-9;

L=23E-3;

dy=[(G*(y(2)-y(1))-myfun(y(1)))/C1;

(G*(y(1)-y(2))+y(3))/C2;

-y(2)/L];

end

调节1/G输出图像nonlinear.m

clear;

global gg;

for gg=1800:1:2200

fs=100000;

[T,Y]=ode45('myode',0:1/fs:0.2,[0;0;0]);

X=Y(10002:end,1);

plot(Y(10002:end,1),Y(10002:end,2));

grid on

axis([-13 13 -3 3]);

xlabel('U1/V');

ylabel('U2/V');

str=[num2str(gg) '.jpg'];

saveas(gcf,str);

end

首次模拟的结果

改变

的由小到大数值，作出

相图，可以观察到一下几种典型图样

之中的几个稳态

图像在屏幕上收缩成一点

图形变化的规律和真实实验观察到的基本一致，但是图像对应的1/G和形状与真实实验不太一致，究其原因，是我的L，C1，C2测量不准确。

参数调整后，模拟的结果

将L、C1、C2换为参考值，再次计算。1/G在1530到2040之间变化。

C1=10E-9;

C2=100E-9;

L=18E-3;

改变参数以后图形和实验更相似,但是变化的规律不变

分岔图分析

分岔图可以直观得看到分岔，这里用U1的最小值来画分岔图，代码如下：

clear;

fs=100000;

for gg=1992:-0.1:1940;

[T,Y]=ode23(@myode,0:1/fs:0.2,[0;0;0]);

data=Y(:,1);

n=length(data);

N=n-round(n/4);

for i=N:n-2

if data(i)

plot(gg,data(i));

hold on;

else

end

end

end

beep

分岔图：

G从1940到1990G从1960到1990(放大)

G=1991为单倍周期G=1980为2倍周期G=1971.6为4倍周期

介于matlab提供的龙格库塔法函数精度很难控制，且计算缓慢，建议使用C语言进行数值计算，并提高精度。

从左边大范围的分岔图上已经清晰显示了1→2→4→8→6→5→3的分岔过程。我们将分岔图区域放大后(如右图)，我们发现分岔图具有自相似性，即分形的最大特点。

处理与分析

利用快速傅立叶变换分析频谱

X=fft(X);

N=length(X);

X(1) = [];

power = abs(X(1:N/2)).^2;

[mp,index] = max(power);

nyquist = 1/2;

freq = (1:N/2)/(N/2)*nyquist*fs;

figure(1);

subplot(211);

plot(freq,power), grid on

title (freq(index))

xlabel('Freqency');

ylabel('Power');

period = 1./freq;

figure(1);

subplot(212);

plot(period,power), grid on

title (period(index))

ylabel('Power');

xlabel('Period');

G=1500、1827、1946、1999采样

G相图功率谱最大周期20倍时长时域图说明

1500

1827功率谱有5个峰，而吸引子有5层圈。最大峰位的频率与其它相比显著小

1946

1971.7

1980

1999

模拟非线性电路同步

随机产生两个混沌电路的初始参数，只需要改变如下一行。

[T,M]=ode45(@myode,0:1/fs:0.19,[rand();rand();rand()/1000]);

未同步，两个电路的U1未同步，两个电路的U2

同步U1，得到两个U2的相图同步U1，C2变大5%，得到两个U2的相图

非线性电路同步的瞬态过程

在极短时间内两个电路达到了同步

同步过程的相图同步过程的时域图

2016 nonlinear circuit with python implementation

Claim: This is a brief documentation of using python to understand the simulation of the Nonlinear circuit problem.

Background: Compared with Matlab , python is weigh more powerful in handling object-oriented simulation despite of simulink.

Also,python is really a hot languange which will be beneficial for our scientific computing.

这里面要求环境装上scipy和numpy,一般windows用pip，而OS用户就可以用brew，linux的话大家都懂

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt

def Nlfunction(V_1):

Ga=-0.00076

Gb=-0.00049

E=15.0

retg=Gb*V_1+(Gb-Ga)/2*(abs(V_1-E)-abs(V_1+E))

return retg

def delta(y,t):

global R

C1=9.91e-9

C2=98.2e-9

l=23e-3

G=1/R

ret_delta =np.array([(G*(y[1]-y[0])-Nlfunction(y[0]))/C1,

(G*(y[0]-y[1])+y[2])/C2,

-y[1]/l])

return ret_delta

time=np.linspace(0,0.2,10000)

r=range(1930,1960,1)

yinit=np.array([0.0,0.0,0.1])

global R

for R in r:

y = odeint(delta,yinit,time)

plt.figure()

plt.plot(y[1000:,0],y[1000:,1])

plt.xlabel('U_1')

plt.ylabel('U_2')

plt.legend()

plt.savefig('%domega'%R)

plt.clf()

张贴一些相图

R=1800

R=1841

R=1859

to be continued —Yufeng Ma

讨论区

罗页，向你致敬！

接下来我苦了，要学习的内容太多！

— 乐永康 2008/12/15 23:14就感觉东西好多好杂啊!

— 罗页 2008/12/16 00:05

很赞马同学的尝试。 — 乐永康 2016/05/17 22:30之前用了ubuntu所以没有用中文写，终于可以切换到windows了，说正题，老师我想问一下这个非线性元件里的E是怎么计算的，我感觉跟学长模拟出来的区间不一样

我怀疑是这个E的值不太对，另外我查文献貌似15V输入到里面会有饱和电压，不知道是什么问题。

乐老师，我不知道怎么把图贴上来啊，为什么点插入图片毫无用处？ — 马雨枫 2016/5/18 21:13编辑栏的倒数第五个按钮是上传附件。图片上传到服务器了，就可以插入到页面了。有可能你没有足够的权限。

电压15V时饱和，那是因为那个负阻的电源电压是15V。余下问题来面谈吧。 — 乐永康 2016/05/19 00:16

非线性负阻的I-V曲线必须包含Gb外侧的另外两个区域，不然模拟的过程中，双吸引子出不来。另外，odeint使用的不是RK4，而是‘Isoda’。应该使用scipy.integrate.RK45这个类才是4阶RK，上面的matlab代码使用的ode23可能只算到2阶，所以精度上不去也是正常的。 — 戴植锐 2019/04/05 20:31