matlab求解速率方程,关于用Matlab求解3能级系统速率方程的困惑

function dy = rate_eq( t,y,flag )

%红宝石激光器参数如下

W13=1e3;

S32=0.5e7;

A31=3e5;

sigma21=2.5e-20;

V=1.7e10;

A21=0.3e2;

S21=0;%很小就认为0吧

NT=1.6e19;

taoL=1e-18;

%光子数密度初值

%phi--y(1);N2--y(2);N3=y(3)

dy=zeros(3,1);

y(1)=max(y(1),100);

% dy=[(y(2)-(NT-y(2)-y(3)))*sigma21*V*y(1)-y(1)/taoL;

%     -(y(2)-(NT-y(2)-y(3)))*sigma21*V*y(1)-y(2)*(A21+S21)+y(3)*S32;

%     (NT-y(2)-y(3))*W13-y(3)*(S32+A31)];

dy=[(2*y(2)-NT+y(3))*sigma21*V*y(1)-y(1)/taoL;

(NT-2*y(2)-y(3))*sigma21*V*y(1)-y(2)*(A21+S21)+y(3)*S32;

(NT-y(2)-y(3))*W13-y(3)*(S32+A31)];

end

clear

clc

%phi--y(1);N2--y(2);N3=y(3)

y0=[100;0;0];%设定初值

%tspan=[0 1.05];

tspan=[0 1e-3];

[t,y]=ode45('rate_eq',tspan,y0);

NT=1.6e19;

N1=NT-y(:,2)-y(:,3);

subplot(3,2,1)

plot(t,N1)

xlabel('t/ms')

ylabel('N1')

subplot(3,2,3)

plot(t,y(:,2))

xlabel('t/ms')

ylabel('N2')

subplot(3,2,5)

plot(t,y(:,3))

xlabel('t/ms')

ylabel('N3')

subplot(3,2,2)

plot(t,2*y(:,2)-NT+y(:,3))

xlabel('t/ms')

ylabel('N2-N1')

subplot(3,2,6)

plot(t,y(:,1))

xlabel('t/ms')

ylabel('\phi')

%end

A.jpg

(56.46 KB, 下载次数: 21)

2013-5-21 21:48 上传

速率方程