fsolve解两个耦合的二元三次非线性方程组，程序可运行但结果...
%%%%%%%%%%%%%%%% 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=; % 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');
可能是初值的问题，研究了一个多月了无法解决，希望高手指点一下。
系数delta改变时，方程的两个解也在改变，这时初值怎么选取？ c1的结果应该是大于或等于0的实数，c2的结果应该在-0.5至0之间，才是合理的。 我试了一下运行不了 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=; % 初值
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'); eastar 发表于 2017-2-22 11:03
我试了一下运行不了
我的机上可以运行，注意M文件中有换行的 yurong 发表于 2017-2-22 11:40
我的机上可以运行，注意M文件中有换行的
我再试试 eastar 发表于 2017-2-22 13:18
我再试试
多谢你了。 yurong 发表于 2017-2-22 13:23
多谢你了。
指点我一下。 学习了
页:
[1]