| 查看: 431 | 回复: 1 | ||
| 当前主题已经存档。 | ||
[资源]
【讨论】一个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);  |
回复此楼» 猜你喜欢
有多少人是今天查系统知道结果的？
已经有19人回复
-
国社科又开始会评了，不知道这次命运如何
已经有9人回复
-
能否申诉？
已经有5人回复
-
为什么到现在没收到通知？
已经有5人回复
-
有没有仍没收到信息的
已经有8人回复
-
面上意见出来了
已经有18人回复
-
找导师
已经有7人回复
-
29号明天会评吗
已经有3人回复
-
麻烦专家们看看评委们的意见（F口面上）
已经有11人回复
-
我就是申请一个面上项目而已，这评审意见是按照杰青的条件评的吧？
已经有11人回复
2楼2009-04-21 09:42:53
投票: