9
$\begingroup$

Bug introduced in 13.0 or earlier. Fixed in 13.1

I am trying to solve the following differential equation with Mathematica:

$$ z\,(1-z)\,\psi''(z)-2 z\,\psi'(z)-\frac{3}{4}\frac{\psi(z)}{z+a^2(1+z)}=0\,. $$

In order to do this, I input

DSolve[(1 - z) z Z''[z] - 2 z Derivative[1][Z][z] - 
   3/4 Z[z]/(z + a^2 (1 + z)) == 0, Z[z], z]

Much to my amusement, Mathematica returns

Z[z] -> Indeterminate

I would like to know why this is the case. Usually, when Mathematica cannot solve a differential equation, it returns the input, but in this case it returns the above. I am currently using Mathematica v13.0.0.

$\endgroup$

4 Answers 4

11
$\begingroup$

[I'm using V13.0.1]

There are three parts to this answer:

  1. What DSolve does in this case.
  2. How to transform the Heun equation to the OP's ODE.
  3. How to get the general solution to the OP's ODE.

First answer: What DSolve does

I can explain, literally, half of it:

ode = (1 - z) z Z''[z] - 2 z Derivative[1][Z][z] - 
    3/4 Z[z]/(z + a^2 (1 + z)) == 0;

coeff = Coefficient[-Z''[z] /. First@Solve[ode, Z''[z]], {Z'[z], 
    Z[z]}] // Simplify
(*  {2/(-1 + z), 3/(4 (-1 + z) z (z + a^2 (1 + z)))}  *)

DSolve`DSolveSecondOrderODE; (* autoload 2nd order fns *)
basis = DSolve`DSolveSecondOrderODEDump`DSolveSpecialOrder2Linear[
  coeff, z]
(*
  {Indeterminate, 
   z HeunG[-(a^2/(1 + a^2)), 
     2 - 3/(4 (1 + a^2)) + 2 (-1 - a^2/(1 + a^2)), 2, 1, 2, 2, z]}
*)

basis . {C[1], C[2]}
(*  Indeterminate  *)

ode /. Z -> 
  Function[z, 
   z HeunG[-(a^2/(1 + a^2)), 
     2 - 3/(4 (1 + a^2)) + 2 (-1 - a^2/(1 + a^2)), 2, 1, 2, 2, z] 
   ] // FullSimplify
(*  True  *)

So what's going on? WRI improved DSolve in V13 (see New in 13: Symbolic & Numeric Computation), and I believe that's what we're seeing. The particular solver finds half the basis for the solution space, but the other basis function comes out Indeterminate. This usually happens if something like 0/0 happens, which could happen because things weren't simplified carefully. (If I have time, I'll see if I can track it down.) We can see from checking with FullSimplify that the second basis function satisfies the ode.

We can also see that the standard linear combination of the basis functions results in Indeterminate (basis . {C[1], C[2]}).

Consider reporting it as a bug. Since this is new functionality, they probably have not run into this edge case. It looks like it could be handled better.

How to transform the Heun equation to the OP's ODE.

@josh's interest prompted me to share the basis for one of my comments below. It is not much help to me, but maybe it will be for others. It shows only the relationship between the OP's equation and the Heun equation.

odeHeun = (-q + z α β) Z[z] +
  ((-1 + z) (-a + z) γ +
      (-1 + z) z (1 + α + β - γ - δ) + 
      z (-a + z) δ) Z'[z] +
  (-1 + z) z (-a + z) Z''[z] == 0;

(* transform Heun to OP's *)
xf = Block[{HeunG = List},
  Thread[
   HeunG[a, q, α, β, γ, δ, z] -> 
    HeunG[-(a^2/(1 + a^2)), -(3/(4 (1 + a^2))), 0, 1, 0, 2, z]]
  ]
(*
{a -> -(a^2/(1 + a^2)), 
 q -> -(3/(4 (1 + a^2))), α -> 0, β -> 1, γ -> 0, δ -> 2,
 z -> z}
*)

(* Check xf transforms the Heun equation to the OP's *)
Solve[odeHeun /. xf, Z''[z]] // FullSimplify
Solve[ode, Z''[z]] // FullSimplify
% === %%
(*  True  *)

(* Extract the general basis and apply xf to get OP's case *)
First@DSolve[odeHeun, Z[z], z] /. (Thread[{C[1], C[2]} -> #] & /@ 
    IdentityMatrix[2]) /. HeunG -> Inactive[HeunG]
% /. xf
% // Activate

Mathematica graphics

We have γ = 0 in solution that evaluates to Indeterminate. The series expansion of the general HeunG function shows this and γ a negative integer are exceptional cases (as implied in the DLMF):

Series[HeunG[a, q, α, β, γ, δ, z], {z, 0, 3}]

Mathematica graphics

Possibly, the current state of knowledge means there is no way to express the complete solution, but I'm not an expert on Heun functions. Maybe the best place to seek an answer to the question whether there is a complete solution is on one of the other mathematics sites. See below.

Update 2: How to get the general solution

At the singular point $z=0$, the OP's ODE has characteristic exponents (or indicial roots) $0$ and $1$. We have one solution above $Z_1$ that solution corresponds to a series solution expanded about the singular point $z=0$ with characteristic exponent $1$. The difference in the characteristic exponents is an integer, so there is a Frobenius solution of the form $$ Z_2 = \sum_{k=0}^\infty b_k z^k + a \log z \cdot Z_1 \,.$$ Whether this may be expressed in terms of higher transcendental functions is unclear to me. (I'd guess not, but I don't know.) However, AsymptoticDSolveValue may be used to get a finite expansion of the Frobenius series. It takes a long time to get more than a few terms, though.

Since the ODE is linear, we may also try to solve it by reduction of order, setting $Z_2(z)=Y(z)Z_1(z)$. This produces a general solution in terms of an unevaluated integral (which is the sole basis for guessing that a solution in terms of higher transcendental functions is not possible).

(* z1sol from spelunking DSolve: *)
z1sol = Z -> 
   Function[z, 
    z HeunG[-(a^2/(1 + a^2)), 
      2 - 3/(4 (1 + a^2)) + 2 (-1 - a^2/(1 + a^2)), 2, 1, 2, 2, z]];

(* DSolve will do the last step in the reduction of order
 * for us, so we need only make the substitution  *)
odeY = odeHeun /. xf /. Z -> Function[z, Y[z] Z[z]] /. z1sol // 
   FullSimplify;
ysol = First@DSolve[{odeY,
     (* irritating way to get the initial point for the
      * integral different than 1, which is a singular point:
      *)
     dummy'[z] == 0, dummy[1/2] == 0},
    {Y, dummy}, z];
gensol = Z -> Function @@ {z,
    (Y[z] + C[2] /. ysol) (Z[z] /. z1sol) // FullSimplify}

(* You can omit the dummy variable and edit the integral 
 * with the following substitution Rule:  *)
(*bettergensol = gensol //. Inactive[Integrate][f_,{x_,1,b_}] :>
    RuleCondition[Inactive[Integrate][f,{x,1/2,b}],True]*)

(* Check: *)
odeHeun /. xf /. gensol // FullSimplify

Mathematica graphics

On the dummy variable. Adding the dummy equation dummy'[z] == 0 makes the system one dimension larger. This can have consequences, sometimes ones you would wish to avoid like DSolve thinking the system is too complicated to try to solve (a problem in previous version, which I assume is still present). However, starting the integral at a point where HeunG is undefined (or Indeterminate in this case) is rather unhelpful. The replacement Rule method for bettergensol is probably the better way to handle the integral. The dummy method was to show that there is a way to get DSolve to do it for you. There does not seem to be any other way to tell DSolve where to start the integral.

Summary/synopsis

  1. This is a bug. The fact that gensol is easily obtained by standard methods shows that DSolve could be improved in its handling of the Heun equation. Even if that were not possible, the Indeterminate result for one of the basis functions should be handled better.

  2. One may obtain the Frobenius series with order = 3; AsymptoticDSolveValue[ode, Z, {z, 0, order}], if that is desired.

  3. HeunG is the solution to the Heun equation corresponding to the characteristic exponent 0 at the singular point z == 0 for which HeunG[..., 0] = 1. It is denoted $H\!\ell$ in DLMF; the solution in z1sol corresponds to eq. (31.3.5). (I'm not an expert, so let me know if I made any blunders.)

  4. HeunG[-(1/2), -(11/8), 2, 1, 2, 2, 2.] is another bug, which I will report to WRI. It is obtained by setting a -> 1, z -> 2. in the HeunG factor in z1sol. It happens for most settings of a as well. Note 0, 1, and the first argument of HeunG are singular points of the Heun equation. When the first argument is equal to 0 or 1 (confluent case) or z is equal to one of singular points, sometimes HeunG does not evaluate or returns Indeterminate.

$\endgroup$
6
  • $\begingroup$ HeunG[-(a^2/(1 + a^2)), -(3/(4 (1 + a^2))), 0, 1, 0, 2, z] appears to be the first solution-basis function, but it evaluates to Indeterminate. $\endgroup$
    – Michael E2
    Apr 10, 2022 at 0:26
  • $\begingroup$ I don't have access to it, but would be interesting if we could compare to Maple. I will see if I can get Maple if nobody adds anything $\endgroup$
    – bmf
    Apr 10, 2022 at 0:29
  • 1
    $\begingroup$ I can only add that the HeunG[] expression in my previous comment is a formal solution of a Heun equation that is equivalent to the OP's ODE. I cannot explain why it evaluates to Indeterminate. $\endgroup$
    – Michael E2
    Apr 10, 2022 at 2:49
  • 1
    $\begingroup$ The situation seems like a bug. As you suggested it should be reported $\endgroup$
    – bmf
    Apr 10, 2022 at 3:04
  • $\begingroup$ @Michael E2: Do I need to load a a package for basis = DSolveDSolveSecondOrderODEDumpDSolveSpecialOrder2Linear[ coeff, z] as Mathematica returns unevaluated on my machine for V. 13. $\endgroup$
    – josh
    Apr 10, 2022 at 13:44
5
$\begingroup$

I don't have access to it, but would be interesting if we could compare to Maple. I will see if I can get Maple if nobody adds anything

Maple 2022 can't solve it. It gives solution in terms of DESol

DESol is a data structure to represent the solution of a differential equation. It is to dsolve as RootOf is to solve.

Here is its attempt

    ode:=z*(1-z)*diff(y(z),z$2)-2*z*diff(y(z),z)-3/4*y(z)/(z+a^2*(1+z))=0
    infolevel[dsolve]:=5;
    dsolve(ode)

It gives

Methods for second order ODEs:
--- Trying classification methods ---
trying a quadrature
checking if the LODE has constant coefficients
checking if the LODE is of Euler type
trying a symmetry of the form [xi=0, eta=F(x)]
   testing BRANCH 1 ->
   testing BRANCH 2 ->
   testing BRANCH 3 ->
   testing BRANCH 4 ->
   testing BRANCH 5 ->
   testing BRANCH 6 ->
checking if the LODE is missing 'y' 
-> Trying a Liouvillian solution using Kovacic's algorithm
<- No Liouvillian solutions exists
-> Trying a solution in terms of special functions:
   -> Bessel
   -> elliptic
   -> Legendre
   -> Whittaker
      -> hyper3: Equivalence to 1F1 under 'a power @ Moebius'
   -> hypergeometric
      -> heuristic approach
      -> hyper3: Equivalence to 2F1, 1F1 or 0F1 under 'a power @ Moebius'
         -> trying extended approach, based on solving the ODE satisfied by the "derivative" of y(z)
            -> ODE succesfully computed; trying to solve: diff(diff(y(z),z),z) = -(5*a^2*z^2+2*a^2*z-a^2+5*z^2-2*z)/(a^2*z+a^2+z)/z/(-1+z)*diff(y(z),z)-1/4*y(z)*(16*a^2*z^2-8*a^2*z-8*a^2+16*z^2-13*z-3)/(a^2*z+a^2+z)/z/(-1+z)^2
         <- extended approach, failed in solving the ODE for the "derivative" of y(z)
   -> Mathieu
      -> Equivalence to the rational form of Mathieu ODE under 'a power @ Moebius'
         -> trying extended approach, based on solving the ODE satisfied by the "derivative" of y(z)
            -> ODE succesfully computed; trying to solve: diff(diff(y(z),z),z) = -(5*a^2*z^2+2*a^2*z-a^2+5*z^2-2*z)/(a^2*z+a^2+z)/z/(-1+z)*diff(y(z),z)-1/4*y(z)*(16*a^2*z^2-8*a^2*z-8*a^2+16*z^2-13*z-3)/(a^2*z+a^2+z)/z/(-1+z)^2
         <- extended approach, failed in solving the ODE for the "derivative" of y(z)
trying a solution in terms of MeijerG functions
-> Heun: Equivalence to the GHE or one of its 4 confluent cases under 'a power @ Moebius'
-> trying a solution of the form r0(x) * Y + r1(x) * Y' where Y = exp(int(r(x), dx)) * 2F1([a1, a2], [b1], f)
-> Tackling the "Normal Form" of the given linear ODE:
   trying a symmetry of the form [xi=0, eta=F(x)]
      testing BRANCH 1 ->
      testing BRANCH 2 ->
      testing BRANCH 3 ->
      testing BRANCH 5 ->
      testing BRANCH 6 ->
-> Tackling the linear ODE using alternative approaches:
   trying differential order: 2; exact nonlinear
   trying symmetries linear in x and y(x)
   trying to convert to a linear ODE with constant coefficients
      testing BRANCH 1 ->
   trying 2nd order, integrating factor of the form mu(x,y)
   Trying an integrating factor as the solution to the adjoint linear ODE
   -> Trying a solution in terms of special functions:
      -> Bessel
      -> elliptic
      -> Legendre
      -> Whittaker
         -> hyper3: Equivalence to 1F1 under 'a power @ Moebius'
      -> hypergeometric
         -> heuristic approach
         -> hyper3: Equivalence to 2F1, 1F1 or 0F1 under 'a power @ Moebius'
            -> trying extended approach, based on solving the ODE satisfied by the "derivative" of y(z)
               -> ODE succesfully computed; trying to solve: diff(diff(y(z),z),z) = -(3*a^2*z^3-9*a^2*z^2+3*a^2*z+3*z^3+3*a^2-9*z^2+6*z)/(a^2*z+a^2+z)/z/(-1+z)/(8*a^2*z^2+8*a^2*z+8*z^2+3*z-3)*diff(y(z),z)-1/4*y(z)*(72*a^2*z^2-48*a^2*z-24*a^2+57*z^2-66*z+9)/(a^2*z+a^2+z)/z/(-1+z)^2/(8*a^2*z^2+8*a^2*z+8*z^2+3*z-3)
            <- extended approach, failed in solving the ODE for the "derivative" of y(z)
      -> Mathieu
         -> Equivalence to the rational form of Mathieu ODE under 'a power @ Moebius'
            -> trying extended approach, based on solving the ODE satisfied by the "derivative" of y(z)
               -> ODE succesfully computed; trying to solve: diff(diff(y(z),z),z) = -(3*a^2*z^3-9*a^2*z^2+3*a^2*z+3*z^3+3*a^2-9*z^2+6*z)/(a^2*z+a^2+z)/z/(-1+z)/(8*a^2*z^2+8*a^2*z+8*z^2+3*z-3)*diff(y(z),z)-1/4*y(z)*(72*a^2*z^2-48*a^2*z-24*a^2+57*z^2-66*z+9)/(a^2*z+a^2+z)/z/(-1+z)^2/(8*a^2*z^2+8*a^2*z+8*z^2+3*z-3)
            <- extended approach, failed in solving the ODE for the "derivative" of y(z)
   -> Tackling the linear ODE using alternative approaches:
      trying 2nd order exact linear
      trying symmetries linear in x and y(x)
      trying to convert to a linear ODE with constant coefficients
         testing BRANCH 1 ->
      trying to convert to an ODE of Bessel type
   trying to convert to an ODE of Bessel type
   -> trying reduction of order to Riccati
      Chini's absolute invariant is: (1/3)*(a^2*z^2+4*a^2*z+a^2+z^2+2*z)^2/(z*(a^2*z+a^2+z)*(-1+z))
      trying Riccati sub-methods:
         -> trying a symmetry pattern of the form [F(x)*G(y), 0]
         -> trying a symmetry pattern of the form [0, F(x)*G(y)]
         -> trying a symmetry pattern of the form [F(x),G(x)*y+H(x)]
--- Trying Lie symmetry methods, 2nd order ---
 -> Computing symmetries using: way = 3
                             [0, y]

 <- successful computation of symmetries.
 -> Computing symmetries using: way = 5

enter image description here

$\endgroup$
4
  • $\begingroup$ Many thanks for taking the time to check this! $\endgroup$
    – bmf
    Apr 10, 2022 at 2:35
  • $\begingroup$ Do you know if it's similar to DifferentialRoot[Function[{Z, z}, {(1 - z) z Z''[z] - 2 z Derivative[1][Z][z] - 3/4 Z[z]/(z + a^2 (1 + z)) == 0, Z[0] == Z0, Z'[0] == Z1}]]? $\endgroup$
    – Michael E2
    Apr 10, 2022 at 2:41
  • $\begingroup$ @MichaelE2 I think DifferentialRoot in Mathematica is like ODESolStruc in Maple but I never was able to confirm this. Maple says For high order ODEs it may happen that dsolve succeeds in reducing the order of the ODE but not in solving the problem to the end. In those cases, the answer is expressed using ODESolStruc. The rule of thumb I use is that if Maple returns ODESolStruc or DESol then the ODE is not solved. And if Mathematica returns DifferentialRoot then also the ode is not solved. $\endgroup$
    – Nasser
    Apr 10, 2022 at 2:55
  • $\begingroup$ For DESol, I read in Maple's docs, "DESol...is to dsolve as RootOf is to solve." I think of DifferentialRoot is to DSolve as Root is to Solve, although each has its limitations. DifferentialRoot can compute numerical values of a solution and series expansions (when the ODE has polynomial coefficients and parameters are sufficiently defined). It has trouble near singular points, so my defining it in my comment at z == 0 (e.g. IC Z[0] == Z0) is a mistake, I think. $\endgroup$
    – Michael E2
    Apr 10, 2022 at 3:13
3
$\begingroup$

An extended comment:

My understanding is this DE can be transformed into the standard Heung DE by a suitable change of variables. DSolve can solve the standard form:

sol = DSolve[ 
  y''[z] + (\[Gamma]/z + \[Delta]/(z - 1) + (
       1 + \[Alpha] + \[Beta] - \[Gamma] - \[Delta])/(z - a)) y'[
      z] + (\[Alpha] \[Beta] z - q)/(z (z - 1) (z - a)) y[z] == 0, 
  y[z], z]

producing two solutions:

(* {{y[z] -> 
   C[1] HeunG[a, q, \[Alpha], \[Beta], \[Gamma], \[Delta], z] + 
    z^(1 - \[Gamma]) C[
     2] HeunG[a, 
      q - (-1 + \[Gamma]) (1 + \[Alpha] + \[Beta] - \[Gamma] + (-1 + 
             a) \[Delta]), 1 + \[Beta] - \[Gamma], 
      1 + \[Alpha] - \[Gamma], 2 - \[Gamma], \[Delta], z]}} *)

I think if we either first transformed the OP's DE into standard form then execute DSolve, or even better, use the two solutions returned by DSolve for the standard form solution and then adjust the constants appropriately to match the OPs original DE. However, I'm not initially able to make this transformation. This might be a work-around to obtaining the solutions.

Edit:

Using Michael's transformation code above exactly as written, running on a clean kernel, I do not obtain an Indeterminate for the final hDev answer but rather 2 HeunG functions:

odeHeun = (-q + z \[Alpha] \[Beta]) Z[
      z] + ((-1 + z) (-a + z) \[Gamma] + (-1 + 
          z) z (1 + \[Alpha] + \[Beta] - \[Gamma] - \[Delta]) + 
       z (-a + z) \[Delta]) Z'[z] + (-1 + z) z (-a + z) Z''[z] == 0;
ode = (1 - z) z Z''[z] - 2 z Derivative[1][Z][z] - 
    3/4 Z[z]/(z + a^2 (1 + z)) == 0;

(*transform Heun to OP's*)
xf = Block[{HeunG = List}, 
  Thread[HeunG[a, q, \[Alpha], \[Beta], \[Gamma], \[Delta], z] -> 
    HeunG[-(a^2/(1 + a^2)), -(3/(4 (1 + a^2))), 0, 1, 0, 2, z]]]
(*{a->-(a^2/(1+a^2)),q->-(3/(4 (1+a^2))),\[Alpha]->0,\[Beta]->1,\
\[Gamma]->0,\[Delta]->2,z->z}*)

(*Check xf transforms the Heun equation to the OP's*)
sVal1 = Solve[odeHeun /. xf, Z''[z]] // FullSimplify
sVal2 = Solve[ode, Z''[z]] // FullSimplify
sVal1 === sVal2
(*True*)

(*Extract the general basis and apply xf to get OP's case*)
hDev = First@
    DSolve[odeHeun, Z[z], z] /. (Thread[{C[1], C[2]} -> #] & /@ 
     IdentityMatrix[2]) /. HeunG -> Inactive[HeunG]
hDev /. xf
hDev // Activate

(* {{Z[z] -> 
       HeunG[a, q, \[Alpha], \[Beta], \[Gamma], \[Delta], z]}, {Z[z] -> 
       z^(1 - \[Gamma])
         HeunG[a, 
         q - (-1 + \[Gamma]) (1 + \[Alpha] + \[Beta] - \[Gamma] + (-1 + 
                a) \[Delta]), 1 + \[Beta] - \[Gamma], 
         1 + \[Alpha] - \[Gamma], 2 - \[Gamma], \[Delta], z]}} *)
$\endgroup$
7
  • 1
    $\begingroup$ I can do this transformation based on my spelunking; there is, I think, an algorithm, but I'm not famliar with it. You end up with the Indeterminate Heun function, so I don't see that it really helped. I left a comment in the hopes that something else to try might suggest itself. Series[HeunG[a, q, α, β, γ, δ, z], {z, 0, 2}] suggests why γ=0 may be a problem. $\endgroup$
    – Michael E2
    Apr 10, 2022 at 14:11
  • 1
    $\begingroup$ I added the transformation to my answer. Take a look if you're interested. $\endgroup$
    – Michael E2
    Apr 10, 2022 at 14:57
  • $\begingroup$ @Michael E2: I ran your transformation code on my machine under ver 13 exactly as written and obtain two Heung functions for the final answer and no Indeterminate although it's a bit tough to follow your code and maybe I'm not understanding something. $\endgroup$
    – josh
    Apr 10, 2022 at 17:57
  • $\begingroup$ You need to plug in the parameter values in the expression at the end of your post. Apply ... /. Block[{HeunG = List}, Thread[HeunG[a, q, \[Alpha], \[Beta], \[Gamma], \[Delta], z] -> HeunG[-(a^2/(1 + a^2)), -(3/(4 (1 + a^2))), 0, 1, 0, 2, z]] ] -- Then you get Indeterminate. I gotta go somewhere, so I'll double-check later. $\endgroup$
    – Michael E2
    Apr 10, 2022 at 18:00
  • $\begingroup$ I though the hDev/.xf construct did this. I'll work on it. Thank! You guys do brilliant work here with DEs I think. $\endgroup$
    – josh
    Apr 10, 2022 at 18:05
3
$\begingroup$

Fixed in 13.1. The following is screen shot

enter image description here

$\endgroup$

Your Answer

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

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