Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How to code the expression with least roundoff error

  1. Jan 29, 2016 #1
    The expression to code is: (z complex, m some positive integer)

    z^m ( ln(z/(z-1)) - sum( 1/(k z^k), {k=1,m} ) )

    The way I code is (in fortran) (in case z<2)


    DO K=1,M; ZK=ZK*Z_1

    But it gives quite big error by this bruteforce way of coding. Can anybody help me out to suggest a better way to code it?

  2. jcsd
  3. Jan 29, 2016 #2


    User Avatar

    Staff: Mentor

    I would do the sum first, then subtract it from the logarithm. Another possibility is to calculate
    z^m \ln \frac{z}{z-1} - \sum_{k=1}^m \frac{z^{m-k}}{k}
    Depending on the value of z and m, one approach might be better than the other.
  4. Jan 29, 2016 #3


    User Avatar
    Homework Helper

    Summing smaller numbers first will reduce truncation error, so this may be better. I assume the common factor z^m can be placed outside the sum, even if z is complex?

    ##z^m ( \ln \frac{z}{z-1} - \sum_{k=m}^1 \frac{z^{-k}}{k})##
  5. Jan 29, 2016 #4
    Thanks for your help. I wonder whether the issue comes from Log[z/(z-1)]. To test it, I evaluated two equivalent expressions with Mathematica making use of the following identity

    ln(z/(z-1))=∫^{1}_{0} 1/(z-x) dx

    Ideally, the difference between left and right hand sides should be zero, but Mathematica tells me otherwise, as shown from the attached figure.

    Can you explain to me why as usually we would assume the elementary functions have been extremely well carried out as built-in functions and should attain the highest machine precision?

    (how to upload a figure?)
  6. Jan 30, 2016 #5


    User Avatar

    Staff: Mentor

    Below the text box, to the right, you should find a botton marked "UPLOAD".
  7. Jan 30, 2016 #6
    Thanks, I had been clicking on the figure icon on top of the textbox, and it really didn't do what I wanted. Here is the figure showing the two otherwise equivalent expressions have error bigger than machine precision and has a strong dependence on the input parameters.

    I understand the precision of the output from the built-in functions will have a dependence on the input parameters, but they should do at least what we ordinary people can achieve. Obviously if |z| is small, we should be able to make it more accurate by using taylor expansion than what gfortran comp.jpg returns us.
  8. Feb 1, 2016 #7
    One thing I might mention is that the lack of accuracy can be in NIntegrate as well as it might fail to reach prescribed precision goal I set. I have written to Mathematica to see whether there is way to detect the condition satisfied or not.
  9. Feb 1, 2016 #8
    Have you set WorkingPrecision as well as PrecisionGoal?
  10. Feb 1, 2016 #9
    I have forgotten all the FORTRAN I once knew (good riddance!) but I do remember some hard-to-track-down bugs due to not forcing double precision - particularly with complex variables.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook