Problem of `Integrate` and `DSolve` with inhomogeneous Legendre differential equation

Question

I have the following inhomogeneous Legendre differential equation $(m=n=2)$:

 deq=(6 - 4/(1 - z^2)) h2[z] - 
  2 z Derivative[1][h2][z] + (1 - z^2) Derivative[2][h2][z] == (
  6 (-1 + z + z^2) μ^2)/(M^4 (-1 + z^2)^2) + (
  3 μ^2 Log[(-1 + z)/(1 + z)])/(M^4 (-1 + z))

where the variable $z>1$. I can not get Mathematica to give me the simplest/correct solutions. I think it has to do with those logarithms. I know that the homogeneous solution is

solH=A*LegendreQ[2,2,3,z]+B*LegendreP[2,2,3,z]

I need mathematica's type 3 Legendre functions because for $z>1$ the type 3 functions are on the right branch. With DSolve however I can only get the expanded forms of type 1:

DSolve[deq[[1]] == 0, h2[z], z]

gives back

 {h2[z] -> -3 (-1 + z^2) C[1] + C[2] (((1 - z^2) (5 z - 3 z^3))/(-1 + z^2)^2 + 
      3 (1 - z^2) (-(1/2) Log[1 - z] + 1/2 Log[1 + z]))}

which is not real for $z>1$. Is there a way to tell DSolve to work with $z>1$ or to search for real solutions under that assumption? $Assumptions={z>1} alone does not do the trick.

Apart from that I have a bigger problem with the particular solution: DSolve[deq, h2[z], z][[1]] gives a ridiculous partial solution which is incredibly long and complex. I want a real particular solution; I want a solution that is real for all $z>1$. I know from literature that there is one (plugging that one in proves that it is indeed one):

{h2p[z] -> -3 μ^2/(16 M^4) (6 z^2 + 3 z - 6 - (4 z^2 + 2 z)/(z^2 - 1)) - 3 μ^2/(32 M^4) (3 z^2 - 8 z - 3 - 8/(z^2 - 1)) Log[(z - 1)/(z + 1)] + 3 μ^2/(16 M^4) (z^2 - 1) Log[(z - 1)/(z + 1)]^2}

Again is there a way to tell DSolve to produce a real solution for $z>1$? I also tried constructing that particular solution using the Variation of Parameters method

WLQ223LP223 = 
 Wronskian[{LegendreQ[2, 2, 3, z], LegendreP[2, 2, 3, z]}, z]
-LegendreQ[2, 2, 3, z]*
    Integrate[deq[[2]]*LegendreP[2, 2, 3, z]/WLQ223LP223, z] + 
   LegendreP[2, 2, 3, z]*
    Integrate[deq[[2]]*LegendreQ[2, 2, 3, z]/WLQ223LP223, z] // 
  FunctionExpand // Simplify

Mathematica finds an analytical solutions but they do not reduce to the particular solution I want: the involved integrals become very lengthy and again complex.

Anyone around here who has some experience with such differential equations in Mathematica?

I send the problem to a friend who has Maple and Maple's dsolve finds the "correct"/simple partial solution, when one prepends assume(z > 1). It does so using the Variation of Parameters method and the solutions of the homogeneous equation.

Answer 1

Personally, I prefer @corey979's approach (+1), but with the following tweaks.

solExpr = 
  Assuming[z > 1 && M ∈ Reals && μ ∈ Reals,
   DSolveValue[deq, h2[z], z] /. {Log[-1 + z^2] -> Log[-1 + z] + Log[1 + z], 
        Log[(-1 + z)/(1 + z)] -> Log[-1 + z] - Log[1 + z]} // Re // 
     ComplexExpand[#, TargetFunctions -> {Re, Im}] & // Simplify
   ];

sol = h2 -> Function @@ {z, solExpr}

deq /. sol // Simplify[#, z > 1 && M ∈ Reals && μ ∈ Reals] &
(*  True  *)

Plot[h2[z] /. sol /. {M -> 1, μ -> 1, C[1] -> -1, C[2] -> 1},
 {z, 1, 5}]

The main tweaks are including all assumptions to make all parameters real in both DSolve and Simplify and the use of TargetFunctions. It's also easier to simplify an algebraic expression; hence DSolveValue[]. The rules

{Log[-1 + z^2] -> Log[-1 + z] + Log[1 + z], 
 Log[(-1 + z)/(1 + z)] -> Log[-1 + z] - Log[1 + z]}

help with simplification; they are not important and may be omitted.

The main reason I prefer this is that the real and imaginary parts are independently solutions of a real ODE, a key point made by corey979. It just seems a simpler approach, when it simplifies down nicely. In this case it does, and we get a general solution that spans the solution space.

Answer 2

About the homogeneous part:

Let us first define the (comparison) solution with the Legendre functions:

solLeg[z_] := B*FunctionExpand@LegendreQ[2, 2, 3, z] + A*LegendreP[2, 2, 3, z]

LegendreQ doesn't want to display anything friendly in this case, so I used FunctionExpand to get a comprehensible result:

FullSimplify@solLeg[z]

Next, let's solve the equation:

solEq = DSolve[deq[[1]] == 0, h2[z], z, Assumptions -> z > 1] /. {C[1] -> A, C[2] -> B}

I used Assumptions -> z > 1 to specify the domain and changed the default constant of integration to A and B just for clarity of the display. Note that there's an $-i\pi$ term before the Logs, but the two functions are quite similar to each other.

Now, the crucial points are: (1) in general, as can be seen above, the solution is complex, even though the domain z>1 is real; this is not a suspicious behaviour by itself, differential equations and their solutions are just like that; (2) the real and imaginary parts of the solution are together and separately the solutions of the equation; one just has to extract the real part to get the real solution.

I define a helping expression:

help = h2[z] /. solEq[[1, 1]]

to make some manipulations on:

sol = Simplify @ PowerExpand @ ComplexExpand @ Re @ help

ComplexExpand @ Re extracts the real (without the imaginary unit $i$) part of help; PowerExpand is used to convert $\log(x^2)$ into $2\log x$, and Simplify as usual. We get

Finally,

sol + FullSimplify@solLeg[z] // FullSimplify

0

meaning that what Mathematica returns - solEq - as the solution of the differential equation is indeed consistent with the solution in a form of Legendre functions solLeg.

---

About the inhomogeneous case:

If one runs

FullSimplify@
 PowerExpand@
  ComplexExpand@
   Re@First[
     h2[z] /. 
       FullSimplify@
        DSolve[deq, h2[z], z, Assumptions -> z > 1] /. {C[1] -> A, 
       C[2] -> B}]

it returns a four-line long expression with Logs and Arg, and of course with the constants of integration A and B (it's uninformative enough to not post it here; it should be noted that it takes about 7 seconds to obtain this expression). The particular real solution from the literature, mentioned by the OP:

has many terms similar to those from the above computation, but no constants of integration in it so I suspect that some initial conditions were employed to obtain it. Once they are provided I expect that a similar treatment as in the homogeneous case could also show that the solutions are equivalent; or maybe it will turn out that DSolve will immediately give the expected solution.

Answer 3

Based on the limited response to question 126132, it appears that there is no general method for causing DSolve to produce results that are free of explicit complex terms. However, it often is possible to transform the DSolve constants of integration to eliminate complex terms, as in 126072. For the solution to the homogeneous ODE in the present question,

Flatten@DSolve[deq[[1]] == 0, h2[z], z];
Collect[h2[z] /. % /. C[1] -> C[1 ] + C[2] Pi I/2, C[_], FullSimplify[#, z > 1] &]
(* (3 - 3 z^2) C[1] + ((z (-5 + 3 z^2) - 3 (-1 + z^2)^2 ArcCoth[z]) C[2])/(-1 + z^2) *)

is free of complex terms. Eliminating complex terms from the corresponding solution of the inhomogeneous ODE is more involved but can be accomplished by the same general approach. Starting as before,

s1 = Collect[h2[z] /. Flatten@DSolve[deq, h2[z], z] /. 
    C[1] -> C[1 ] + C[2] Pi I/2, C[_], FullSimplify[#, z > 1] &]
(* (1/(16 M^4 (-1 + z^2))) 3 μ^2 (-6 + z (5 + z (16 - 3 z (1 + 2 z))) + 
   I π (1 + z (-5 + z (-2 + z (3 + z)))) + (2 - 3 I π (-1 + z^2)^2 + 
   2 z (4 + z (-10 + z (-4 + 5 z)))) ArcCoth[z] + 
   4 (-1 + z^2)^2 ArcCoth[z]^2 - 7 (-1 + z^2)^2 ArcTanh[z]) + 
   (3 - 3 z^2) C[1] + ((z (-5 + 3 z^2) - 3 (-1 + z^2)^2 ArcCoth[z]) C[2])/(-1 + z^2) *)

Because ArcTanh[z] is complex for z > 1, it first should be eliminated, which can be accomplished using the identity

FullSimplify[ArcCoth[z] - Pi I/2 == ArcTanh[z], z > 1]
(* True *)

s2 = FullSimplify[s1 /. ArcTanh[z] -> ArcCoth[z] - Pi I/2]
(* (1/(32 (-1 + z^2)))((1/(M^4)) 3 μ^2 (-2 (6 + z (-5 + z (-16 + 3 z + 6 z^2))) + 
   I π (9 + z (-10 + 3 z (-6 + z (2 + 3 z)))) + 2 ArcCoth[z] (-5 - 3 I π 
   (-1 + z^2)^2 + z (8 + z (-6 + z (-8 + 3 z))) + 4 (-1 + z^2)^2 ArcCoth[z])) - 
   96 (-1 + z^2)^2 C[1] + 32 (z (-5 + 3 z^2) - 3 (-1 + z^2)^2 ArcCoth[z]) C[2]) *)

Finally, the explicit complex numbers are eliminated by finding the values of {C[1], C[2]} that cancel those complex numbers.

t1 = Coefficient[s2, Pi] Pi;
t2 = Coefficient[s2, {C[1], C[2]}].{C[1], C[2]};
{C[1], C[2]} /. Flatten@Solve[
    Thread[CoefficientList[t1 - t2, ArcCoth[z]] == 0], {C[1], C[2]}];
Thread[Rule[{C[1], C[2]}, {C[1], C[2]} - %]]
(* {C[1] -> (9 I π μ^2)/(32 M^4) + C[1], C[2] -> -((3 I π μ^2)/(16 M^4)) + C[2]} *)

s3 = Collect[s2 /. %, C[_], Simplify]
(* (1/(16 M^4 (-1 + z^2))) 3 μ^2 (-6 + 5 z + 16 z^2 - 3 z^3 - 
   6 z^4 + (-5 + 8 z - 6 z^2 - 8 z^3 + 3 z^4) ArcCoth[z] + 4 (-1 + z^2)^2 ArcCoth[z]^2) + 
   (3 - 3 z^2) C[1] + ((z (-5 + 3 z^2) - 3 (-1 + z^2)^2 ArcCoth[z]) C[2])/(-1 + z^2) *)

free of complex terms, as desired.