fsolve解两个耦合的二元三次非线性方程组，程序可运行但结果...

yurong

%%%%%%%%%%%%%%%% 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');

可能是初值的问题，研究了一个多月了无法解决，希望高手指点一下。

yurong

系数delta改变时，方程的两个解也在改变，这时初值怎么选取？

yurong

c1的结果应该是大于或等于0的实数，c2的结果应该在-0.5至0之间，才是合理的。

eastar

我试了一下运行不了

yurong

eastar 发表于 2017-2-22 11:03
我试了一下运行不了

M-file文件有点错误修改了一下
%%%%%%%%%%%%%% 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)*conj(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)*conj(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];                               % 初值

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');

yurong

eastar 发表于 2017-2-22 11:03
我试了一下运行不了

我的机上可以运行，注意M文件中有换行的

eastar

yurong 发表于 2017-2-22 11:40
我的机上可以运行，注意M文件中有换行的

我再试试

yurong

eastar 发表于 2017-2-22 13:18
我再试试

多谢你了。

yurong

yurong 发表于 2017-2-22 13:23
多谢你了。

指点我一下。

dh492510085

学习了