# Solving functional equation

Can I use Mathematica to solve a functional equation?

For instance how to solve the following functional equation by using Mathematica?

$$H(x+1)=xH(x)+\frac{1}{\Gamma(1-x)}$$

Solution is the Hadamard's gamma function

$$H(x) = \frac{\Psi(1-x/2)-\Psi(1/2-x/2)}{2\Gamma(1-x)}$$

• There is a typo: you probably meant $H(x+1)=xH(x)...$ – BlacKow Jan 10 '17 at 16:56
• surely there is more to this, that equation alone doesn't uniquely define the result. – george2079 Jan 10 '17 at 18:06
• Qualifying a bit, it is unique on the integers, but you can define H(x) on 0<x<1 as anything you like so long as H(1)==1 and H(0) is finite. – george2079 Jan 10 '17 at 19:14

RSolve can solve some functional equations. For example, if you replace $1/\Gamma[1-x]$ with $x^2$, Mathematica can solve the equation in under half a second on my machine:

RSolve[hh[x + 1] == x hh[x] + x^2, hh[x], x]

(* {hh[x] -> C[1] Pochhammer[1, -1 + x] + ((E x Gamma[-1 + x, 1])/Gamma[-1 + x]
- (1 + E (-2 + x) Gamma[x, 1])/Gamma[x]) Pochhammer[1, -1 + x]} *)


Note the presence of the undetermined constant C[1]; in general, the solutions to arbitrary functional equations will have some number of undetermined constants (as referred to by @george2079 in the comments.) Boundary/initial conditions can also be included that allow Mathematica to solve for these constraints, for example:

RSolve[{hh[x + 1] == x hh[x] + x^2, hh[1] == 1}, hh[x], x]

(*{hh[x] -> -1 + Gamma[x] - E x Gamma[-1 + x, 1] + E x^2 Gamma[-1 + x, 1] + 2 E Gamma[x, 1] - E x Gamma[x, 1]} *)


Unfortunately, Mathematica has a much harder time with your original equation. The code to ask Mathematica to solve this equation would be

RSolve[hh[x + 1] == x hh[x] + 1/Gamma[1 - x], hh[x], x]


but I issued this command when I started writing this answer and it hasn't stopped thinking since. I'll update this answer if it comes up with a response. However, as noted above, if Mathematica obtains an answer it will most likely be a general answer that has the Hadamard gamma function as a special case. (EDIT: Mathematica ran for about four hours without producing any output.)

I don't know how to solve the functional equation; however, you can use Mathematica to verify the equation.

Working from the definition of H[x]

Clear[H];

(H[x_] = 1/Gamma[1 - x] D[Log[Gamma[1/2 - x/2]/Gamma[1 - x/2]], x] //


For positive integer values of x

seq = Table[Limit[H[x], x -> n], {n, 10}]

(*  {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880}  *)

FindSequenceFunction[seq][x] // FunctionExpand

(*  Gamma[x]  *)

Show[
Plot[H[x], {x, -2, 5}, PlotRange -> All],
DiscretePlot[Gamma[x], {x, 1, 5}]]


Verifying the functional equation

H[x + 1] == x H[x] + 1/Gamma[1 - x] // FullSimplify

(*  True  *)