Problem solving laser rate equations with ode45 command

6 views (last 30 days)
Hi there, I have a question about modelling rate equations. Basically the modelling is solving differential equations hence I'm trying to use ODE45.The equations are as follows:
It should also be noted that the two values of and are as follows
I have written the above equations in a program as follows.
------Main program-----------
clear all
tspan = [0 2d-9]; % time interval, up to 2 ns
y0 = [0,0,0,0];
[t,y] = ode45('rate_eq_program_1',tspan,y0);
size(t);
t=t*1d9;
y_max = max(y);
y1 = y_max(1);
y2 = y_max(2);
y3 = y_max(3);
y4 = y_max(4);
d = plot(t, y(:,1)/y1,'-.', t, y(:,2)/y2,'--', t, y(:,3)/y3,'-*', t, y(:,4)/y4,'-.'); % divided to normalize
xlabel('time [ns]','FontSize',14); % size of x label
ylabel('Arbitrary units','FontSize',14); % size of y label
set(gca,'FontSize',14); % size of tick marks on both axis
legend('carrier density','nanocavity carrier density','forward field','backward field') % legend inside the plot
pause
close all
-------function--------
function y = rate_eq_program_1(t,y)
%
param_rate_eq_fano % input of needed parameters
%
current1 = 4d-2; % bias current (step function) [A]; 400mA
sigmaa=(2.*epsilon0.*ref_index.*ref_indexg)./(hbar.*w_s).*((1+(abs(r_R)).^2).*(1-r_R))./(conf.*V_g.*g_n.*(y(1)-N_0));
ydot(1) = current1./(e.*V_a)-y(1)./tau1-(V_g.*g_n.*(y(1)-N_0).*sigmaa.*(abs(y(3))).^2)./V_m;
ydot(2) = -y(2)./tau1-(conf_NC.*V_g.*g_n.*(y(2)-N_0).*rho.*(abs(y(4))).^2)./V_NC;
ydot(3) = 1/2.*(1-1i.*henry).*conf.*V_g.*g_n.*(y(1)-N_ss).*y(3)+gamma_L.*((y(4)./r_R-y(3)));
ydot(4) = (-1i.*delta_w-gamma_T).*y(4)-p.*gamma_c.*y(3)+1/2.*(1-1i.*henry).*conf_NC.*V_g.*g_n.*(y(2)-N_0).*y(4);
ydot = ydot'; % must return column vector
-------param_rate_eq_laser--------
c = 3d10; % velocity of light [cm/s]
e = 1.6021892d-19; % elementary charge [C]
h_Planck = 6.626176d-34; % Planck constant [J s]
hbar = h_Planck/(2.0*pi); % Dirac constant [J s]
% Geometrical dimensions
L = 5d-4;
w = 9d-5;
d = 80d-8;
%V = L*w*d;
V_a = 5.26d-7;
conf = 0.5;
conf_NC =0.3;
V_m = V_a/conf;
V_NC=2.4d-7;
ref_index = 3.5;
ref_indexg=3.5;
V_g = c/ref_index;
tau1 = 5d-10;
g_n = 5d-16;
N_0=1d18;
N_ss=3d18;
henry_i=100000;
henry=1;
gamma_L=8.5*10^12;
gamma_c=383*10^9;
gamma_T=387*10^9;
w_r=193*10^12;
w_s=196*10^12;
w_c=250*10^12;
p=-1;
epsilon0=8.85*10^-14;
rho=(2*epsilon0*ref_index*c)./gamma_c*hbar*w_r;
r_L=1;
delta_w=(w_c-w_s);
And I get the following answer at the output.
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
Even if the given parameters have errors for any reason, this should not cause zero output.
Thanks in advance for helping!!

Accepted Answer

Star Strider
Star Strider on 3 Nov 2019
Edited: Star Strider on 3 Nov 2019
Change the function declaration line of your differential equation system function to:
function ydot = rate_eq_fano_program_1(t,y)
That should work.
The problem is that it is returning ‘y’ (that is the input to the function), and is not returning ‘ydot’ that it should be returning.
Also, the problem with the plot call is that:
y_max =
232.9712e+012 0.0000e+000 0.0000e+000 0.0000e+000
so only ‘y(:,1)’ plots. The rest (all normalised by ‘ymax’) are Inf, and will not plot.
  8 Comments

Sign in to comment.

More Answers (1)

Steven Lord
Steven Lord on 3 Nov 2019
Your initial condition vector is:
y0 = [0,0,0,0];
What happens when you evaluate your ODE function at t = 0 and this y0 vector? I'm betting that results in an all-zero vector.
What happens when you evaluate your ODE function at t = 1e-9 (for example) and the all zero vector? I'm betting that too results in an all-zero vector.
  2 Comments
Hassan Sultan
Hassan Sultan on 25 Oct 2020
Put some sponteneous emission to the initial fields. in the initial conditions, y0=[ 0 0 1e-6 1e-6];

Sign in to comment.

Categories

Find more on Graphics Objects in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by