Problem solving laser rate equations with ode45 command

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!!

Answer 1

[Star Strider]

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.

[mohammad heydari]

Thanks for the answer.

But my main problem is that I have to have four columns of numbers even if they have errors.

ode has to solve four unknowns, but that doesn't happen.

[Star Strider]

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]

very thanks :Star Strider , Steven Lord

The curve I put below is my favorite curve, which is not a logarithm curve. This curve is extracted from a valid article.

It has only been normalized to N۰. I have tried to select the most correct values for the parameters but this difference can be considerable.

[Star Strider]

Our pleasure!

I know essentially nothing about laser physics, so I cannot help you with respect to the parameter choices. There are methods to estimate them from data, if you have data that you want to fit to your differential equaations.

[mohammad heydari]

Thank you very much.

I have only two questions about your suggested edits

1. What effect does eps have on the initial condition?

2. are abs and semilogy just a matter of scaling to see results because I draw with plot command and getting different response from semilogy.

[Star Strider]

As always, my pleasure!

1. For the present purposes, the eps funciton returns a very small non-zero number. It prevents the zero values from propagating through the integration, producing zeros in the output matrix.

2. The abs function returns the amplitude of the complex elements of ‘y’, and to be compatible with the semilogy plot. The semilogy plot simply presents the results in a way that makes the individual vectors visible, and is the reason I used it. It is not necessary to use semilogy. However to use semilogy to its best advantage, use abs(y) or real(y), or real(y) and imag(y), depending on what you want in the plot. You would expect to get different results from plot and semilogy. Choose the one that presents your data as you want them to be presented.

[mohammad heydari]

Thank you very much for your help

Best regards

[Star Strider]

As always, my pleasure!

Answer 2

[Steven Lord]

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.

[mohammad heydari]

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]

Put some sponteneous emission to the initial fields. in the initial conditions, y0=[ 0 0 1e-6 1e-6];