6
eq = x^2 f''[x] + x f'[x] - 4 f[x] == 0;
bc = f[∞] == 0; (* boundary condition *)

sol = DSolve[ {eq, bc}, f[x], x]

(* {{f[x] -> 0}} *)

The correct solution is f[x]=c/x^2. Why does DSolve give a wrong answer?

gwr
  • 13,452
  • 2
  • 47
  • 78
KDH
  • 91
  • 1
  • Where is the actual equation? I see no ODE here. – zhk Mar 08 '18 at 09:33
  • I have an image there. The equation is x^2 f '' (x) + x f'(x) - 4 f(x) = 0 with boundary condition f(Infinity)=0. DSolve gives f(x)=0 but the correct solution is f(x) = c / x^2 – KDH Mar 08 '18 at 09:36
  • Clearly, the answer is correct, trivially so. – TimRias Mar 08 '18 at 09:39
  • f(x)=0 is a solution of the ODE but f(x)=c/x^2 is the general solution. DSolve missed the general solution. Why? – KDH Mar 08 '18 at 09:42

5 Answers5

8

The boundary condition f[Infinity]==0is the problem:

You might solve your problem using the general solution 

sol = DSolve[{x^2 f''[x] + x f'[x] - 4 f[x] == 0}, f, x][[1]]
TrigToExp[ f[x] /. sol]
(* C[1]/(2 x^2) + 1/2 x^2 C[1] - (I C[2])/(2 x^2) + 1/2 I x^2 C[2] *)
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55
6

Alternatively to gwr's extension, one could examine the general solution from Ulrich Neumann symbolically:

sf = C[1]/(2 x^2) + 1/2 x^2 C[1] - (I C[2])/(2 x^2) + 1/2 I x^2 C[2]
Limit[sf, x->Infinity]

(C[1] + I C[2]) Infinity

We need this expression to equal 0, so set C[2]->I C[1]:

sf /. C[2] -> I C[1]

C[1]/x^2

To get to the heart of why this happens, I imagine it's because DSolve evaluates Infinity preemptively. Even after taking the limiting case specifically and seeing the term containing Infinity, a bit of a logical leap is required to say that the only way this will work is if (C[1] + I C[2]) is 0. Even after that leap that term is Indefinite, not 0 exactly.

eyorble
  • 9,383
  • 1
  • 23
  • 37
4

Building upon Ulrich Neumann's answer you may approach the solution by including the boundary condition for arbitrarily large x:

eq = x^2 f''[x] + x f'[x] - 4 f[x] == 0;

sol = DSolve[ {eq, f[1.*^60] == 0}, f[x], x] // RightComposition[
    First,
    TrigToExp, 
    FullSimplify,
    N
]

$\text{f}[x] \rightarrow 0. +\frac{1. \text{C[1]}}{x^2}$

gwr
  • 13,452
  • 2
  • 47
  • 78
  • +1 for introducing me to RightComposition! – Thies Heidecke Mar 08 '18 at 13:01
  • @ThiesHeidecke You are welcome - I must admit that I tend to use it excessively -- feeling vindicated, though, by looking at more structured code ... ;-) – gwr Mar 08 '18 at 13:07
  • 1
    Through my previous exposition to functional programming with Haskell i was already a fan of point free style and liked to use Composition a lot, it's good to know there are more options that feel more natural in some situations! – Thies Heidecke Mar 08 '18 at 13:14
  • How is RightComposition here better than sol=... // First // TrigToExp // FullSimplify // N? – Ruslan Mar 08 '18 at 18:26
  • 1
    @Ruslan Essentially, there is no difference but the block-writing with , makes for better readability and Composition[ g, h ] can be immediately treated as a function while ... // g // h can not. – gwr Mar 08 '18 at 19:00
3

Further to the other answers: Mathematica can find the general solution so long as you ask it to solve for the boundary condition $f(a) = 0$ at a finite (but arbitrary) $a$, and only then take the limit $a \to \infty$:

sol = DSolve[{x^2 f''[x] + x f'[x] - 4 f[x] == 0, f[a] == 0}, f[x], x][[1]]
Limit[f[x] /. sol, a -> \[Infinity]]

(* {f[x] -> C[1] Cosh[2 Log[x]] - C[1] Coth[2 Log[a]] Sinh[2 Log[x]]} *) 
(* C[1]/x^2 *)

See also my answer to DSolve not satisfying initial conditions, which dealt with a similar problem of DSolve missing certain solutions.

Michael Seifert
  • 15,208
  • 31
  • 68
1

Carrying the above solutions further:

eq = x^2 f''[x] + x f'[x] - 4 f[x] == 0;

DSolve[eq, f[x], x] // Flatten
(* {f[x] -> C[1]*Cosh[2*Log[x]] + I*C[2]*Sinh[2*Log[x]]} *)

f[x_] = f[x] /. % /. {C[1] -> c1, C[2] -> c2/I}
(* c1*Cosh[2*Log[x]] + c2*Sinh[2*Log[x]] *)

f[x] // TrigToExp // Collect[#, x^_] &
(* x^2*(c1/2 + c2/2) + (c1/2 - c2/2)/x^2 *)

f[x_] = % /. {c1/2 + c2/2 -> c1, c1/2 - c2/2 -> c2}
(* c1*x^2 + c2/x^2 *)

Applying the bc, for f to be 0 at x == infinity, c1 must be 0.

Bill Watts
  • 8,217
  • 1
  • 11
  • 28