24小时热门版块排行榜    

查看: 433  |  回复: 1
当前主题已经存档。

freebjx

金虫 (小有名气)


[资源] 【讨论】一个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);

回复此楼

» 猜你喜欢

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

freebjx

金虫 (小有名气)


结果应该是双峰,高手请看看有什么问题
2楼2009-04-21 09:42:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 freebjx 的主题更新
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-08-31 11:51 by xiangy672
[基金申请] 29号明天会评吗 +4 笨笨唐 2026-08-28 4/200 2026-08-31 09:30 by huixian257
[基金申请] 国社科又开始会评了,不知道这次命运如何 +6 雨打竹帘 2026-08-30 10/500 2026-08-31 09:29 by huixian257
[文学芳草园] 梦想 +4 myrtle 2026-08-26 4/200 2026-08-31 08:43 by 逸域隐帆
[基金申请] 面上意见出来了 +10 黄鸟于飞Chao 2026-08-29 19/950 2026-08-31 08:01 by blueearth171
[基金申请] 为什么到现在没收到通知? +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[考博] 找导师 +6 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 国合现在查不到了吗? +10 chengyan1220 2026-08-24 20/1000 2026-08-26 08:57 by peasantsprig
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
信息提示
请填处理意见