Why does the Part::partw error appear when using DSolve?

Question

Bug introduced in 9.0 or earlier and fixed in 11.1.0

---

I solved the following ordinary differential equation for $F(x)$ using DSolve in Mathematica

$$ \begin{align} \begin{split} & A x^2 (x^2 + 1) F'''(x) - \left(3 A (b - 2) x^3 + a x^2 +A (b - 2)x +a \right)F''(x) \\ & + 3 (b - 1) x \left(A (b - 2) x + \frac{2}{3} a \right)F'(x) -(b - 1) b \left(A(b - 2) x + a\right) F(x) = 0 \end{split} \end {align} $$ Mathematica Input

DSolve[A x^2 (x^2 + 1) F'''[x] - (3 A (b - 2) x^3 + a x^2 + A (b - 2) x + a) F''[x] + 3 (b - 1) x (A (b - 2) x + 2/3 a) F'[x] - (b - 1) b (A (b - 2) x + a) F[x] == 0, F, x];

After some time, Mathematica actually produced the closed form solution below. $$ \begin{align} \begin{split} F(x) = & C_1 ( x + i)^b + C_2 ( x - i) ^b + C_3\left[a \left(b-1+x \left(1+\sqrt {-(b - 1 )^2}-b \right) \right. \right. \\ & \left. \left. +\sqrt {-(b - 1)^2} \right)+A (b - 2) \left(b-1+x \sqrt {-(b - 1)^2} \right) \right]^b \end{split} \end{align} $$

However, it also gave me an error message

Part::partw: Part 1 of {} does not exist. >>

From another answer, I found out that

Part::partw - This error occurs when you're trying to obtain a part of an expression at a position longer than the Length of the expression (or subexpression).

Questions

• Why does the error message come up in DSolve?
• Is the above solution a valid general solution?
• What algorithm might Mathematica use for this type of ODE?

Any thoughts are much appreciated!

Answer 1

The computation

eq = A x^2 (x^2 + 1) F'''[x] - (3 A (b - 2) x^3 + a x^2 + A (b - 2) x + a) F''[x] + 
    3 (b - 1) x (A (b - 2) x + 2/3 a) F'[x] - (b - 1) b (A (b - 2) x + a) F[x];
s = DSolve[eq == 0, F, x][[1, 1]];

after further simplification yields

s[[2, 2]] = FullSimplify[s[[2, 2]]]; s
(* F -> Function[{x}, (-I + x)^b C[1] + (I + x)^b C[2] + 
   (A (-2 + b) (-1 + b + Sqrt[-(-1 + b)^2] x) + a (-1 + Sqrt[-(-1 + b)^2] + 
   b + x + Sqrt[-(-1 + b)^2] x - b x))^b C[3]] *)

which is the solution given in the question. (The error message in the question is, of course, also produced.) Indeed, it is a valid solution,

Simplify[eq /. s]
(* 0 *)

However, it is not the most general solution, as can be seen by forming the Wronskian of the three solutions contained in F.

Wronskian[List @@ (F /. s)[x], x]
(* 0 *)

In other words, the term multiplying C[3] is a linear combination of the terms multiplying C[1] and C[2]. So there must be a third independent solution. In principle, it can be obtained using the procedure described here. Begin by defining

w = Wronskian[eq == 0, F, x] // FullSimplify
(* E^(-(a/(A x))) x^(-2 + b) (1 + x^2)^(-2 + b) *)
y1 = (x - I)^b; y2 = (x + I)^b;
w12 = Wronskian[{y1, y2}, x] // FullSimplify
(* -2 I b (-I + x)^(-1 + b) (I + x)^(-1 + b) *)

Then, the third solution is given by,

Integrate[Simplify[y2 (y1 w/w12^2 /. x -> xp) - y1 (y2 w/w12^2 /. x -> xp)], xp]

Unfortunately, this command returns unevaluated. Alternatively, the third solution, here designated y3[x], can be obtained by solving

DSolve[(w12 y3''[x] - D[w12, x] y3'[x] + 
    (D[y1, x] D[y2, {x, 2}] - D[y2, x] D[y1, {x, 2}]) y3[x]) == w, y3, x]

This computation ran for about one half hour, at which point the kernel crashed. It appears, therefore, that the third solution cannot be obtained in symbolic form.

To address the specific issues posed in the question,

1. Why the error message? A bug in DSolve.
2. Is the solution in the question a valid, general solution. Valid, but not general.
3. How was the solution obtained? Unclear (to me) how y1 was obtained. However, once obtained, y1 must be the complex conjugate of y2. The third term is not linearly independent and, in this context, must be considered an error.

By the way, an alternative approach,

SetOptions[Solve, Method -> Reduce];
r = DSolve[eq == 0, F, x][[1, 1]]
(* F -> Function[{x}, (-I + x)^b C[1] + (I + x)^b C[2] + 
   (-1 + Sqrt[-(-1 + b)^2] - (1 + Sqrt[-(-1 + b)^2]) x + b (1 + x))^b C[3]] *)

produces a different solution, also valid but not general, along with the same error message.

Version 11.1 update

Version 11.1 returns three independent solutions, although in terms of DifferentialRoot and no error message.

SetSystemOptions["HolonomicOptions" -> {"VariableSymbol" -> x, "FunctionSymbol" -> y}];
s = DSolve[eq == 0, F, x][[1, 1]]
(* F -> DifferentialRoot[Function[{y, x}, 
   {-((-1 + b)*b*(a - 2*A*x + A*b*x)*y[x]) + (-1 + b)*x*(2*a - 6*A*x + 3*A*b*x)*y'[x] + 
    (-a + 2*A*x - A*b*x - a*x^2 + 6*A*x^3 - 3*A*b*x^3)*y''[x] + A*x^2*(1 + x^2)*y'''[x] 
    == 0, y[1] == C[1], y'[1] == C[2], y'[1] == C[3]}]] *)

(The first line of code provides a cleaner format, as described in the answer to question 123540.) Because verifying the solution with, for instance,

FullSimplify[eq /. s]

runs seemingly forever, the following can be used to lend credence to the result.

rule = Thread[{C[1], C[2], C[3], A, a, b} -> Round[100 RandomReal[{1, 5}, 6]]/100]
sn = First@NDSolve[{eq == 0, F[1] == C[1], F'[1] == C[2], F''[1] == C[3]} /. 
    rule, F, {x, 1, 5}];
f = (F /. s /. rule);
Plot[{(F /. sn)[x], f[x]}, {x, 1, 5}]
(* {C[1] -> 84/25, C[2] -> 417/100, C[3] -> 401/100, 
    A -> 92/25, a -> 237/100, b -> 223/100} *)