# how to simplify large expression with lots of special functions in it (BesselY, Hypergeometric, MeijerG etc…)

I saw this DE in Maple forum. When solving it using Mathematica 9.01, even though the result was correct (both solutions gave the same numerical answer for some random values), Mathematica's answer was much larger in size.

Using Simplify or FullSimplify did not help much (expression was still order of magnitude larger than Maple's result). I think special simplification rules might be needed?

I'd like to ask if there are other special simplifications to use, or commands I might be missing to reduce this result more. Below is the ODE and Maple and Mathematica results.

Mathematica:

ClearAll[y, x];
ode = y'''[x] + (1/x) y'[x] - (1/x^2) y[x] == 0;
ic = {y[1] == 1, y'[1] == 0, y''[1] == 1};
sol = y[x] /. First@DSolve[{ode, ic}, y[x], x];
Simplify[sol]


Result is too large fit in the margin of this answer, so I zoomed in

Maple:

ode:=diff(y(x),x\$3)+(1/x)*diff(y(x),x)-y(x)/x^2=0;
dsolve({ode, y(1)= 1, D(y)(1)=0, (D@@2)(y)(1)=1});
simplify(%,wronskian);

y(x) = 2*x-(BesselY(0, 2)-2*BesselY(1, 2))*Pi*sqrt(x)*
BesselJ(1, 2*sqrt(x))+Pi*(-2*BesselJ(1, 2)+BesselJ(0, 2))
*sqrt(x)*BesselY(1, 2*sqrt(x))


Update

Here is Maple result written using Mathematica notation in case someone wants to use to help find simplification to Mathematica solution. Also tried Simplify[mma == maple] suggested in the comments below.

ClearAll[y, x];
ode = y'''[x] + (1/x) y'[x] - (1/x^2) y[x] == 0;
ic = {y[1] == 1, y'[1] == 0, y''[1] == 1};
mma = y[x] /. First@DSolve[{ode, ic}, y[x], x];
maple = 2 x - (BesselY[0, 2] - 2 BesselY[1, 2]) Pi Sqrt[x] BesselJ[1,
2 Sqrt[x]] + Pi (-2 BesselJ[1, 2] + BesselJ[0, 2]) Sqrt[x] BesselY[1,
2 Sqrt[x]];

N[mma /. x -> 3]
(*3.33870606057436*)

N[maple /. x -> 3]
(*3.33870606057435*)

Simplify[mma == maple]
(* too large to post *)


Related question here

-
Interesting question, but what is your experience simplifying this expression with FullSimplify? Perhaps you could do somthenig like this FullSimplify[expr1 == expr2], where expr1 is the result of Mathematica and expr2 - that of Maple? If it didn't work I would start simplifying a few terms of expr1 to get an insight what we can do more. – Artes Nov 18 '13 at 23:57
Try FunctionExpand. It might be able to expand hypergeometric and MeijerG functions. – Vladimir Reshetnikov Nov 18 '13 at 23:57
@Artes my experience simplifying this expression with FullSimplify is that I let it run for few minutes and nothing happened, so I stopped the kernel. I was planning to leave it running overnight under FullSimplify when done using the computer. – Nasser Nov 19 '13 at 0:04
@Artes To use FullSimplify[expr1 == expr2] one needs to know the objective expr2 (which in this case comes from Maple) beforehand. Mathematica should be able to find expr2, or something similar, by itself. If Maple can do it, so can Mathematica! – becko Nov 19 '13 at 1:14

Using FunctionExpand as suggested by Vladimir can improve the simplification, but for some inputs FunctionExpand is very slow in this case MeijerG is the one causing problems. By replacing all occurances of MeijerG with a temporary symbol it is possible to get improved results quicker:

safeApply[f_, expr_, bad_List] := Module[{
temp, cleanRules
},
(* Replace anything bad with temp[i], apply function and unreplace *)
f[expr /. cleanRules] /. Reverse[cleanRules, 2]
]

LeafCount[sol]
(*       5097 *)

LeafCount[Simplify[sol]]//AbsoluteTiming
(* 5.5s  1992 *)

LeafCount[res = Simplify@safeApply[FunctionExpand, sol, {_MeijerG}]]//AbsoluteTiming
(* 8.5s   434 *)

LeafCount[res2 = safeApply[
Simplify[#, TransformationFunctions -> {Automatic, FunctionExpand}] &,
res,
{_MeijerG}]] // AbsoluteTiming
(* 9.5s   380 *)


Still a bit to go to reach LeafCount of 55 in the maple result.

-
very nice. These types of simplifications should be build into Mathematica I think. – Nasser Nov 19 '13 at 2:10