I have the following code:

(*Input parameters*)
fp = 6.0; (*Pump Frequency*)
ϕ0 = 207/100;
LL = 1/10;
a = 10;(*μm*)
r = 1/2;
Cj = 329/10^6;(*nF*)
CC = 39/10^6;(*nF*)

ωp = 2.0*π*fp;
detu[Δ_] := 2*π*Δ;
ωs[Δ_] := ωp - detu[Δ];
ωi[Δ_] := ωp + detu[Δ];
ωth = 3.0*ωp;
zchar = Sqrt[LL/CC];
I0 = ϕ0/(2.0*π*LL);
Ip = r*I0;
ω0 = 1.0/Sqrt[LL*CC];
RR = 10^4/(CC*ω0);
ωpt = ωp/ω0;
ωst[Δ_] := ωs[Δ]/ω0;
ωit[Δ_] := ωi[Δ]/ω0;
ωtht = ωth/ω0;
β = Cj/CC;
δ0 = 1.0/(RR*CC*ω0);
kpt = ωpt/Sqrt[1.0 - β*ωpt^2];
kst[Δ_] := ωst[Δ]/Sqrt[1.0 - β*(ωst[Δ])^2];
kit[Δ_] := ωit[Δ]/Sqrt[1.0 - β*(ωit[Δ])^2];
ktht = ωtht/Sqrt[1.0 - β*ωtht^2];
Ap0t = (Ip*zchar)/(ωp*LL*I0);
αpt = kpt^5/(16.0*ωpt^2)*Ap0t^2;
αtht = (kpt^2*ktht^3)/(8.0*ωtht^2)*Ap0t^2;

Γst[Δ_] := (kst[Δ]*δ0)/ωst[Δ];
Γit[Δ_] := (kit[Δ]*δ0)/ωit[Δ];

C2 = (kpt^3*(2.0*kpt - ktht)*ktht)/(16.0*ωpt^2);
D2 = -((kpt^4*ktht)/(16.0*ωtht^2));
U1[Δ_] := (kpt^2*kst[Δ]^3)/(8.0 ωst[Δ]^2);
U2[Δ_] := ((2.0*kpt - kit[Δ])*kpt^2*kst[Δ]*kit[Δ])/(16.0*(ωst[Δ])^2);
U3[Δ_] := ((kpt + kit[Δ] - ktht)*kpt*ktht*kst[Δ]*kit[Δ])/(8.0*(ωst[Δ])^2);

V1[Δ_] := (kpt^2*kit[Δ]^3)/(8.0 ωit[Δ]^2);
V2[Δ_] := ((2.0*kpt - kst[Δ])*kpt^2*kst[Δ]*kit[Δ])/(16.0*(ωit[Δ])^2);
V3[Δ_] := ((kpt + kst[Δ] - ktht)*kpt*ktht*kst[Δ]*kit[Δ])/(8.0*(ωit[Δ])^2);

ΔkL[Δ_] := 2.0*kpt - kst[Δ] - kit[Δ];
Δkth = 3.0*kpt - ktht + 3.0*αpt - αtht;

Ap[x_] := Ap0t*E^(I*αpt*x) + I*(C2*D2*Ap0t^5*E^(I*αpt*x))/Δkth*(x + (E^(-I*Δkth*x) - 1.0)/(I*Δkth));
Ath[x_] := D2/Δkth*Ap0t^3*E^(I*αtht*x)*(E^(I*Δkth*x) - 1.0)

(*Matrix Exponential*)
MM[xp_, Δ_] := Simplify[{{I*U1[Δ]*ComplexExpand[Abs[Ap[xp]]^2] - Γst[Δ]/2, I*U2[Δ]*E^(I*ΔkL[Δ]*xp)*ComplexExpand[(Ap[xp])^2] + I*U3[Δ]*Conjugate[Ap[xp]]*Ath[xp]*E^(I*(ktht - kpt - kst[Δ] - kit[Δ])*xp)}, {-I*V2[Δ]*E^(-I*ΔkL[Δ]*xp)*ComplexExpand[(Ap[xp])^2] - I*V3[Δ]*Ap[xp]*Conjugate[Ath[xp]]*E^(-I*(ktht - kpt - kst[Δ] - kit[Δ])*xp), -I*V1[Δ]*ComplexExpand[Abs[Ap[xp]]^2] - Γit[Δ]/2}}, Assumptions -> {xp >= 0, xp ∈ Reals}];
PP[xx_, Δ_] := Integrate[MM[y, Δ], {y, 0, xx}];
Smat[xx_, Δ_] := MatrixExp[PP[xx, Δ]];

Upon calling Smat[x,0], I get the following error:

enter image description here

Based on my preliminary search of the problem (here, and here), it seems that the issue usually stems from not rationalizing the input parameters. But my input parameters were all rationalized. What is going on here?

Edit: Here is a screenshot of me printing the matrix PP[x,0], which Mathematica has no problem printing for me. enter image description here

  • $\begingroup$ (1) It takes some time to produce the 2x2 matrix you are trying to exponentiate. It has a leaf count around 48K. I've rune Eigensystem on it for several minutes with no result thus far (and it is possible that MatrixExp is coughing on a different-but-related matrix). I'm using a development kernel for this computation so perhaps something got worse in 13.3 (the lack of speed) or perhaps better (no error messages thus far). If you have details about input size and time needed please add them to the post. $\endgroup$ May 16, 2023 at 21:08
  • $\begingroup$ (2) I do not know if this will make a difference to the outcome, but it does seem that several definitions in the code are in terms of inexact numbers. $\endgroup$ May 16, 2023 at 21:09
  • $\begingroup$ @DanielLichtblau is it possible to give any unofficial very rough estimate on a time frame when 13.3 might be released? $\endgroup$
    – Nasser
    May 16, 2023 at 22:19
  • $\begingroup$ @Nasser Not something I could say. Many moving parts... $\endgroup$ May 16, 2023 at 22:38
  • 1
    $\begingroup$ Could you please add the message name for the error message? $\endgroup$
    – Michael E2
    May 16, 2023 at 23:54

3 Answers 3


As I had thought, the issue is that Eigensystem (and, by extension, MatrixExp) cannot deduce there is a null vector corresponding to each eigenvalue. This in turn is because the presence of symbolic input forces the eigensystem computation to compute a null space by row reduction, and it is effectively impossible to determine that a large nonpolynomial symbolic expression containing approximate numbers is zero.

There are two reasonable ways I know of to proceed. I think the better one is to simply compute the exponential for different values of the variable xx and perhaps interpolate from that. The "direct" way, which I do not think is numerically stable and might suffer various other ills, is to provide Eigensystem with a non-default zero test.

Here is the matrix in question.

pmat = {{(-0.00005121161524020216 + 0.0026178649909336586*I)*xx + 
   (0. + 3.993655784881471*^-12*I)*xx^3 + 
   (0. + 4.000111699067347*^-9*I)*xx*
    Cos[0.07739713430827709*xx] + (0. + 0.0000590771943669321*I)*
  (1.5788343611162837*^-6*E^((0. + 0.0773971343082771*I)*xx) - 
    0.47712990394682414*E^((0. + 0.15217182721996186*I)*xx) - 
    0.05012997734894152*E^((0. + 0.1547942686165542*I)*xx) - 
    0.00001530065012797118*E^((0. + 0.22694652013164662*I)*xx) + 
    0.0017007732791776591*E^((0. + 0.2321914029248313*I)*xx) - 
    3.285851077155409*^-9*E^((0. + 0.30434365443992367*I)*xx) + 
    3.1763715703476475*^-9*E^((0. + 0.3095885372331083*I)*xx) - 
    (0. + 3.4842579271608996*^-6*I)*
     E^((0. + 0.1547942686165542*I)*xx)*xx + 
    (0. + 1.1418783646506823*^-7*I)*
     E^((0. + 0.2321914029248313*I)*xx)*xx + 
    (6.462222647503057*^-9*E^((0. + 0.306966095836516*I)*xx) + 
      E^((0. + 0.07477469291168475*I)*xx)*
       (-1.1557957707646625*^-24 - 
        (0. + 8.642427382755163*^-26*I)*xx) + 
      E^((0. + 0.22956896152823897*I)*xx)*
       (0.000029622594647121442 - 
        (0. + 1.9345090229089154*^-9*I)*xx) + 
      E^((0. + 0.15217182721996186*I)*xx)*(0.5255432008849645 - 
        (0. + 0.00006927475025946675*I)*xx - 
     Cos[0.0026224413965923324*xx] + 
    ((0. - 1.0947950680776157*^-10*I)*
       E^((0. + 0.306966095836516*I)*xx) + 
      E^((0. + 0.07477469291168475*I)*xx)*
       ((0. - 1.1557957707646625*^-24*I) + 
        8.642427382755163*^-26*xx) + 
      E^((0. + 0.22956896152823897*I)*xx)*
       ((0. - 9.787056088209194*^-7*I) + 1.9345090229089154*^-9*
         xx) + E^((0. + 0.15217182721996186*I)*xx)*
       ((0. + 0.5255432008849645*I) + 0.00006927475025946675*
         xx - (0. + 2.2843155561479443*^-9*I)*xx^2))*
   E^((0. + 0.15217182721996186*I)*xx)}, 
 {0.5740137499557736 - 0.05012997734894152/
    E^((0. + 0.0026224413965923454*I)*xx) - 
   0.5255432008849648*E^((0. + 0.0026224413965923454*I)*xx) - 
   0.000015300650127971185/E^((0. + 0.07477469291168473*I)*xx) + 
   1.5788343611162841*^-6*E^((0. + 0.07477469291168476*I)*xx) + 
   0.0016864755101056526/E^((0. + 0.0800195757048694*I)*xx) - 
   0.000014297769072006634*E^((0. + 0.08001957570486945*I)*xx) - 
   3.285851077155408*^-9/E^((0. + 0.15217182721996186*I)*xx) - 
   3.1763715703476467*^-9*E^((0. + 0.1574167100131465*I)*xx) + 
   ((0. + 3.4842579271609005*^-6*I)*xx)/
    E^((0. + 0.0026224413965923454*I)*xx) + 
   (0. + 0.00006927475025946678*I)*
    E^((0. + 0.0026224413965923454*I)*xx)*xx - 
   ((0. + 1.1418783646506826*^-7*I)*xx)/
    E^((0. + 0.0800195757048694*I)*xx) + 
   2.284315556147945*^-9*E^((0. + 0.0026224413965923454*I)*xx)*
    xx^2 + E^((0. + 0.0026224413965923454*I)*xx)*
    (9.78705608820919*^-7 + (0. + 1.9345090229089154*^-9*I)*xx)*
    Cos[0.07739713430827709*xx] + 
   (1.1319599590726744*^-26/E^((0. + 0.0026224413965923454*I)*
        xx) + 1.094795068077615*^-10*
      E^((0. + 0.0026224413965923454*I)*xx))*
    Cos[0.15479426861655418*xx] - 
   (0. + 0.000029622594647121452*I)*
    E^((0. + 0.0026224413965923454*I)*xx)*
    Sin[0.07739713430827709*xx] - 1.9345090229089154*^-9*
    E^((0. + 0.0026224413965923454*I)*xx)*xx*
    Sin[0.07739713430827709*xx] + 
   (0. + 6.679677322243623*^-25*I)*Sin[0.15479426861655418*xx] + 
   ((0. + 6.681595028036278*^-25*I)*Sin[0.15479426861655418*xx])/
    E^((0. + 0.0026224413965923454*I)*xx) - 
   (0. + 6.462222647503055*^-9*I)*
    E^((0. + 0.0026224413965923454*I)*xx)*
  (-0.00005121161524020216 - 0.0026178649909336586*I)*xx - 
   (0. + 3.993655784881471*^-12*I)*xx^3 - 
   (0. + 4.000111699067347*^-9*I)*xx*
    Cos[0.07739713430827709*xx] - (0. + 0.0000590771943669321*I)*

The zero assessment I'll use just plugs in a random real value for xx and checks if the result is beneath some threshold.

esys = Eigensystem[pmat,

This takes around a half second on my desktop machine and gives a full result for the eigenvectors.


Your problem is to get the matrix exponential of a 2x2 matrix with a free parameter t

      MatrixPower[ E,  t {{a,b},{c,d}}

A trivial problem: Diagonalize, calculate the power of the eigenvalues as a diagonal matrix, transform back using the eigensystem.

As a result, clear from the beginning: All matrix elements have the form

          a_ik Exp[t a_1 ] + b_ik Exp[t a_2]

Thats the kernel of any solution of time evolutions resulting from linear ODES.

You ran into the fact, that Wolfram Language's

Integrate, NIntegrate

with definite integrals is a definite no go.

Whenever you have a non-trivial definite Integrate, the system begins a never ending search on the hard disk you can follow in the Task Manager.

In contradistinction, indefinite Integrate is functional as long as standard solving strategies are working. But by any measure, the probability to write an indefinite integral yielding a closed finite algebraic expression as a result, is zero, unfortunately.

Except if the probability measure is defined on the set of integrals from exercises.

So as in this case.

Define your integrals as functions of the upper limit

PP[xx_, \[CapitalDelta]_, y_] := Integrate[MM[xx, \[CapitalDelta]], y]

Smat[xx_, \[CapitalDelta]_, y] := MatrixExp[PP[xx, \[CapitalDelta], y]];


  • 2
    $\begingroup$ In no way does this address the question that was posed. $\endgroup$ May 17, 2023 at 13:23
  • $\begingroup$ @Roland F Thank you for your response. I think I understand what you are trying to say but I'd like to clarify a few things. (1) The free parameter in my MM matrix would be xp, but the analogy to your t is not straightforward as it its not simply just a prefactor, it enters into an exponent in the matrix elements as well. (2) It is true that diagonalizing it might be easier, but it is unclear how to revert back to the eigen basis post integration, since PP is defined as the integral of MM (3) There is no concepts of probability here $\endgroup$
    – kowalski
    May 17, 2023 at 15:13
  • $\begingroup$ @Roland F Can you elaborate on what you might do to resolve the error/warnings in my case? $\endgroup$
    – kowalski
    May 17, 2023 at 15:14

I fixed the problems by inactivating Exp, because I know that Mathemtica does not Simplify

   MM[xp_, \[CapitalDelta]_] :=
     Activate[With[{Exp = Inactive[Exp]},
     (Chop[# , 10^-16] &)@
      ComplexExpand[Abs[Ap[xp]]^2] - 
       ComplexExpand[(Ap[xp])^2] + 
       Exp[(I*(ktht - kpt - kst[\[CapitalDelta]] - 
            kit[\[CapitalDelta]])*xp)]}, {-I*V2[\[CapitalDelta]]*
       ComplexExpand[(Ap[xp])^2] - 
       Exp[(-I*(ktht - kpt - kst[\[CapitalDelta]] - 
            kit[\[CapitalDelta]])*xp)], -I*V1[\[CapitalDelta]]*
        Abs[Ap[xp]]^2] - \[CapitalGamma]it[\[CapitalDelta]]/2}}, 
   Assumptions -> {NumericQ[xp] == True, xp >= 0, 
     xp \[Element] Reals}]]]];

Now there are complaints about my Assumptions, but in no time

         Smat[RandomReal[], RandomReal[]]

    {0.999991 + 0.000373264 I,   3.06243*10^-8 + 0.000258902 I},
 {-1.30875*10^-8 - 0.000191847 I,   0.999991 - 0.000536456 I}}

Its of course not sound, mathematically, insisting on xp to be numeric and the other conditons: it is a Pattern object. It is better style, to Simplify expressions first and then splice them into a function definition.

My standard in such cases is

  SetDelayed@@{ x_,  FullSimplify[expression]  }

  • $\begingroup$ Thank you for your response. Just a few things: (1) If I simply run the MM as defined in your post with Inactive[Exp], I get a warning saying FullSimplify: Warning: one or more assumptions evaluated to False. This is a warning that I get even when I ask to simply print PP[x,0], which is something that I didn't get in my definition of MM. (2) I cannot use RandomReal[] for the arguments of Smat, because well, they are not random numbers. I don't quite get the overall purpose of your approach. Could you elaborate? $\endgroup$
    – kowalski
    May 17, 2023 at 22:14
  • $\begingroup$ Design and conquer. First, write the mathematical expressions as functions with all constants as variables. Second, make all Simplify and complex expanding, TrigToExp and FullSimplify and define th result as a new function. Evaluate the functions with the constants and further FullSimplify. Check for occurencies of 0. and 1.*I etc and Rationalize eg by Rationalize[ expr, 2^-64] depending on the overall mean numeric entities. Finally never do a search for definite integrals. Indefinite integral evaluated at both boundaries are sound, and fast. If they dont deliver, resort to NIntegrate. $\endgroup$
    – Roland F
    May 18, 2023 at 7:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.