3

For example, C=0,

solve('psi(x)=0')

ans =
-226.83295306016122662496413158295

psi(ans)

???Error using ==> psi
Input must be single or double.

I cannot get the right answer

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
user3114396
  • 31
  • 2
  • What version of Matlab are you using? I have no trouble in R2013a. The symbolic version of `psi` is only supported for R2011b+. – horchler Dec 18 '13 at 15:48
  • Thank you very much, when I tried it in R2013a, no trouble happened. However, I don't know how to solve the equation digamma(x)=C, in which C is a variable, not a const. Can you help me with this? – user3114396 Dec 19 '13 at 07:51
  • You essentially want an analytical formula for the inverse of the digamma function? I'm not aware one exists. Neither Matlab nor Mathematica can solve for it even with assumptions. The appraise place to ask about this would be [Math.StackExchange](http://math.stackexchange.com). But I think that you're stuck using numeric or symbolic numeric methods unless you go to an approximate analytic solution. – horchler Dec 19 '13 at 15:11
  • No...I just want to know how to solve the equation, for example, you can solve psi(x)=10 in matlab through x=solve('psi(x)=0'), but how to solve x=solve('psi(x)=C'), in which C is a variable, not a const. – user3114396 Dec 20 '13 at 02:25
  • That's exactly the definition of the analytical [inverse](http://en.wikipedia.org/wiki/Inverse_function) of the digamma function. I can't add anything to my previous comment. – horchler Dec 20 '13 at 14:41

1 Answers1

1

Interesting...this seems like a bug in solve to me...Whatever value I try to solve for, I always get a weird value of around -227. Even when I try to trick MATLAB by giving an approximation of the digamma, I get the same result or worse:

>> solve('(gamma(x+0.01)-gamma(x))/0.01/gamma(x)=0')
ans = 
    matrix([[-226.83790783643886637282996154237]])

>> solve('(gammaln(x+0.01)-gammaln(x-0.01))/0.02 = 0')
??? Error using ==> mupadmex
Error in MuPAD command: cannot differentiate equation [numeric::fsolve]

The following numerical approach works: 

%// value of the digamma to solve for
Y = -10; 

%// Solve using numerical scheme 
X = fsolve(@(x)psi(max(0,x)) - Y, exp(Y))

%// Check solution: psi(X) ≈ Y
psi(X)
Rody Oldenhuis
  • 37,726
  • 7
  • 50
  • 96
  • +1 Same here (2010b). I always get a value around `-226`. Interestingly, `solve('exp(psi(x))=1')` gives around `-14` – Luis Mendo Dec 18 '13 at 13:02
  • @LuisMendo: hmmmmmmmmmmmmmmm....I really don't see what could be going on...The offset is not an "obvious" value, and doesn't seem to follow some simple transformation rule...The fact that there is no warning or error just makes it smell like "bug". – Rody Oldenhuis Dec 18 '13 at 13:28
  • 2
    @LuisMendo, et al.: Seems to work perfectly fine in R2013a. According to the [release notes](http://www.mathworks.com/help/symbolic/release-notes.html#bs17m4_-1), `sym/psi` was only added in R2011b, so you may be using a pre-release version with plenty of issues. – horchler Dec 18 '13 at 15:41
  • And `solve` doesn't like `gammaln` simply because it's a numeric function. There is no `sym/gammaln`. The MathWorks in their laziness didn't bother to make function names match up so you need to use [`lngamma`](http://www.mathworks.com/help/symbolic/mupad_ref/lngamma.html) instead. The same annoyance holds for the incomplete gamma function, `igamma`, but not only is the name different, the arguments are switched, and it's the unregularized version. – horchler Dec 18 '13 at 15:46