Laser rate equation modelling using ODE45

Hi there, I have a question about modelling rate equations. Basically the modelling is solving differential equations hence I'm trying to use ODE45.

However, the solving seems running forever (3 days already). I was wondering if I'm using the wrong solver.

Can anyone give me some clues? Thanks a lot for your viewing and kindly help!

The code I wrote is as follows:

Function

function dy = rate_eq( t,y )
%Simplified rate equation for Yb ions.
dy = zeros(5,1);
h = 6.62606957E-34;
% Planck constant (m^2kg/s)
c = 299792458;
% Speed of light (m/s)
lambdaS = 1030E-9;
% Wavelength of signal light (m)
lambdaP = 975E-9;
% Wavelength of pump light (m)
vS = c/lambdaS;
% Frequency of signal (Hz)
vP = c/lambdaP;
% Frequency of pump (Hz)
A21 = 1/(600E-6);
% Spontaneous emission transition time from L2 -> L1 of 600 us
A20 = 1/(600E-6);
% Spontaneous emission transition time from L2 -> L1 of 600 us
W10 = 1/(500E-12);
% Fast decay transition time from L1 -> L0 of 500 ps
sigma21 = 0.5E-20;
% Stimulation emission cross section from L2 -> L1 of 0.5E-20 cm^2
sigma02 = 2.5E-20;
% Absorption cross section from L0 -> L2 of 2.5E-20 cm^2
dy(1) = sigma02*y(4)*y(3)/(h*vP)-A20*y(1)-A21*y(1)-sigma21*y(5)*y(1)/(h*vS)
% Change of population of upper lasing level (L2)
dy(2) = A21*y(1)+sigma21*y(5)*y(1)/(h*vS)-W10*y(2)
% Change of population of lower lasing level (L2)
dy(3) = -sigma02*y(4)*y(3)/(h*vP)+A20*y(1)+W10*y(2)
% Change of population of ground lasing level (L2)
dy(4) = (-sigma02*y(4)*y(3)/(h*vP)+A20*y(1))*(h*vP)
% Change of intensity of pump beam
dy(5) = (sigma21*y(5)*y(1)/(h*vS)+A21*y(1))*(h*vS)
% Change of intensity of signal beam
end

Main Script

TSPAN = [0 1];
Y0 =[0 0 8.75E20 1 0.2];
%y(1) = population of Upper Lasing level
%y(2) = population of Lower Lasing level
%y(3) = population of Ground Level
%y(4) = Optical intensity of Pump W/cm^2
%y(5) = Optical intensity of Signal W/cm^2
[T,Y] = ode45(@rate_eq,TSPAN,Y0);
subplot(3,1,1)
plot(T,Y(:,3),T,Y(:,2),T,Y(:,1))
legend('Ground level','Lower lasing level','Upper lasing level')
title('Population different levels')
subplot(3,1,2)
plot(T ,Y(:,4))
title('Pump Intensiry') % Pump intensity
subplot(3,1,3)
plot(T ,Y(:,5))
title('Singal Intensity') % Signal intensity

The equations correspond to the energy levels and transitions of a photoluminescence system showing below:

where

S0.2 is the absorption from ground level (L0) to upper lasing level (L2).

A2.0 is the spontaneous emission from upper lasing level (L2) to ground level (L0).

A2.1 is the spontaneous emission from upper lasing level (L2) to lower lasing level (L1).

SE2.1 is the stimulated emission from upper lasing level (L2) to lower lasing level (L1).

W1.0 is the fast non-radiative decay from lower lasing level (L1) to ground level (L0).

The rate equations are:

where

n2, n1 and n0 are the populations of the state in upper lasing level (L2), lower lasing level (L1) and ground level (L0), respectively.

σ02 and σ21 are the absorption and stimulation emission cross section of S0.2 and SE2.1 transitions.

A20 and A21 is the Einstein coefficient of transition A2.1 and A2.0.

W10 is the rate of the fast transition W1.0.

h is the planck's constant, vp and vs are the frequencies of the pump laser and signal laser. Ip and Is are the intensities of the pump and signal laser.

Answer 1

[Star Strider]

I’m surprised that you’re getting nothing out, since you left the semicolons off the ends of the lines in your ODE functions. Every time ‘rate_eq’ is evaluated, the results of those lines should display to the Command Window. Also, I suggest using the ‘dot operator’ to ‘vectorize’ them so that they do element-by-element operations instead of matrix operations (since it seems you are not doing matrix operations with them), for example:

dy(1) = sigma02.*y(4).*y(3)./(h.*vP)-A20.*y(1)-A21.*y(1)-sigma21.*y(5).*y(1)./(h.*vS);

and so with the others. It may not always be necessary, but I do it as a general rule.

Whenever ode45 gets stuck, I use ode15s. I suggest changing to it.

[Xingzhao]

Hi Star Strider,

I'm sorry about the misleading text. I meant the result when ode45 function finished its calculation. I'll try your advise and see what happens.

Thanks a lot!

[Star Strider]

My pleasure!

With ode15s it ran in less than half a second:

[T,Y] = ode15s(@rate_eq,TSPAN,Y0);

I also suggest changing to:

subplot(3,1,1)
semilogy(T,Y(:,3),T,Y(:,2),T,Y(:,1))

and adjusting the y-axis limits to get the best resolution, or using a linear y-scale and scaling the various values so ‘Ground Level’ doesn’t dominate the plot.

[Xingzhao]

Yep, that works! Thanks a lot!

The 'Ground Level' dominates the plot is probably due to the low input power. Once I enhance the input power it should be OK.

Really appreciate your help!

[Star Strider]

My pleasure!

I’m going to go back later and see if I can understand what you’re doing. It seems fascinating.

[Xingzhao]

It's a model of a piece of solid state laser material. Basically it describes this process: the material received energy from pump laser, stored the energy for a while, and then emitted light. Part of the light is from spontaneous emission (i.e. A2.0 and A2.1), part of the light is from stimulated emission (i.e. SE2.1).

From this modeling I want to see the amplification of the signal laser, which is related to transition A2.1 and SE2.1. The energy is stored in the dopants of the material, the population of the dopants in different energy state is n2, n1 and n0.

[Star Strider]

Thank you! Solid-state physics is far from my areas of expertise, so I appreciate the opportunity to learn something new.

[David Zare]

hi Xingzhao, did you get your answer?

Answer 2

[Mustafa Codur]

hey, I started it but it didnt give answer does it take so much time?

[Bahri Eren Uzuner]

Chill out bro, why so serious?