使用ode在MATLAB中求解环形激光方程
我写了一个Matlab代码，它是使用环形激光器模拟非门实现的。 Si变量是在我的代码中输入的。 但是我的代码只适用于Si = 0 。 对于任何其他非零值，它不显示输出。
%-----------Rate Equation MATLAB code ---------
function dydt = requations(t,y)
dydt = zeros(size(y));
%prompt = 'Sinput???';
%Si = input(prompt);
Si = 0; %MY INPUT
q = 1.6e-19; % charge of electron
tau_e = 1e-9; % carrier lifetime
No = 3.3e18; % # No of carriers at transparency
a = 1e-15; % Linear gain coefficient
Vg = 7.5e9; % group velocity
Vp = 3e-11; %Photon reservoir volume
V = 1e-11; %Carrier reservoir Volume
tau_p = 1.7e-12; % photon lifetime
beta = 1e-5; % spontateous emission coefficient
eps = 7.5e-17; % Nonlinear gain suppression coefficient
Ni = 0.8; %Internal quantum efficiency
w = 2*pi*10e5;
Io = 10e-3;
%I = Io*sin(w*t);
I = 2.5*Io; %for test purposes
tp = 1/tau_p;
te = 1/tau_e;
Aint = 6; %Internal losses inside cavity waveguides
%G = ((a/Vp)* (N-(V*No)))/(1-eps*(Sc + Scc));
alpha = -2; %alpha factor
L = 76e-4 ;%size of the ring
wcb = 2*pi*70;
%R = 0.25;
Wcb = wcb*1000000;
%r = 1/R;
tpcw = Vg*(Aint + ((1/L)*log(4)));
Tpcw = 1/tpcw;
%------Rate equations-------
N = y(1);
Sc = y(2); %Clock wise photon number
yc = y(3);
Scc = y(4); %anti clockwise photon number
ycc = y(5);
G = ((a/Vp)* (N-(V*No)))/(1-eps*(Sc + Scc));
dydt(1) = (Ni*I)/q - y(1)*te - Vg*G*(y(2) + y(4)); %dN/dt
dydt(2) = (Vg*G-Tpcw)*y(2) + beta*y(1)*te; %dSc/dt
dydt(3) = -(alpha/2)*(Vg*G-Tpcw); %dyc/dt
dydt(4) = (Vg*G-Tpcw)*y(4) + beta*y(1)*te + ((2*Vg)/L)*cos(y(5))*(sqrt(Si*y(4))); %dScc/dt
dydt(5) = -Wcb - ((alpha/2)*(Vg*G-Tpcw)) - ((Vg/L)*sin(y(5))*(sqrt(Si/y(4)))); %dycc/dt
以下是Ode文件
%------Rate equations for requation file------
format bank;
close all;
clear all;
clc;
%time interval
ti=0;
tf=200;
tspan=[ti tf];
x0 = [3.75e7, 2.25e6, 0, 2.25e6, 0]; %initial vectors
%options= odeset('RelTol',100, 'AbsTol',[3.75e7, 2.25e6]);
[t,y]= ode23t(@requations,tspan,x0);
%Plotting the graphs:
figure
subplot(5,1,1), plot(t,y(:,1),'r'),grid on;
title('Laserrate equations'),ylabel('N');
subplot(5,1,2), plot(t,y(:,2),'b'),grid on;
ylabel('Scw'); xlabel('t');
subplot(5,1,3), plot(t,y(:,3),'g'),grid on;
ylabel('ycw');xlabel('t');
subplot(5,1,4), plot(t,y(:,3),'g'),grid on;
ylabel('Sccw');xlabel('t');
subplot(5,1,5), plot(t,y(:,3),'g'),grid on;
ylabel('yccw');xlabel('t');