How to integrate a set of coupled ODE's?
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
NUR SIDDIQ
el 18 de Ag. de 2019
Comentada: Star Strider
el 18 de Ag. de 2019
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
0 comentarios
Respuesta aceptada
Star Strider
el 18 de Ag. de 2019
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:
8 comentarios
Star Strider
el 18 de Ag. de 2019
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.)
Más respuestas (0)
Ver también
Categorías
Más información sobre Particle & Nuclear Physics en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!