Problem solving laser rate equations with ode45 command

Question

mohammad heydari 2019-11-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/488970-problem-solving-laser-rate-equations-with-ode45-command

评论： Hassan Sultan 2020-10-25

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.

Even if the given parameters have errors for any reason, this should not cause zero output.

Thanks in advance for helping!!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2019-11-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/488970-problem-solving-laser-rate-equations-with-ode45-command#answer_399580

编辑：Star Strider 2019-11-3

在 MATLAB Online 中打开

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 个评论
显示 6更早的评论隐藏 6更早的评论

Star Strider 2019-11-3

在 MATLAB Online 中打开

As always, my pleasure!

Steven Lord pointed to another problem (that I did not initially appreciate), that being the zero initial conditions. Correct for that by adding eps to the existing ‘y0’ to get:

tspan = [0 2e-9];                   % time interval, up to 2 ns
y0 = [0,0,0,0]+eps;     
[t,y] = ode45(@rate_eq_fano_program_1,tspan,y0);
size(t);
t=t*1e9;
y = abs(y);                                         % Change As Necessary To Get The Desired Resullt
y_max = max(y);
y1 = y_max(1);
y2 = y_max(2);
y3 = y_max(3);
y4 = y_max(4);
d = semilogy(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', 'Location','SE')  % legend inside the plot

I use a semilogy plot and abs(y) here so the full range of the data (all four columns) are easily visible and distinguishable. Change abs(y) to real(y) or plot the real and imaginary parts separately, depending on the desired result.

mohammad heydari 2019-11-3

Thank you very much for your help

Best regards

Star Strider 2019-11-3

As always, my pleasure!

请先登录，再进行评论。

Answer 2

Steven Lord 2019-11-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/488970-problem-solving-laser-rate-equations-with-ode45-command#answer_399581

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

mohammad heydari 2019-11-3

No problem in terms of timing. Surprisingly, the ordinary differential equations do not work properly and the 4 unknowns are not achieved while everything is structurally correct.

Hassan Sultan 2020-10-25