Approximating logarithm using harmonic mean

Question

Here is a function to approximate log10(x+1) for (x+1) < ~1.2:

a = 1.097
b = 0.085
c = 2.31
ans = 1 / (a - b*x + c/x)

It should look like that:

equation of approximation

It works by adjusting harmonic mean to match log10, but the problem is in values of a, b, c.

The question is how to get just right a, b and c and how to make better approximation.

I made this code that can give a pretty good approximation for a, b, c, but my code wasn't able to make it any better.

import numpy as np
a = 1
b = 0.01
c = 2
def mlg(t):
 x = t
 if t == 0:
 x = 0.00000001
 x2 = x*x
 o = a - (b * x) + (c / x)
 return 1/o
def mlg0(t):
 x = t
 if t == 0:
 x = 0.00000001
 x2 = x*x
 o = a - (b * x) + (c / x)
 return o
for i in range(9000):
 n1 = np.random.uniform(0,1.19,1000)
 for i in range(1000):
 n = n1[i]
 o = np.log10(n+1)
 u = mlg(n) - o
 e = u ** 2
 de_da = 0 - 2 * (u) / (mlg0(n) ** 2)
 de_db = de_da * n 
 de_dc = de_da / n
 a -= de_da * 0.00001
 b -= de_db * 0.00001
 c -= de_dc * 0.00001
print(a,b,c)

How could the code be changed to generate better values?

I've used a method alike back propagation in NN, but it could not give me values any better.

Here is how the error is calculated:

equation of error calculation

• mkrieger1: By what metric do you consider values for a, b, c "better"?
• mkrieger1: "My code wasn't able to make it better" - what makes you think that it could be better?
• Eternal Dreamer: For what purpose did you set the first loop ? To range 9000. If it's genuinely useful you should use another for your variable i. I also don't understand why you subtract a, b and c this way, although I'm not into math like you. Could you detail a little more ? Maybe add comments to your code
• luaenjoyer: @Mike there is no purpose of "i" so there is no problem of using it 2 times. i subtracted b just because the true value of b in standart form (ax^2 +bx +c) is -0.085, so thats the same as + (-0.085)
• luaenjoyer: @mkrieger1 i was talking about "better approximation on (0<x<1) " so thats integral from 0 to 1 of | approx(x) - log10(x) | dx

Answer 1

Here are two approaches.

Method 1: series expansion in x (better for negative and positive x)

Method 2: fit the curve that passes through 3 points (here, x=0, 1⁄2 and 1)

Method 1.

If you expand them by Taylor series as powers of x then

enter image description here

Equating coefficients of x, x^2 and x^3 gives

enter image description here

In code:

import math
import numpy as np
import matplotlib.pyplot as plt
c = math.log( 10 )
a = c / 2
b = c / 12
print( "a, b, c = ", a, b, c )
x = np.linspace( -0.2, 0.2, 50 )
y = np.log10( 1 + x )
fit = x / ( a * x - b * x ** 2 + c )
plt.plot( x, y , 'b-', label="Original" )
plt.plot( x, fit, 'ro', label="Fit" )
plt.legend()
plt.show()

Output:

a, b, c = 1.151292546497023 0.19188209108283716 2.302585092994046

enter image description here

Method 2.

Fit to three points. Here we require

enter image description here

If we require this to fit at x=0, 1⁄2 and 1 we get (including L’Hopital’s rule for the limit at x=0)

enter image description here

This time I have used your interval x in [0,1] to plot the fit

import math
import numpy as np
import matplotlib.pyplot as plt
c = math.log( 10 )
a = c * ( 2/math.log(1.5) - 1/math.log(2) - 3 )
b = 2 * c * ( 1/math.log(1.5) - 1/math.log(2) - 1 )
print( "a, b, c = ", a, b, c )
x = np.linspace( 0.0, 1.0, 50 )
y = np.log10( 1 + x )
fit = x / ( a * x - b * x ** 2 + c )
plt.plot( x, y , 'b-', label="Original" )
plt.plot( x, fit, 'ro', label="Fit" )
plt.legend()
plt.show()

Output:

a, b, c = 1.1280638006656465 0.1087207987723298 2.302585092994046

enter image description here

• luaenjoyer: In case of your a, b, c total error ( integral from 0 to 1 of | approx_function(x) - log10(x) | dx ) is 0.00130902544352. In first case when a is 1.097, b is 0.085, and c is 2.31 - total error is 0.0000687849876695. the question is how to get better abc or how even starting abc`s were calculated. i used area under the curve (approx - true) to calculate error here.
• lastchance: The level of fit depends on your error metric. You are expanding in x, so it would be more natural to take a symmetric interval. What happens if you compute integral errors from x=-0.5 to 0.5 rather than 0 to 1? Also, integral errors are more usually evaluated with the 2-norm than the abs() value here: it behaves better analytically.
• luaenjoyer: i see, for smaller x your functions behave way better, thanks a lot
• lastchance: I have included a second method, @luaenjoyer. The second method simply fits the required parametric curve through three points. To accommodate your interval [0,1] I have chosen the fit points at 0, 1/2, 1 (the value x=0 requiring a limit taken from L'Hopital's rule or equivalent).
• luaenjoyer: 0.0000537050434564 total error (in case of second method) is ~22% better result than starting abc`s, man thats amazing. I am not that good at math so my best result using something like back propagation but without neural network is 0.000291286020986. By the way do you as university lecturer think this is "good" task for students? Just asking

Answer 2

The Taylor series is almost never the right way to do it unless the series that you obtain is very well behaved and strongly convergent (eg exp(x)). Log(1+x) is neither with alternating very slowly decreasing coefficients.

Since log10 is just natural log with a scaling factor I'll treat that.

log10(1+x) = log(1+x)/log(10)

The optimum approximation for most functions is typically a rational polynomial. The Pade approximants to match the first 3 terms of the Taylor series for log(1+x) are starting from the original Taylor series:

[3,0] x-x^2/2 + x^3/3
[2,1] x(6+x)/(6x+4)
[1,2] 12x/(12+6x-x^2) = 1/(1/x+1/2-x/12)

Of these options (OP's is equivalent to the last one) give or take the scale factor of log(10) the middle expression is by far the best simple rational Pade approximation for log(1+x). Both of the equations with a non trivial divisor benefit from extending the equivalent polynomial expansion with more terms and extending the series convergence . It is worth looking at the raw numbers:

Function	x=-0.5	x = 1
log(1+x)	-0.693147	0.693147
Poly[3,0]	-0.66667	0.83333
Pade[2,1]	-0.6875	0.7
Pade[1,2]	-0.68571	0.70588

It is often the case that the rational polynomial with about the same number of terms in numerator and denominator provides the best compromise fit. It can more accurately model both poles and zeroes.

When optimised for equal absolute ripple error of 0.00024 on the range [0,1] the optimal coefficients are slightly different to the analytical result:

5.998550 0.698519 6 3.658462

An equal relative error approximation is also possible on the same range

6 0.780703 6 3.780244

It is also possible to optimise it for a region either side of zero which will obtain the best possible accuracy with a simple formula. In the OP's case -0.2 to +0.2 gives a coefficient set for [2,1] with minimax ripple 5.8e-5

 5.999871522 1.019979324 6 4.024891596
 1.490691459 0.253417837 1.490723379 1

The second set are in canonical form to avoid one multiplication. There are tables of these approximations dating way back. Hart, Computer Approximations 1968 is one such bible and although dated still useful.

An even better one is possible by considering the series for log((1-y)/(1+y)) where y = x/(2+x). It all depends how much accuracy you need.

log(1+x) = log((1+y)/(1-y)) ~= 2y(15 - 4y^2 )/(15 - 9y^2)

which matches the Taylor expansion 2y + 2y^3/3 + 2y^5/5 + ... an already better behaved series to begin with.

• luaenjoyer: Man, I'm only in 10th grade, it is very interesting but i didnt understand anything. I got too many questions. Could you tell me what to do with that 4 coefficients you calculated. Its like i`m just too stupid for that
• Martin Brown: They go where 6,1,6,4 were in the original [2,1] Pade approximation (and that trick itself is typically used at degree level and above). IOW for your original |x| <= 0.2 approx you want x*(5.999871522 +1.019979324*x)/(6 + 4.024891596*x). That's for natural logarithm base e (divide by log(10) = 2.302585 to get base 10). What you did originally was find one of those Pade forms [1,2] which gives a much better answer than the more obvious Taylor series expansion. [3,0] So well done! [N,D] is shorthand for highest power of x in the numerator and denominator respectively.
• Stef: OP appears to be optimising their approximation of log10(1+x) over [0,1], so there is an argument for using the Taylor series at x=0.5 instead, ie using the Taylor series of log10(1.5+u) at u=0.
• Martin Brown: @Stef they specified the range of 1-1.2 on the first line (their derived approximation is good for 1 +/- 0.2). What you suggest will work for any of these approximations since `log(1.5+u) = log(1.5)+log(1+u/1.5). The Taylor series for log(1+x) has just about the worst convergence possible (and diverges for x=-1). Rational approximations always win out for alternating series.