Laser rate equation(ODE) simulation problem

[oronno]

hello...

Can anyone tell me how can I simulate the laser rate equation with matlab? i allready wrote a code whics is as follows,

%have to write one editor

function dy=rate_equation(t,y)
dy = zeros(2,1);
I =50e-3;

%for t=1:10;

q =1.602e-19;
v =9e-11;
eps =3.4e-17;
N0 =1.8e-18;
tn =3e-9;
g0 =3e-6;
lamda=0.44;
bita=4e-4;
tp=1e-12;
dy(1)=(I/(q*v))-((g0*(y(1)-N0)*y(2))/(1+eps*y(2)))-(y(1)/tn);
dy(2)=(lamda*g0*(y(1)-N0)*y(2))/(1+eps*y(2))-(y(2)/tp)-((lamda*bita*y(2))/tn);
end

%have to write other editor
%step size=0.01;
clc
tspan=[0,10];
y0=[0,1];
[T,Y]= ode45(@try1,tspan,y0);
plot(T,Y(:,1),'-',T,Y(:,2),'.');
title('plot of carrier and photon densities');
xlabel('time');
ylabel('densities');

 

[matematikawan]

So what really is the problem? You have written the code. Have you run it?

 

[oronno]

bro...i wrote the code and tried to run it. but every time i try to debug it...the MATLAB gets busy and no result comes out. seems there's something wrong with the code. so, i need someone to help me out

 

[matematikawan]

Actually I don't have slightest idea the physical meaning of your equation and I'm also quite new to Matlab. Only able to write small script file.

From your equations, it looks as if you have three dependent variables N, S, and I. Do you assume I(t) to be constant here?

In your script file you use the anonymous function @try1 whereas your function m-file is rate_equation.

Hopefully you have type your equations commands correctly.

Is there any chance you dimensionless all your variables?

 

[oronno]

thanx bro for yoyr reply. here S and N are two dependent variables. they are dependent on each other. and I here is constant. or we can vary I in any range ( like 50-70 mA , or 20-50 mA). and you are right on the @try1. it should be @rate equation. actually I tried with different different m files and mistakenly I pested @tri1 here.

and I didnt understand what did you mean by your last line (dimentionless variables).

thanx for your concern and I am waiting for your reply.

 

[matematikawan]

May be dimensionless the variables are not important in this problem. Dimensionless doesn't make the problem simplier, it only us to control what values to generate. I think.
I used to work on diffusion problem before, solving pde for material concentration c(x,t). We didn't solve the equation numerically but analytically. A lot of linearization has to been done.
We knew, such as, the maximun length material could diffused x₀, etc.
Rather than work with the variable x, we used the variable $$y=\frac{x}{x_0}$$ which is dimensionless. We then generate only the value 0<y<1. Similarly were done to other quantities t, c etc.

In your program, what make you chose t=[0 10] ? May be your intuition.
Your equations are nonlinear. Value need to be careful chosen so that the solution converges. If you think you have done everything correct but still don't get any result, probably the subroutine ode45 is not suitable for solving your problem. We don't know the algorithm behind ode45. Try other subroutine.

 

[billboss]

Hi!
I'm working successfully with rate equations for lasers with matlab.
Your code, oronno, is ok. The ode45 routine of matlab, which uses Runge-Kutta method, is good for our kind of non-linear differential equations. So, is it ok? No, I think there's a kind of numerical instability of algoritm when you choose wrongs initial conditions(a good idea is 0 for carriers and 0 for photons). Then, maybe there's a problem that regard your formulation of the rate equations. I tried one like your but i had no answers from MATLAB like you.
So, i searched for another formulation and ...WOW .
I append the current version,that's a little bit redundant, you can modify it!
open sources, open minds

This is the "main" script (you can convert in a function).

clc
clear all
close all

tau_s = 3e-9;
N0 = 1e24;
A =1e-12;
P0 = 1/(A*tau_s);

TSPAN = [0 10];
Y0 =[0 0];

[T,Y] = ODE45(@rate_equation,TSPAN,Y0);


subplot(2,1,1)
plot(T*tau_s ,Y(:,1)*N0)
title('Densità dei portatori nel livello laser superiore') % carriers density in high laser level

subplot(2,1,2)
plot(T*tau_s ,Y(:,2)*P0)
title('Densità dei fotoni nel mezzo attivo') % photons density in activer region


And this is the function.

function dy = rate_equation(t,y)
dy = zeros(2,1);

tau_s = 3e-9; % carriers lifetime
tau_p = 1e-12; % photons lifetime
A = 1e-12; % linear gain costant
N0 = 1e24; % trasparency carries density
V = 3.75e-14; % modal volume
gamma = 1e-5; % gain compression factor
q = 1.6e-19; % electron charge

I0 = N0*q*V/tau_s; % trasparency current
tau_norm = tau_s/tau_p;
eta = A*tau_p*N0; % efficiency

I = 2.5*I0; % pumping current ( try: from I0 to 3*I0 for example ...and see what happens!)

dy(1)= I/I0 -y(2)*(y(1) - 1) -y(1);
dy(2) = tau_norm*(y(2)*(eta*(y(1) - 1) -1) + gamma*eta*y(1));

% this is a very simple (dimensionless) and efficient version of the rate equations (note you
% can add another equation, dy(3)(not coupled), that gives you the phase(t) of your laser!

Sorry for my english!
I hope the world will live as one.

 

[oronno]

thanx to both of u my brothers. I am out of my stations right now... i will go back to my place tonight and i will try with your code. and offcourse I'll let u know. thank u once again for your concern. i am really greatful to u.

 

[billboss]

Some explanations:

The code is based on a rate equation dN/dt which shows only a linear dipendence from N (or from (N - N0)), where N0 is the trasparency density:

dN/dt = I/(qV) - A(N - N0) P - N/tau_s

where A is the linear gain of the laser medium, tau_s is the carrier lifetime(before spontaneous de-excitation). Here, I/(qV) represents the positive pumping rate, while A(N -N0)P and N/tau_s represent ,respectively, the negative absorption rate(which depends from incident photons density P) and stimulated emission rate(which depends from density and lifetime of the excited-level carriers).
Then, the second rate equation dP/dt will show similar terms:

dP/dt = A(N -N0)P -P/tau_p + gamma N/tau_s

where tau_p is the photon lifetime in the laser and gamma is the spontaneous emission coefficient(dimensionless).

In these equations:

= A
[A] = m^3/s
[V] = m^3
[tau_s] =[ tau_p] = s
[N] = [N0] = m^(-3)
[gamma] = 1

For Matlab implementation, it's a good idea to normalize the variables and, of course, the other parameters. The equations become:

dn/dT = I/I0 - p(n-1) - n

dp/dT = tau_norm * (p (eta(n -1) -1) + gamma*eta*n)


where

n = N/N0
T = t/tau_s
I0 = N0*q*V/tau_s
P0 = 1/(A*tau_s)
p = P/P0
tau_norm =tau_s/tau_p
eta = A*tau_p*N0

You can try the values of the parameters written in the code:

tau_s = 3 x 10^(-9)
tau_p = 1 x 10^(-12)
A = 1 x 10^(-12) [m^3/s]
N0 = 1 x 10^24 [m^(-3)]
V = 4 x 10^(-14) [m^3]
gamma = 1 x 10^(-5) []

I hope it's all clear. bye

 

[oronno]

thanx bro...thank u very much. i am greatful to u. you helped me a lot. i have one more question dear. can u give me a reference of this equation? I mean, from the paper u got these parameters. please let me know if u have any. that would be very kind of you. thanks again.

 

[oronno]

bro...I am waiting for your reply. can u please give me the link of the equation or the reference of this equation? please help me out. thank u

 

[snatcos]

thanks to everyone : it's helpfull to share code. billboss, your code is just i was trying to obtain in matlab. I'm new with MATLAB programming and your code is an excellent starting point for the model i want to implement. thanks again

 

[snatcos]

... just one thing : from your rate eqs. gamma is the spontaneous emission coefficient(dimensionless) as you say in youor second post and not the compression factor, as said in the first

 

[meow44]

simulation rachat de credit​

thanks for the post! it's very informative

 

[fahad305]

the rate equation is
DN/Dt = I/qv - g0(N-N0)S/(1+epsilon*S) -N/Tow_s
i wanted to know why are we ignoring the term (1+epsilon*S)?

 

[AiRAVATA]

May be dimensionless the variables are not important in this problem. Dimensionless doesn't make the problem simplier, it only us to control what values to generate.

Thas far from truth. You should always adimensionalize your problem, as it will give you a general perspective on the system and (in nonlinear cases) important values to work with.

Imagine you want to model the behavior of a string of infinite length, tension $N$ and constant density $\rho$, with initial conditions $u(x,0) = L_0 f(x)$ and $u_t(x,0) = v_0 g(x)$, where $f(x),\,g(x)$ are dimensionless functions. Then the problem reads

$$\rho u_{tt}(x,t) - N u_{xx}(x,t) = 0,$$

$$u(x,0) = L_0 f(x), \qquad u_t(x,0) = v_0 g(x).$$

If you haven't stuidied this equation earlier, then there is no reason to believe that it will behave the same way for all parameters $\rho,\,N,\,L_0,\,v_0$. Now, if you take the dimensionless variables

$$t = T t^*,\qquad x = X x^*, \qquad u(x,t) = L_0 u^*(x^*,t^*),$$

then the problem transforms in

$$\frac{\rho}{T^2} u^*_{t^*t^*}(x^*,t^*) - \frac{N}{X^2} u^*_{x^*x^*}(x^*,t^*) = 0,$$

$$u^*(x^*,0) = f(x), \qquad u^*_t^*(x^*,0) = \frac{T v_0}{L_0} g(x).$$

Now, take

$$T = \frac{L_0}{v_0}, \hbox{ and } X = \frac{1}{c} \sqrt{\frac{NL_0^2}{\rho v_0^2}},$$

where $c$ is a dimensionless parameter, then, by loosing the asterisks, the system becomes

$$u_{tt} - c^2 u_{xx} = 0, \qquad u(x,0) = f(x), \qquad u_t(x,0) = g(x).$$

This is very important, as you can study any infinite string regardless of it's physicall parameters with the same equation. Furthemore, there is only one parameter to manipulate, instead of four. Once you've solved the wave equation, you only need to see how the solution varies with different c's, to know how the string will behave for different values of density, tension, etc.

I've portrayed a simple example here, but when dealing with nonlinear PDE's (or ODE's), the response of the system need not to be the same for different values of the dimensionless parameters, and it is crucial for the study of such systems, to analize the behavior of the PDE's for different values of the parameters. A classic example is the multiplicity of solutions. There are nonlinear PDE's depending on a parameter c that have two solutions for $c \in (0,c^*)$, one solution for $c=c^*$, and none for $c \in (c^*,\infty)$.

 

[matematikawan]

Wow. Thats cool AiRAVATA.

For analysis purposes I agree it is better to dimensionalize the equation. But when the programmer are only interested in the profile output, is it worthful to dimensionalize the equation first ?

 

[billboss]

Good morning!
Some of you, probably, hate me, because I gave my contribute, and then I flyed away from here...
But my university is a kind of nightmare, and I had a very bad time this year...
Anyway, I finished my work on laser DFB and I want to share it with us.
Sorry, the language of the document is italian...
Maybe, I'll translate it ... if there will be many questions!
Bye

 

[matematikawan]

At the beginning of the year Matlab was new to me. Please forgive me because my earlier posts were really an embarrassment. Didn't contribute to solving the problem. I just simply want to reply.
But today if I see a MATLAB code I can confidently interprete it.

 

[katla]

Hi all, I've another problem. How L-I curve of laser diode can be obtained using this relations? I tried to define different amounts of I in script program but i could not use this in Function program. Please help me ...

 

[billboss]

Hi, Katla, try to read pages #28-31 of my little work (file "Paragrafi+Appendici.pdf"): I used only the relations which are there to obtain the L-I curve...
If you need a translation from italian, contact me here: info@billboss.homelinux.org
What you have to do in Matlab is essentially:
- create a working directory
- create the function "rate_equation" and the script (page 30) in the working directory
- run the script (eventually edit the values of parameters and/or the range of current values you want to sweep)

 

[katla]

Dear billboss, thanks a lot for your attention and help :shy: