# Mathematica not able to confirm its own solution to differential equation

I type the following into Mathematica:

DSolve[q''[x] + 2 x/(x^2 - 1) q'[x] - 4*q[x]/(x^2 - 1) == 0, q[x], x]


It gives me the result

q[x] -> C[1] LegendreP[1/2 (-1 + Sqrt[17]), x] +
C[2] LegendreQ[1/2 (-1 + Sqrt[17]), x]


Suppose I take C[2]=0 and C[1] =1, and then I try the following command:

FullSimplify[
D[D[LegendreP[1/2 (-1 + Sqrt[17]), x], x], x] +
2*x/(x^2 - 1)*D[LegendreP[1/2 (-1 + Sqrt[17]), x], x] -
4*LegendreP[1/2 (-1 + Sqrt[17]), x]/(x^2 - 1)]


This is evaluating the LHS of the diff eq at the solution, as provided by Mathematica, and yet the result I get is NOT zero:

I get

(1/(2 (-1 +
x^2)^2))((9 + Sqrt[17]) LegendreP[1/2 (-1 + Sqrt[17]),
x] - (19 + 3 Sqrt[17]) x LegendreP[1/2 (1 + Sqrt[17]), x] +
2 (5 + Sqrt[17]) LegendreP[1/2 (3 + Sqrt[17]), x])


I was wondering if someone knew how I could make it fully simplify the result to give me zero, as numerically evaluating it gives extremely small numbers.

-

As the other answers have mentioned, you can confirm the solution without specifying values for C[1] and C[2]

lhs = q''[x] + 2 x/(x^2 - 1) q'[x] - 4*q[x]/(x^2 - 1);
sol = DSolve[lhs == 0, q, x][[1, 1]];
FullSimplify[lhs /. sol]


0

When C[2] is set to zero, FullSimplify is unable to find the right transformations to reduce the expression to zero, and the simplest form it can find is, as noted in the question:

result = FullSimplify[lhs /. sol /. {C[1] -> 1, C[2] -> 0}]


$\frac{1}{2 \left(-1+x^2\right)^2}\left(\left(9+\sqrt{17}\right) \text{LegendreP}\left[\frac{1}{2} \left(-1+\sqrt{17}\right),x\right]-\left(19+3 \sqrt{17}\right) x \text{LegendreP}\left[\frac{1}{2} \left(1+\sqrt{17}\right),x\right]+2 \left(5+\sqrt{17}\right) \text{LegendreP}\left[\frac{1}{2} \left(3+\sqrt{17}\right),x\right]\right)$

This result is zero, which can be confirmed by evaluating for numerical values of $x$. Be aware though, that using machine precision may lead to numerical errors:

result /. x -> 0.1


1.81242*10^-15

A good idea when you see unexpected results like this is to make Mathematica use arbitrary-precision by specifying the precision of the number:

result /. x -> 0.120


0.*10^-18

-
thanks, it worked. I would have thought it was easier for mathematica to simplify something for one particular value of parameter, rather than doing it for all possible values!, but I guess not. – user7173 Apr 29 '13 at 17:29

You can do :

sol[x_] = q[x] /.
First@DSolve[q''[x] + 2 x/(x^2 - 1) q'[x] - 4*q[x]/(x^2 - 1) == 0, q[x], x]

FullSimplify[Derivative[2][sol][x] + 2 x/(x^2 - 1) Derivative[1][sol][x] -
4*sol[x]/(x^2 - 1)]
(* 0 *)

-

On the other hand, we can use DifferentialRootReduce[] on LegendreP[1/2 (-1 + Sqrt[17]), x] to see what linear differential equation is satisfied by it:

Operate[FullSimplify, DifferentialRootReduce[LegendreP[1/2 (-1 + Sqrt[17]), x], x]]
DifferentialRoot[Function[{y, x}, {-4 y[x] + 2 x y'[x] + (-1 + x^2) y''[x] == 0,
y[0] == Sqrt[π]/(Gamma[3/4 - Sqrt[17]/4] Gamma[(3 + Sqrt[17])/4]),
y'[0] == -((1 + Sqrt[17]) Sqrt[π])/
(2 Gamma[1/4 - Sqrt[17]/4] Gamma[(5 + Sqrt[17])/4])}]][x]


where all the formal symbols were replaced with their normal equivalents for clarity. (Note the use of Operate[] so that FullSimplify[] is applied only to the DifferentialRoot[] object.) In any case, the ODE thus produced (the first element in the list within the Function[]) can then be seen to be the same as the one in the OP.

Apparently, the reason why Mathematica could not recognize the complicated expression in Simon's answer as $0$ is due to its inability to see that one can use the recurrence relation for the Legendre functions in this case.

expr = FullSimplify[q''[x] + 2 x/(x^2 - 1) q'[x] - 4*q[x]/(x^2 - 1) /.
q -> Function[x, LegendreP[1/2 (-1 + Sqrt[17]), x]]];

FullSimplify[expr /. Activate[Function[{n, z}, Inactivate[LegendreP[n, z] ==
(((2 n - 1) z) LegendreP[n - 1, z] -
(n - 1) LegendreP[n - 2, z])/n]][(3 + Sqrt[17])/2, x] /. Equal -> Rule]]
0


where I took the required form of the recurrence from the results of MathematicalFunctionData[LegendreP, "DifferenceEquations"].

-

I think you may not chose one parameter to 0. E.g. try c2=2 and c1=1:

FullSimplify[
D[LegendreP[1/2 (-1 + Sqrt[17]), x] +
2 LegendreQ[1/2 (-1 + Sqrt[17]), x], {x, 2}] +
2*x/(x^2 - 1)*
D[LegendreP[1/2 (-1 + Sqrt[17]), x] +
2 LegendreQ[1/2 (-1 + Sqrt[17]), x], x] -
4*(LegendreP[1/2 (-1 + Sqrt[17]), x] +
2 LegendreQ[1/2 (-1 + Sqrt[17]), x])/(x^2 - 1)]


As output you will get 0. But I admit that this is really odd. Actually I agree with you and would say that you may set one of your constants to 0. Also, if I replace one of the constants with non-integers, I get back values which are only almost 0.

But if you try

 DSolve[D[q[x], {x, 2}] + 2 x/(x^2 - 1) *D[q[x], x] -  4*q[x]/(x^2 - 1) == 0, q, x]

FullSimplify[D[q[x], {x, 2}] + 2 x/(x^2 - 1)*D[q[x], x] - 4*q[x]/(x^2 - 1) /. %]
`

your result is confirmed. So, I guess if you work with real numbers, you get numerical problems.

-