I am trying to solve the following nontrivial equation:

    fun[z_] := (Log[Gamma[n]] + n Log[Gamma[z + 1]] + z n Log[n] - Log[Gamma[n (z + 1)]])/n
    fun[z] /. n -> 8
    m = Table[i, {i, -1, 1, 0.1}] 
    sols = Table[FindRoot[fun'[z] + 0.5 == m[[i]], {z, 1}], {i, Length[m]}]
    sols2 = Table[FindRoot[fun'[z] + 0.5 == m[[i]], {z, sols[[j]]}], {i, Length[m]}, 
        {j,   Length[sols]}]

Basically, I get roots from the equation that are stored in sols. However, I also need all values from sols as starting values for the next iteration that will be executed in sols2. However, this doesn't work and produce the following errors:

    {{FindRoot[\!\(\*SuperscriptBox["fun", "\[Prime]",MultilineFunction->None]\)[z] +0.5 == m[[i]], {z, sols[[j]]}], 
    FindRoot[\!\(\*SuperscriptBox["fun", "\[Prime]",MultilineFunction->None]\)[z] + 0.5 == m[[i]], {z, sols[[j]]}], 
    FindRoot[\!\(\*SuperscriptBox["fun", "\[Prime]", MultilineFunction->None]\)[z] + 0.5 == m[[i]], {z, sols[[j]]}], 

Any suggestion on how to fix this?

  • 1
    $\begingroup$ Are you aware that Mathematica has LogGamma[] built in? $\endgroup$ – J. M.'s discontentment May 16 '12 at 0:25
  • $\begingroup$ @J.M.Thank you for your information. $\endgroup$ – Wayan May 16 '12 at 2:01

Let's polish up your code a bit...

With[{n = 8},
 fun[z_] := (LogGamma[n] + n LogGamma[z + 1] + z n Log[n] - LogGamma[n (z + 1)])/n;
 m = Range[-1, 1, 1/10];
 Map[(z /. First[FindRoot[fun'[z] + 1/2 == #, {z, 1}]]) &, m]

Some notes:

  • You don't really need or want the last step in your question here; a far better route would have been to tweak the options of FindRoot[], as well as finding better initial estimates for your roots.

  • I have elected to use the built-in LogGamma[] function here, which avoids all the problems of a construction like Log[Gamma[z]]. (There's something strange about having to deal with large numbers just to get a reasonably-sized result...)

  • Range[] is as useful as Table[] for generating a list of numbers in arithmetic progression, and should be always kept in mind.

  • Map[] is useful if you want to run some function through the elements of a list.

| improve this answer | |
  • $\begingroup$ Do you show Range and Map rather than Table just to teach it as an option, or do you find is superior in this case? (+1 by the way) $\endgroup$ – Mr.Wizard May 16 '12 at 6:16
  • $\begingroup$ A bit of both, actually. After all, in this specific case there is the much more compact Table[z /. First[FindRoot[fun'[z] + 1/2 == k, {z, 1}]], {k, -1, 1, 1/10}], but the approach here also works if m had been some random array of numbers... $\endgroup$ – J. M.'s discontentment May 16 '12 at 6:28
  • 1
    $\begingroup$ One could also use Table[. . . , {k, m}] in that case. (Telling anyone reading, not you J.M.) $\endgroup$ – Mr.Wizard May 16 '12 at 6:44
  • $\begingroup$ Oh, it's fine. I actually forgot that Table[] can work like a foreach() construct these days. $\endgroup$ – J. M.'s discontentment May 16 '12 at 10:02

Syntactically correct version of sols2 is

 sols2 = Table[FindRoot[fun'[z] + 0.5 == m[[i]], {z, sols[[j, 1, 2]]}], 
 {i, Length[m]}, {j, Length@sols}]

using fun[z] as in Mark's answer.

However, this does not serve your purpose for two reasons:

First, the solution of equation i from sols is used as starting value for the search for the root of equation j (for different i and j) in the second round which is not likely to improve the results. Perhaps, you meant to use a single index as in:

 sols2 = Table[FindRoot[fun'[z] + 0.5 == m[[i]], {z, sols[[i, 1, 2]]}], 
 {i, Length[m]}]

Second, considering the error message you get for both sols and sols2

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

you need to consider playing with the options for FindRoot as J.M. suggested in his answer and comments.

EDIT: Inpecting the plots could suggest an explanation of the numerical results.

Plotting three examples from the 21 equations:

 Row[Plot[fun'[z] + 0.5 - m[[#]], {z, -3, 3}] & /@ {1, 15, 21}]

enter image description here

So, as the region z<=-1 seems hopeless, one needs to restrict root search to the appropriate region (z>-1).

Next, plotting all 21 functions

 Plot[Evaluate[(fun'[z] + 0.5 - m[[#]]) & /@ Range@Length@m], {z, -1, 5}]

enter image description here

we get a hint as to why the first 14 equations are solved without any issue while the last 7 are not. In particular, equation 15 simplifies to

  FindRoot[fun'[z] == 0, {z, 1}] 

Note that Limit[fun'[z] ,z->Infinity]=0 and fun'[z] is strictly increasing on (0,Infinity) (This reference and citations therein on complete monotonicity properties of Gamma and related functions). These observations together with the fact that

 fun'[0] = - 0.513416 < 0

explain why the equation fun'[z]==0 does not have a solution in [0,Infinity]. By the same reasoning, equations 15-21 (which differ from equation 15 by a constant) do not have a solution in [0,Infinity].

| improve this answer | |
  • $\begingroup$ Nice analysis! :) $\endgroup$ – J. M.'s discontentment May 16 '12 at 10:04
  • $\begingroup$ @kguler Thank you. Your deep analysis helps me in understanding more about my problem. $\endgroup$ – Wayan May 16 '12 at 11:44

Your fun depends upon n and, though you plugged in n->8 at line 2, your subsequent code doesn't refer to that line. Thus, FindRoot receives input that doesn't return numerical output. That's the source of the error.

I don't see the purpose behind the two steps, as in both steps you're solving equations of the form fun'[z]+0.5==c, where c is constant. These are easily solved.

fun[n_, z_] := (Log[Gamma[n]] + n Log[Gamma[z + 1]] + 
  z n Log[n] - Log[Gamma[n (z + 1)]])/n;
f[z_] = fun[8, z];
FindRoot[f'[z] + 1/2 == -3/10, {z, 0}]

{z -> -0.321733}

| improve this answer | |
  • $\begingroup$ It would seem that OP is trying to polish the results in sols... $\endgroup$ – J. M.'s discontentment May 16 '12 at 0:44
  • $\begingroup$ @Mark. I need the second step,i.e.sols2 to refine the root to certain range. $\endgroup$ – Wayan May 16 '12 at 2:16
  • $\begingroup$ @Wayan: the refinement doesn't work very well here in this case. You'd be better off changing FindRoot[]'s options, or even trying to find a better starting estimate for FindRoot[]. $\endgroup$ – J. M.'s discontentment May 16 '12 at 2:34
  • $\begingroup$ @J.M. Thank you for your suggestion. $\endgroup$ – Wayan May 16 '12 at 2:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.