How to integrate a set of coupled ODE's?

Hi everyone

I want to calculate pulse energy (Eout) of a q-switched laser, the equation is below:

Meanwhile phi(t) is depend to 3 coupled equation in laser rate equation [y(1)]. I try to use this code to solve the integral equation:

g = @(t) y(1);
phi = integral (g,0,inf,'ArrayValued',1)

Meanwhile after running the code, Matlab still continue calculating more than 2 hours! Anyone know the problem? What I want to get is single value of phi.

This is my full code, you can run by yourself.

Thank you

clear
clc
ti = 0; 
tf = 12E-6; 
tspan=[ti tf]; 
y0=[0.1; 0; 0];
[t,y] = ode45(@rate_eq,tspan,y0); 
%y(1) = Photon Density
%y(2) = Inverted Population Density
%y(3) = Photocarrier Density

subplot(3,1,1);
plot(t,y(:,1)) 
title('IPhoton Density'); 
xlabel('Time'); 
ylabel('Photon Density');
subplot(3,1,2);
plot(t,y(:,2)); 
title('Inverted Population Density '); 
xlabel('Time'); 
ylabel('Inverted Population Density'); 
subplot(3,1,3)
plot(t,y(:,3)); 
title('Photocarrier Density'); 
xlabel('Time'); 
ylabel('Photocarrier Density'); 

 

function dy = rate_eq( t,y )
%Rate equation for Q-switched fiber laser.
dy = zeros(3,1);
% Planck constant (m^2kg/s)
h = 6.62606957E-34;
% Speed of light
c = 299792458;
% Dissipative optical loss
delta = 0.4;
%Inversion reduction factor
gamma = 1.8;
% Length of total fiber (m)
lr = 15;
% Length of active fiber (m)
L = 3;
% Thickness of the SA (m)
Lsa = 1E-5;
% Output coupling ratio
R = 0.95;
% Stimulated emission area (m2)
sigmaes = 2E-25;
% Spontataneous decay time(t)
tg = 1E-2;
% Saturation photocarrier density 
Nsa = 2.4E27;
% SA carrier recombination time (s) 
tsa = 0.1E-9;
% nonsaturable loss (s) 
lambdans = 0.4;
% modulation depth 
lambdas = 0.1;
% Effective doping area of active fiber (m2) 
A = 1.26E-11;
% Pumping power (W)
Pp = 2000E-3;
% Wavelength of signal light (m)
lambdaP = 1530E-9;
% Frequency of signal (Hz)
v = c/lambdaP;
% Round trip transit time (Hz)
tr = lr/c;
% Pumprate of active medium
Wp = Pp/(h*v*A*L);
% Saturable absorption
lambdana = (lambdas/(1+(y(3)/Nsa)))+lambdans;
%lambdana = 0.5;
% Change of Photon density
dy(1) = (y(1)/tr)*(2*sigmaes*y(2)*L-2*lambdana*Lsa+log(R)-delta);
% Change of Population Inversion density
dy(2) = Wp-gamma*sigmaes*c*y(1)*y(2)-(y(2)/tg);
% Change of Photocarrier density
dy(3) = c*y(1)*lambdana-(y(3)/tsa);
%Pulse Energy
g = @(t) y(1);
phi = integral (g,0,inf,'ArrayValued',1)
end

Answer 1

[Star Strider]

Your system is ‘stiff’, so use one of the stiff solvers instead of ode45:

[t,y] = ode15s(@rate_eq,tspan,y0); 

With this change, you code ran with:

Elapsed time is 19.170238 seconds.

(on my machine), and produced this plot:

[NUR SIDDIQ]

Thanks for your answer Star Strider

But I get an error

Error: File: ode15s.m Line: 618 Column: 13

"rate" previously appeared to be used as a function or command, conflicting with its use here as the name of a

variable.

A possible cause of this error is that you forgot to initialize the variable, or you have initialized it implicitly

using load or eval.

Error in laser (line 7)

[t,y] = ode15s(@rate_eq,tspan,y0);

I used Matlab R2019a with student license.

[Star Strider]

My pleasure.

You need to either run all your code inside a separate function (that is what I did), or save your ‘rate_eq’ function as a separate .m file as rate_eq.m on your MATLAB user path.

Either of those should work.

[NUR SIDDIQ]

I try both of your suggestion, run all code and save as separate .m file.

But there is a problem

>> laser
Error: File: laser.m Line: 1092 Column: 1
The function "rate_eq" was closed with an 'end', but at least one other function definition was not. All functions
in a script must be closed with an 'end'.

But however, I just succesfull using

[t,y] = ode23s(@rate_eq,tspan,y0); 

But thats not the point of my question, i want to calculate the single value of phi (after integrated.

Here the result

Warning: Minimum step size reached near x = 3.24519e+30. There may be a singularity, or the tolerances may be too
tight for this problem. 
> In integralCalc/checkSpacing (line 467)
  In integralCalc/iterateArrayValued (line 215)
  In integralCalc/vadapt (line 130)
  In integralCalc (line 83)
  In integral (line 88)
  In laser>rate_eq (line 84)
  In odenumjac (line 131)
  In ode23s (line 337)
  In laser (line 7) 

And with phi value that changes overtime

Do you have an idea how to solve this problem? Thank you

[Star Strider]

The function I use to run code such as yours does not allow me to use persistent variables inside ‘rate_eq’, so I simply created another function to calculate ’phi’, using the solved values for ‘t’ and ‘y’:

function phi = phi_rate_eq( t,y )
%Rate equation for Q-switched fiber laser.
dy = zeros(3,1);
% Planck constant (m^2kg/s)
h = 6.62606957E-34;
% Speed of light
c = 299792458;
% Dissipative optical loss
delta = 0.4;
%Inversion reduction factor
gamma = 1.8;
% Length of total fiber (m)
lr = 15;
% Length of active fiber (m)
L = 3;
% Thickness of the SA (m)
Lsa = 1E-5;
% Output coupling ratio
R = 0.95;
% Stimulated emission area (m2)
sigmaes = 2E-25;
% Spontataneous decay time(t)
tg = 1E-2;
% Saturation photocarrier density 
Nsa = 2.4E27;
% SA carrier recombination time (s) 
tsa = 0.1E-9;
% nonsaturable loss (s) 
lambdans = 0.4;
% modulation depth 
lambdas = 0.1;
% Effective doping area of active fiber (m2) 
A = 1.26E-11;
% Pumping power (W)
Pp = 2000E-3;
% Wavelength of signal light (m)
lambdaP = 1530E-9;
% Frequency of signal (Hz)
v = c/lambdaP;
% Round trip transit time (Hz)
tr = lr/c;
% Pumprate of active medium
Wp = Pp/(h*v*A*L);
% Saturable absorption
lambdana = (lambdas/(1+(y(3)/Nsa)))+lambdans;
%lambdana = 0.5;
% Change of Photon density
dy(1) = (y(1)/tr)*(2*sigmaes*y(2)*L-2*lambdana*Lsa+log(R)-delta);
% Change of Population Inversion density
dy(2) = Wp-gamma*sigmaes*c*y(1)*y(2)-(y(2)/tg);
% Change of Photocarrier density
dy(3) = c*y(1)*lambdana-(y(3)/tsa);
%Pulse Energy
g = @(t) y(1);
phi = integral (g,0,inf,'ArrayValued',1)
end

for k = 1:numel(t)
    phiv(k) = phi_rate_eq( t(k), y(k,:) );
end

figure
plot(t, phiv)
grid
xlabel('t')
ylabel('\phi(t)')

[NUR SIDDIQ]

You are really help me Sir! Thank you

But still there is a problem, this is the problem:

Error using odearguments (line 95)
PHI_RATE_EQ returns a vector of length 1, but the length of initial conditions vector is 3. The vector returned
by PHI_RATE_EQ and the initial conditions vector must have the same number of elements.

Error in ode23s (line 121)
    = odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in coba4 (line 7)
[t,y] = ode23s(@phi_rate_eq,tspan,y0);
 
>> 

This is my full code

clear
clc
ti = 0; 
tf = 10E-6; 
tspan=[ti tf]; 
y0=[0.01; 0; 0];
[t,y] = ode23s(@phi_rate_eq,tspan,y0); 
%y(1) = Photon Density
%y(2) = Inverted Population Density
%y(3) = Photocarrier Density
 
for k = 1:numel(t)
    phiv(k) = phi_rate_eq( t(k), y(k,:) );
end
figure
plot(t, phiv)
grid
xlabel('t')
ylabel('\phi(t)')

function phi = phi_rate_eq( t,y )
%Rate equation for Q-switched fiber laser.
dy = zeros(3,1);
% Planck constant (m^2kg/s)
h = 6.62606957E-34;
% Speed of light
c = 299792458;
% Dissipative optical loss
delta = 0.4;
%Inversion reduction factor
gamma = 1.8;
% Length of total fiber (m)
lr = 15;
% Length of active fiber (m)
L = 3;
% Thickness of the SA (m)
Lsa = 1E-5;
% Output coupling ratio
R = 0.95;
% Stimulated emission area (m2)
sigmaes = 2E-25;
% Spontataneous decay time(t)
tg = 1E-2;
% Saturation photocarrier density 
Nsa = 2.4E27;
% SA carrier recombination time (s) 
tsa = 0.1E-9;
% nonsaturable loss (s) 
lambdans = 0.4;
% modulation depth 
lambdas = 0.1;
% Effective doping area of active fiber (m2) 
A = 1.26E-11;
% Pumping power (W)
Pp = 2000E-3;
% Wavelength of signal light (m)
lambdaP = 1530E-9;
% Frequency of signal (Hz)
v = c/lambdaP;
% Round trip transit time (Hz)
tr = lr/c;
% Pumprate of active medium
Wp = Pp/(h*v*A*L);
% Saturable absorption
lambdana = (lambdas/(1+(y(3)/Nsa)))+lambdans;
%lambdana = 0.5;
% Change of Photon density
dy(1) = (y(1)/tr)*(2*sigmaes*y(2)*L-2*lambdana*Lsa+log(R)-delta);
% Change of Population Inversion density
dy(2) = Wp-gamma*sigmaes*c*y(1)*y(2)-(y(2)/tg);
% Change of Photocarrier density
dy(3) = c*y(1)*lambdana-(y(3)/tsa);
%Pulse Energy
g = @(t) y(1);
phi = integral (g,0,inf,'ArrayValued',1)
end

[Star Strider]

Do not use ‘phi_rate_eq’ with the ODE solver. Your equations have already been solved, so just use ‘phi_rate_eq’ in the for loop with the solved variables to calculate and return ‘phi’ (as ‘phiv’).

Do this instead:

ti = 0; 
tf = 10E-6; 
tspan=[ti tf]; 
y0=[0.01; 0; 0];
[t,y] = ode23s(@rate_eq,tspan,y0); 
%y(1) = Photon Density
%y(2) = Inverted Population Density
%y(3) = Photocarrier Density
 
for k = 1:numel(t)
    phiv(k) = phi_rate_eq( t(k), y(k,:) );
end
figure
plot(t, phiv)
grid
xlabel('t')
ylabel('\phi(t)')

That should do what you want.

[NUR SIDDIQ]

Thank you very much!

Thats what I want :)

I just wonder why R2018b is better than R2019a for a ode15s

My friend also success like you when using command ode15s, but I dont

[Star Strider]

As always, my pleasure!

I have never had any problems with ode15s. If you are having problems with it, be sure you are using the latest MATLAB update, using ‘Check for Updates’. In R2018b and earlier, this was in the ‘Add-Ons’ tab in the top toolstrip. In R2019a, it moved to ‘Help’. (Update 5 to R2019a was just released. I installed it yesterday.)