5
$\begingroup$

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!

$\endgroup$
5
  • $\begingroup$ Please take a look at this page: How to copy code from Mathematica so it looks good on this site and re-paste your code for the differential equation. All those \[Prime] are an artifact and should be fixed. It will be a lot easier for people to help you if they can run your code. $\endgroup$
    – MarcoB
    Jun 7, 2016 at 21:12
  • $\begingroup$ @MarcoB: Thanks for the hint. Typed up the input now, as I currently only have the Wolfram CDF Player available at home. Hopefully it should be working now. $\endgroup$
    – TimD
    Jun 7, 2016 at 21:33
  • $\begingroup$ Tim, thanks for the update. I ran your code and observed the same Part error. I am not sure where it comes from, except that it seems to come from inside DSolve. I'm afraid that I am not sure about points 2. and 3. $\endgroup$
    – MarcoB
    Jun 7, 2016 at 23:07
  • 2
    $\begingroup$ @TimD You might want to report this to Wolfram Support; perhaps it's a symptom of an internal problem within DSolve. $\endgroup$
    – MarcoB
    Jun 9, 2016 at 18:23
  • $\begingroup$ Can somebody with access to multiple versions please add the customary header? $\endgroup$ Jun 10, 2016 at 4:46

1 Answer 1

5
$\begingroup$

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} *)

enter image description here

$\endgroup$
3
  • $\begingroup$ Put another way: Mathematica has yet again returned a generic solution. $\endgroup$ Jun 10, 2016 at 4:36
  • $\begingroup$ Thanks for the well-written answer! I have reported the bug to the Wolfram Support. One more question regarding your answer: Can we conclude the linear dependence for any value of the parameters a, b and A from just knowing that Wronskian[List @@ (F /. s)[x], x] is zero? This article suggests that there are additional conditons. $\endgroup$
    – TimD
    Jun 10, 2016 at 6:11
  • $\begingroup$ @TimD Note footnote 1 (on page 2) of the reference. W = = 0 is sufficient for linear independence, if the functions are analytic. The three functions in the solution (erroneously) supplied by DSolve are analytic except, perhaps, at a few isolated points. Moreover, a bit of algebra is sufficient to show that the third solution is proportional to either the first or second solution, depending on the value of b. So, the three functions are, indeed, dependent. $\endgroup$
    – bbgodfrey
    Jun 10, 2016 at 21:38

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.