2
$\begingroup$

I have an expression like the following (anyway to input latex?),

$$ z^{m}\left( \ln\left(\frac{z}{z-1}\right)-\sum_{k=1}^m \frac{z^{-k}}{k} \right)$$

for m integer, z complex

It seems if I code it the following way, it gives noticeable round-off error in case $|z| \approx 1$.

    Z_1=1/Z

    ZV1=LOG(Z/(Z-1))
    ZK=1.0_q
    DO K=1,M;          ZK=ZK*Z_1
      ZV1=ZV1-ZK/K
    ENDDO
    ZV1=ZV1/ZK

Any suggestions on a better way to code this expression? Thanks,

Added note:

I have compared with gfortran built-in function Log and Mathematica NIntegrate the following two equivalent expressions,

$$ \ln\left(\frac{z}{z-1}\right)$$

and

$$ \int_{0}^{1}\frac{1}{z-x} dx $$

for $z$ to be complex, and found that these two can differ as big as $10^{-15}$ at small $|z|$. (I pushed NIntegrate to a PrecisionGoal of $10^{-17}$). This difference seems to be responsible for the error of the expression I was seeking help (I can upload a figure to support my claim made here if possible).

Do you think this is likely the case? Which result should I trust, from the built-in function or from the numerically converged results? If use Taylor series to help improve numerical convergence, how to treat the case with $|z|\sim 1$, which does not validate any series expansion?

Notes added:

Besides the possibility that Mathematica didn't do a good enough job with Log function, it could equally be that NIntegrate does not reach the prescribed precision goal for some input z values. I am trying to clarify this issue with Mathematica and see how they respond.

$\endgroup$
  • $\begingroup$ You can write inline latex enclosed in pairs of $, and display latex in double $, so $z$ and $$z^m(\cdots)$$. $\endgroup$ – Steve Jan 29 '16 at 11:18
  • 1
    $\begingroup$ It might be best to calculate your sum and your $z^m\ln(z/(z-1))$ terms separately then subtract them after. For the sum try to add the smallest terms first; it might mean instead summing the $k=M$ term first. (I may have misunderstood the issue) $\endgroup$ – Steve Jan 29 '16 at 11:28
  • $\begingroup$ This site supports MathJax for posting mathematical expressions in $\LaTeX$. That Help Page links to a homegrown tutorial at Meta.Math.SE. $\endgroup$ – hardmath Jan 29 '16 at 13:48
  • $\begingroup$ I've tried converting you expression to $\LaTeX$. Please check it over to see if I introduced some mistakes. $\endgroup$ – hardmath Jan 29 '16 at 13:59
  • 1
    $\begingroup$ I guess it could be related to the built-in log function. But what confuses me is that usually we could trust the quality of these internal functions. Maybe not the case? More information is provided in the revised post. $\endgroup$ – bsmile Jan 30 '16 at 6:59
3
$\begingroup$

If I properly converted the expression to be calculated into math notation, then what we are dealing with is subtracting off the $m$ leading terms (ignoring a constant of zero) in the power series expansion of $\ln(1 - z^{-1})$ (corresponding to an expansion of $\ln(1+z)$). This would cause subtractive cancellation and thus amplification of rounding errors relative to the magnitude of the difference.

$$ \ln(1 - z^{-1}) = \sum_{k=1}^\infty \frac{z^{-k}}{k} $$

The difference between this power series and its truncation at the $k=m$ index is then scaled by $z^m$. But multiplying by $z^m$ for $|z| \approx 1$ should not have substantial impact on the relative error of the final result.

My suggestion is to implement a summation of the tail of the power series, the portion left after the leading terms are subtracted off:

$$ \sum_{k=m+1}^\infty \frac{z^{m-k}}{k} $$

I will turn my hand to implementing this for complex $z$ in a neighborhood of $1$, but note that the Question asks for something different: a computation valid for $z$ in a neighborhood of the unit circle!

Here lies a serious difficulty, in that the natural logarithm is a multi-valued function in the complex plane, and any valid branch cut will intersect the unit circle.

The principal branch of the natural logarithm is formed to agree with the real natural logarithm on the positive real axis. The branch cut is then chosen (as a matter of convention) to be the origin and the negative real axis (in the complex plane).

Therefore some clarification is needed of where in the neighborhood of the unit circle (annulus) the computation is expected to be valid (which branch of the natural logarithm are we going to compare to for accuracy).

$\endgroup$
  • $\begingroup$ Thanks for your comment. The first series expansion works for |z| >>1, which is not the case I am asking here. For the second expression, I agree that multiplying a constant should not affect the relative error, which is why I didn't code it that way. Actually, the expression shows more significant error as |z| gets smaller, as I can pain-takingly implement a Taylor expansion to help with that, I thus focus my question around the region of |z|~1. $\endgroup$ – bsmile Jan 30 '16 at 7:14
  • $\begingroup$ I am not very familiar with branch cut. But different branch cuts should bring up a multiple of $\pi$, instead of a tiny number, thus I don't know whether a choice of branch cut is relevant here. To avoid confusion, let I stick with the branch cut convention of Mathematica. $\endgroup$ – bsmile Jan 30 '16 at 7:15
  • $\begingroup$ If z is very close to 1 or 0, then I would not mind it much as the expression will be divergent anyway. My concern is as you said when z falls nearby the convergence circle yet away from 0 or 1. The expression would still be well-defined, yet we don't have the convenience of series expansion. $\endgroup$ – bsmile Jan 30 '16 at 7:16
  • $\begingroup$ Thanks for the clarification. I'm pretty sure Mathematica will use the standard branch cut, so I'll post a picture to illustrate the region of subtractive cancellation. $\endgroup$ – hardmath Jan 30 '16 at 10:56
1
$\begingroup$

You could try asking Herbie, a piece of software that was designed to find more accurate formulas http://herbie.uwplse.org/

$\endgroup$

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.