System of equations: Two different versions of Matlab, with the same script, give different solutions

Good afternoon,

I am writing because I noticed that the same script gives different results when using two Matlab versions. In particular, I tested it with Matlab versions 2017b and 2019a. The script calculates the solution of a system of two laser rate equations in the steady state as follows:

%RATE EQUATIONS IN THE STEADY STATE
PARAMETERS_ISSUE        % Runs the parameter's script. 
syms Curr n nph         % Definition of the parameters related to the injected current Curr, the carrier number n and the photon number nph

%  System of equations
eqn1=eta_inj.*Curr./q-n.*vs.*A./V-gamma_m.*((n-N0).*nph+n)-n./tau_l-C.*n.^3==0;
eqn2=gamma_m.*((n-N0).*nph+n)-nph./tau_p==0;
sol=solve([eqn1, eqn2],[n nph]);  %Solution of the equations

% Plot of the solutions

Curr=logspace(-6,-3,100);        %Vector with current values (x-axis)

%All the different solution:
carr1=sol.n(1);     
carr2=sol.n(2);
carr3=sol.n(3);
carr4=sol.n(4);

%Plots ONLY OF THE CARRIER NUMBER AS AN EXAMPLE
figure (1)  
hold on
loglog(Curr,subs(carr1,'Curr',Curr),'Linewidth',2)
loglog(Curr,subs(carr2,'Curr',Curr),'Linewidth',2)
loglog(Curr,subs(carr3,'Curr',Curr),'Linewidth',2)
loglog(Curr,subs(carr4,'Curr',Curr),'Linewidth',2)
ylabel('Carrier number') 
title('Carrier number')
legend('Solution 1','Solution 2','Solution 3','Solution 4')
xlabel('Current (A)') 
ax = gca;
ax.FontSize = 17.5;
hold off

Where the script "PARAMETERS_ISSUE " in line 2 contains all the parameters, and it is given by:

% CONSTANTS AND PARAMETERS
hbar=1.0546.*10.^-34;           % Planck's constant, J*s
q=1.6.*10.^-19;                 % Unit charge, C
m0=9.1*10^-31;
Mm0=25.30*q;
c=3*10.^8;                        % Speed of light, m/s
epsilon_0=8.85.*10.^-12;        % Vacuum permittivity, F/m
epsilon_r=11.6;                 % Dielectric constant of the InGaAs, adim.
lambdacav=1550.*10.^-9;          % Cavity wavelength, m
omega_cav=2*pi*c/lambdacav;      % Cavity angular frequency, rad/s
delta_lambda=100.*10.^-9;
delta_omega=delta_lambda*omega_cav/lambdacav;   % Gain medium homogeneous broadening in ang. freq., rad/s
tau_p=0.19.*10.^-12;       % Cavity photon lifetime (it can be calculated from Q/omega_cav, with Q=q factor), s
beta=0.40;
dif=sqrt(3*q^3/(2*omega_cav^2*m0)*Mm0/q)              % Dipole moment of the lasing transition, C*m
C=8.5.*10.^-41;                 % Auger recombination coefficient, m^6/s
eta_inj=1;                      % Injection efficiency, adim.
%    1 - CYLINDRICAL NANOSTRUCTURE: Characterized by radius and height. 
diam=0.28.*10.^-6;                     % Diameter, m
radius=diam/2;                       % Radius, m
height=1.29.*10.^-6;                 %Height of the cylindrical nanopillar, m
V=pi.*(diam/2).^2.*height       ;     %Volume of the cylindrical nanopillar, m^3 
A=2.*pi.*radius.*(height+radius) ;   %Total area of the cylindrical, m^2
Veff=0.75*V ;                        % Mode volume, m^3

vs=7.*10.^2;                     % Surface velocity non-passivated InGaAs, m/s

N0=(1.1*10^24).*V;                   % Number of carriers at transparency (calculated from the quasi-Fermi energies, note that it is NOT a density as it is multiplicated by the Volume), adim.      
gamma_m=2.*omega_cav.*dif.^2./(hbar.*epsilon_0.*epsilon_r.*Veff*delta_omega);
tau_l=(gamma_m.*(1-beta)./beta).^(-1)

When I run the script with Matlab 2017b, I obtain four pair of solutions that are correct and expected. The following picture displays the plot of the solutions of the variable n (for semplicity I am not showing nph):

On the other hand, when I run the same script with Matlab 2019a, I obtain the following:

Where there is a jump from one solution to the other that makes this not useful.

Moreover, in Matlab 2019a I get the warning "Warning: Imaginary parts of complex X and/or Y arguments ignored", which I think is the cause of these jumps, but unfortunately I could not solve this problem by checking other questions in the support section.

Could you please help me to let this script run in any version of Matlab without these jumps?

Thanks in advance, and kind regards.

Answer 1

[Walter Roberson]

• The exact solutions to those equations are badly non-linear, and will give incorrect answers for less than about 60 digits (default is 32)
• MATLAB has never defined what the index order means for a root() with an index, which is the structure being generated when you solve the two equations for n and nph in terms of Curr. Most of the time the index ordering appears to have to do with a sorting by absolute value and phase angle, but not always so.
• In any set of higher order polynomials with complex roots, it is common for the index order to appear to switch when you use root() and vary a parameter such as Curr. "switch" in the sense that you will be following one root, and then there will be a jump and you seem to follow another root. This is confusing, but it is not a "bug" precisely because the meaning of the order was never defined. All that you can really say about indexing with root() is that for any given release and any given set of parameter values, the result will be self-consistant -- but as you vary the parameter values then the order can look inconsistent to us.
• You are using a 4th order polynomial. You can deal with the situation by using 'MaxDegree', 4 in solve(), which will force MATLAB to generate the full closed-form formulas for the four roots. These expressions are slow to evaluate, and have the precision problems I mentioned above, but they will not have the jumps.
• As far as I can tell, your equations really do generate complex roots for two of the four solutions. Including for the first value of Curr, 1e-6 . In some cases the complex components might look like they are just round-off problems (such as 1e-27), but they presist even when evaluated to 1000 digits. In other cases, the complex component are about 1e+25 when evaluated to enough precision.
• When you do not evaluate to enough precision, the values you get out can be wildly wrong. Including appearing to be real-valued when the reality is that the true complex component is large.

[DPel]

Hello Walter,

Thanks for this feedback. Indeed, by adding 'MaxDegree' it works and I observe no "jumps" among the solutions, but as you mentioned the script is really slow. This aspect triggered my curiosity, as I do not understand how an older version of Matlab can handle this problem with short computation times (10 seconds), while this is not true for more recent versions. Do you have a clue about this, and a piece of advice to shorten a bit the computation time for my specific problem?

Thanks again

Answer 2

[Matt J]

Seems to me that you are getting the same solutions for each fixed Current level, but presented to you in a different order depending on the Matlab version. Can't you just sort() them into the order that you prefer?

[Matt J]

E.g.,

solutions=[subs(carr1,'Curr',Curr);...
            subs(carr2,'Curr',Curr);...
            subs(carr3,'Curr',Curr);...
            subs(carr4,'Curr',Curr)].';
        
loglog(Curr(:),sort(solutions,2))        

[DPel]

Hello Matt,

Thanks for your feedback. Nice idea, but unfortunately this script does not join properly the solutions..

[Matt J]

Yet, I find that my approach reproduces your original plot quite well,

%Plots ONLY OF THE CARRIER NUMBER AS AN EXAMPLE
figure (1)  
hold on
loglog(Curr,subs(carr1,'Curr',Curr),'Linewidth',2)
loglog(Curr,subs(carr2,'Curr',Curr),'Linewidth',2)
loglog(Curr,subs(carr3,'Curr',Curr),'Linewidth',2)
loglog(Curr,subs(carr4,'Curr',Curr),'Linewidth',2)
ylabel('Carrier number') 
title('Carrier number')
legend('Solution 1','Solution 2','Solution 3','Solution 4')
xlabel('Current (A)') 
ax = gca;
ax.FontSize = 17.5;
hold off

h=findobj(ax,'Type','line','-or','Type','functionline');

figure(2)

solutions=vertcat(h.YData).';   
plot(h(1).XData,sort(solutions,2,'descend'))  

ylabel('Carrier number') 
title('Carrier number')
legend('Solution 1','Solution 2','Solution 3','Solution 4')
xlabel('Current (A)') 

[DPel]

Yes it actually works well. My knowledge of the sort command is limited so I was not able to generate such a script. Thanks a lot!