声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2190|回复: 9

[编程技巧] fsolve解两个耦合的二元三次非线性方程组,程序可运行但结果...

[复制链接]
发表于 2017-2-22 09:25 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

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

回复
分享到:

使用道具 举报

 楼主| 发表于 2017-2-22 09:29 | 显示全部楼层
系数delta改变时,方程的两个解也在改变,这时初值怎么选取?
 楼主| 发表于 2017-2-22 10:58 | 显示全部楼层
c1的结果应该是大于或等于0的实数,c2的结果应该在-0.5至0之间,才是合理的。
发表于 2017-2-22 11:03 | 显示全部楼层
我试了一下运行不了
 楼主| 发表于 2017-2-22 11:37 | 显示全部楼层
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');
 楼主| 发表于 2017-2-22 11:40 | 显示全部楼层
eastar 发表于 2017-2-22 11:03
我试了一下运行不了

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

点评

我再试试  详情 回复 发表于 2017-2-22 13:18
发表于 2017-2-22 13:18 | 显示全部楼层
yurong 发表于 2017-2-22 11:40
我的机上可以运行,注意M文件中有换行的

我再试试
 楼主| 发表于 2017-2-22 13:23 | 显示全部楼层

多谢你了。
 楼主| 发表于 2017-2-22 15:25 | 显示全部楼层

指点我一下。
发表于 2017-6-11 18:55 | 显示全部楼层
学习了
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-6 03:45 , Processed in 0.074105 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表