|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
%%%%%%%%%%%%%%%% M-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dx=SF_g1(x)
global delta;
global Omega;
global Sp;
global g;
global J;
global kat;
global kbt;
global kae;
global gamma;
global gamma_d;
global P;
global chi;
delt_a=1i*delta+0.5*kat;
delt_b=1i*2*delta+0.5*kbt;
delt_q=1i*(2*delta+Omega)+0.5*P+0.5*gamma+gamma_d;
dx(1)=4*chi^2*abs(delt_q)^2*g^4*x(1)^3 ...
+(2*chi*g^2*((delt_a*delt_b+conj(delt_a)*conj(delt_b))*abs(delt_q)^2-2*(delt_a*conj(delt_q)+conj(delt_a)*delt_q)*J^2*x(2)))*x(1)^2 ...
+abs(delt_a)^2*(abs(delt_b)^2*abs(delt_q)^2-2*(delt_b*delt_q+conj(delt_b)*conj(delt_q))*J^2*x(2)+4*J^4*x(2)^2)*x(1) ...
-kae*abs(delt_b*delt_q-2*J^2*x(2))^2*Sp^2;
dx(2)=-4*(gamma+P)*J^4*x(2)^3 ...
+2*J^2*((gamma+P)*(delt_b*delt_q+conj(delt_b*delt_q)+(P-gamma)*J^2))*x(2)^2 ...
-((gamma+P)*abs(delt_b)^2*abs(delt_q)^2+2*(delt_q+conj(delt_q))*J^2*g^2*x(1)^2+(P-gamma)*(delt_b*delt_q+conj(delt_b*delt_q))*J^2)*x(2) ...
+0.5*(P-gamma)*abs(delt_b)^2*abs(delt_q)^2;
%%%%%%%%%%%%%%%%%%%%% running-file %%%%%%%%%%%%%%%%%%%%%%%
clear;
global delta;
global Omega;
global Sp;
global g;
global J;
global kat;
global kbt;
global kae;
global gamma;
global gamma_d;
global P;
global chi;
chi=1;
Omega=0;
c0=2.997925e8;
namda_p=1550e-9;
omega_p=2*pi*c0/namda_p;
hbar=6.6260693e-34/(2*pi);
Pp=50e-6;
Sp=sqrt(Pp/(hbar*omega_p)); % \varepsilon_p
g=2*pi*2*10^9;
J=2*pi*10*10^9;
P=0;
Qai=7*10^5;
Qat=1.3*10^4;
lambda=1550*10^(-9);
kai=2*pi*c0/(lambda*Qai); % kai=pi*c/(lambda*Qi)
kat=2*pi*c0/(lambda*Qat); % kae=pi*c/(lambda*Qe)
kae=kat-kai;
Qbi=7*10^5;
Qbt=1.3*10^4;
lambdb=775*10^(-9);
kbi=2*pi*c0/(lambdb*Qbi); % kbi=pi*c/(lambda*Qi)
kbt=2*pi*c0/(lambdb*Qbt); % kbe=pi*c/(lambda*Qe)
kbe=kbt-kbi;
gamma=2*pi*0.16*10^9;
gamma_d=2*pi*1*10^9;
x0=[500;-0.5]; % x0(i)对应与xi的初值
f00=-50:0.1:50;
for j=1:length(f00)
delta=f00(j)*2*pi*1*10^9;
delt_a=1i*delta+0.5*kat;
delt_b=1i*2*delta+0.5*kbt;
delt_q=1i*(2*delta+Omega)+0.5*P+0.5*gamma+gamma_d;
y=fsolve(@SF_g1,x0,optimset('Display','off'));
c1(j)=y(1); % c1>=0
c2(j)=y(2); % -0.5=<c2<=0
end
plot(f00,c1,'-b',f00,c2,'-r');
可能是初值的问题,研究了一个多月了无法解决,希望高手指点一下。
|
|