MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Solving an equation

  • To: mathgroup at
  • Subject: [mg51442] Re: Solving an equation
  • From: "David W. Cantrell" <DWCantrell at>
  • Date: Sun, 17 Oct 2004 03:06:10 -0400 (EDT)
  • References: <cklle4$eqe$> <ckqn9i$ngq$>
  • Sender: owner-wri-mathgroup at

"David W. Cantrell" <DWCantrell at> wrote:
> jujio77 at (Scott) wrote:
> > I have an equation, Gamma[a+I b] = some complex number.  I need to
> > solve this equation for a and b.
> Typically, there will be infinitely many solutions. See below.
> > I can write a+I b as z, but it can't
> > be solved through NSolve.  What I have been doing is a double do loop
> > for a and b and getting some number.  Then I compare this to the
> > number I have.  Then I narrow down my possibilities for a and b, and
> > go through the process again.  Does anyone know of a better way to do
> > this problem?  Is there a way to have Mathematica compare each result
> > of the do loop to a given value, and given certain conditions spit out
> > an answer for a and b?
> Certainly look at the response already given by Daniel Lichtblau, but let
> me mention a very different approach which might be of use. Some time
> ago, I developed an approximation for the inverse (principal branch) of
> the real gamma function; see "Inverse Gamma Function" at
> <>. [At the time I
> wrote that, I was unaware that Daniel had previously suggested something
> similar.] Bizarrely, until I saw your question, I had never considered
> using my approximation in the complex domain, but it works nicely!
> The main "adaptation" I've made to my previous approximation is to allow
> different braches of the logarithm to be used, since we're now working
> with complex numbers and multiple solutions. AIG stands for
> ApproximateInverseGamma.
> L[z_, n_] := Log[z/Sqrt[2*Pi]] + 2*n*Pi*I;
> AIG[z_, n_] := L[z, n]/ProductLog[L[z, n]/E] + 1/2
> For our first examples, let's use the same value, 3 + 5*I, used by Daniel
> earlier in this thread.
> In[5]:= N[AIG[3 + 5*I, 0]]
> Out[5]= 4.040766662437202 + 0.809458889448528*I
> In[6]:= Gamma[%]
> Out[6]= 2.9542538687291082 + 4.952066867619615*I
> To get this solution far more precisely
> In[7]:= FindRoot[Gamma[z] == 3 + 5*I, {z, 4}]
> Out[7]= {z -> 4.048806909342241 + 0.8061219366736553*I}
> but, _considering its simplicity_, AIG works reasonably well here.
> Note that 4.04... + 0.809...*I is not an approximation of one of the two
> solutions mentioned by Daniel. To get an approximation of one of those,
> we use a different branch:
> In[8]:= N[AIG[3 + 5*I, 1]]
> Out[8]= 5.244221150676627 + 4.373221814339161*I
> The other solution which Daniel gave apparently cannot be approximated by
> AIG. Nonetheless, AIG provides a very easy way to approximate as many
> solutions as desired by using different branches.

For a more complete answer, I'll now mention how still other solutions
(such as the other one given by Daniel) can be very easily approximated.
This will give another infinite family of approximate solutions.

Toward the end of the article at the above link, I had said
"Nice _algebraic_ asymptotic approximations can be obtained fairly
easily for all nonprincipal branches. (For all nonpositive n, the gamma
function has a simple pole at n.) For example, for the -1 branch, the
function 1/(g+x), where g is the Euler gamma constant (approximately
0.577216), works well when |x| is large.
For all other branches, the simplest nice asymptotic approximations are
rational functions having both numerator and denominator of degree 1. As
examples: For the -2 branch, we have the approximation (2-g+x)/(g-1-x),
and for the -3 branch, 4(g-2+2x)/(3-2g-4x)."

So now let's see how to get an approximation for that other solution
previously mentioned by Daniel.

Of course, we could use something like
In[20]:= FindRoot[Gamma[z] == 3 + 5*I, {z, 0.1}]
Out[20]= {z -> 0.09103227139407653 - 0.13362798410019328*I}

But an alternative, if we can settle for a simple approximation, is to use

AIGNP[z_, n_] := -n + 1/(EulerGamma + (-1)^n*(n!*z + StirlingS1[n+1,2]/n!))

where AIGNP stands for ApproximateInverseGammaNearPole at -n for
nonnegative integer n.

We then get
In[25]:= N[AIGNP[3 + 5*I, 0]]
Out[25]= 0.09464416872334132 - 0.13228747941025595*I

[By the way, substantially better, but messier, approximations than AIGNP
are readily available near the poles. For example, instead of using AIGNP
above, we could have used
In[26]:= N[2/(EulerGamma + z + Sqrt[-EulerGamma^2 - Pi^2/3 + 2*EulerGamma*z
+ z^2])]/. z -> 3 + 5*I
Out[26]= 0.0905892579464156 - 0.1332971242252317*I
In[27]:= Gamma[%]
Out[27]= 3.0051527268465104 + 5.020664268483219*I
However, for simplicity, we'll stay with just AIGNP for the rest of this

Note that AIGNP agrees with the three simple approximations mentioned in
the quotation above:

In[28]:= Table[Together[AIGNP[z, n]], {n, 0, 2}]
{1/(EulerGamma + z), (2 - EulerGamma + z)/(-1 + EulerGamma - z),
-((4*(-2 + EulerGamma + 2*z))/(-3 + 2*EulerGamma + 4*z))}

AIGNP allows us to get as many approximate solutions as desired. For

In[29]:= Table[N[AIGNP[3 + 5*I, n]], {n, 0, 7}]
{0.09464416872334132 - 0.13228747941025595*I,
 -1.0932246259350227 + 0.136182442140838*I,
 -1.9596335538408043 - 0.07950508472241272*I,
 -3.015152775376722 + 0.02360721247808417*I,
 -3.9963605517971357 - 0.006195343055295479*I,
 -5.000736925510258 + 0.0012224158911233105*I,
 -5.999877501029459 - 0.00020434212160696532*I,
 -7.000017508100925 + 0.000029176278728846562*I}
In[30]:= Gamma[%]
{3.0978502145116456 + 4.8903488180579275*I,
 3.130381565322433 + 4.81893763567967*I,
 3.0336452404650402 + 4.921088107165388*I,
 3.005898092059356 + 4.99124348680135*I,
 3.0004015914098887 + 4.999305766120647*I,
 3.000018510086481 + 4.999969393030145*I,
 3.0000005650151262 + 4.999999056961428*I,
 3.00000001254731 + 4.999999979090844*I}

Apparently, as n increases, so does the accuracy of the approximation.

So AIG and AIGNP give approximations of solutions in infinite disjoint
families. But do we now have a complete set of approximations, or are
there other solutions lying outside of those families? I suspect the

David Cantrell

  • Prev by Date: Re: Re: Sorting a list of pairs on the second elements
  • Next by Date: Re: Re: Re: Sorting a list of pairs on the second elements
  • Previous by thread: Re: Solving an equation
  • Next by thread: Re: Solving an equation