24小时热门版块排行榜

返回列表

当前主题已经存档。

freebjx

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1284
帖子: 228
在线: 9.5小时
虫号: 612208

[资源] 【讨论】一个matlab解复变量方程实例

忙了较长一段时间了，终于有点结果，首先得感谢fspdlh，大部分工作是他完成的。下面是程序和结果供大家讨论，还望大家多指点，结果的数据有点发散，理论上应是光滑曲线，请高手指点。其中S.a,S.b分别有三个数据，我取了 S.a（1）,S.b（1）。其它的也试过
function equation()
clear
clc
%--------------------------------------------------------------------------
kappa = 0.44;
beta=9.6;
kt=1.2;
gamma1=0.16;
gamma2 =2.4;
g0 = 0.6;
wc =0;
wa=0;
E=0.296e4;
s_j=E/(i*(sqrt(2*kappa)));
ns=gamma1*gamma2/(4*g0^2);
nloop =4001;
wlist = linspace(-30,30,nloop);
pt1 = zeros(nloop,1);
pr1=zeros(nloop,1);
c0=E/(sqrt(2*ns)*kt);
c1=4*g0^2/(2*kt*gamma2);
%--------------------------------------------------------------------------
fun1='2*a^3+2*a*b^2+a*c2+a*c1-2*b*c3*a^2-2*c3*b^3-b*c3*c2+b*c4-2*c0*a^2-2*c0*b^2-c0*c2=0';
fun2='2*b*a^2+2*b^3+b*c2+b*c1+2*c3*a^3+2*a*c3*b^2+a*c3*c2-a*c4=0';
fun1=subfun(fun1,'c0',c0);
fun1=subfun(fun1,'c1',c1);
fun2=subfun(fun2,'c0',c0);
fun2=subfun(fun2,'c1',c1);
for k = 1:nloop
wl=wlist(k);
c2=(wa-wl)^2/kt^2+1;
c3=((wc-wl)-beta)/kt ;
c4=c1*(wa-wl)/gamma2;
temp_fun1=subfun(fun1,'c2',c2);
temp_fun1=subfun(temp_fun1,'c3',c3);
temp_fun1=subfun(temp_fun1,'c4',c4);
temp_fun2=subfun(fun2,'c2',c2);
temp_fun2=subfun(temp_fun2,'c3',c3);
temp_fun2=subfun(temp_fun2,'c4',c4);
S=solve(temp_fun1,temp_fun2,'a','b');

a1=S.a(1);
b1=S.b(1);
x1=real(a1)+i*(b1);
asw1=sqrt(ns)*x1;
x2=c0/(1+i*((wc-x1)+beta)/kt );
asw2=sqrt(ns)*x2;
ac=sqrt(2)*(asw1+asw2)/2;
acc=sqrt(2)*(asw1-asw2)/2;
%----------------------------------------------------------------------
t=s_j+(i*sqrt(2*kappa))*ac;
pt(k)=t'*t;
r=(i*sqrt(2*kappa))*acc;
pr(k)=r'*r;
end

%--------------------------------------------------------------------------
pt=eval(pt);
pr=eval(pr);
plot(wlist,real(pt),wlist,real(pr));
xlabel('dwl'); ylabel('cavity transmission/reflection');
legend('transmission','reflection');

%==========================================================================
function fun=subfun(fun,str,para)
para=cat(2,'(',num2str(para),')');
fun=strrep(fun,str,para);