Approximations for the Logarithm Function

[Thomas Klemm]

Numerical Integration of the Reciprocal Function

We use several numerical methods to calculate the integral in this formula for the logarithm:

\(
\begin{aligned}
\log(1+x) = \int_{1}^{1+x}\frac{1}{t}\,dt
\end{aligned}
\)

These approximations are then compared to its Taylor series:

\(
\begin{aligned}
\log(1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \frac{x^5}{5} - \frac{x^6}{6} + \frac{x^7}{7} + \mathcal{O}(n^8) \\
\end{aligned}
\)

Anyone familiar with the HP-41C/HP-42S will recognise the LN1+X function.

The difference in the coefficients of the first deviating term is calculated and can be used to compare the methods.

Also for the example \( x = 0.2 \) the approximation is compared with the true value:

Code:

from math import log1p

x = 0.2
L = log1p(x)
L

0.18232155679395462

Furthermore, for each approximation, a program for the HP-42S is provided with the result of the calculation for the same value.

Conclusions

Adding intervals in the case of the trapezoidal or midpoint rule brings us closer to the true coefficients of the third-order term.
However, Simpson's 1/3 rule reduces the error to a 5th order term, while its complexity is similar to the others.
We can improve this a bit with Simpson's 3/8 rule, but we can't change the order of magnitude.
Using Gaussian quadrature, the constants are not integers anymore, but for 2 points the formula is similar.
Adding one more point reduces the error to a 7th order term.

When integrating a function numerically, choose your sampling points carefully.

References

• Trapezoidal rule
• Midpoint rule
• Simpson's rule
• Gaussian quadrature
  

Trapezoidal Rule

\(
\begin{aligned}
\int_{a}^{b}f(x)\,dx
&\approx \sum_{k=1}^{N}\frac{f(x_{k-1})+f(x_{k})}{2}\Delta x_{k} \\
&= \frac{\Delta x}{2}\left(f(x_{0})+2f(x_{1})+2f(x_{2})+2f(x_{3})+2f(x_{4})+\cdots +2f(x_{N-1})+f(x_{N})\right) \\
\end{aligned}
\)

1 Interval

\(
\int_{a}^{b}f(x)\,dx \approx (b-a)\cdot {\tfrac {1}{2}}(f(a)+f(b))
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx \frac{x}{2}\,\left[1 + \frac{1}{1+x}\right] \\
&= x - \frac{x^2}{2} + \frac{x^3}{2} + \mathcal{O}(x^4) \\
\end{aligned}
\)

Difference

\(
\begin{aligned}
\frac{1}{2} - \frac{1}{3} = \frac{1}{6} \approx 0.166667
\end{aligned}
\)

Example

Code:

I = x/2*(1+1/(1+x))
I, L, (I-L)/L

(0.18333333333333335, 0.18232155679395462, 0.005549407086964257)


Program

Code:

00 { 13-Byte Prgm }
01 1
02 RCL+ ST Y
03 1/X
04 1
05 +
06 ×
07 2
08 ÷
09 END

0.183333333333
0.182321556794


2 Intervals

\(
\begin{aligned}
\int_{a}^{b}f(x)\,dx \approx \frac{1}{2}\frac{b-a}{2} \left[f(a) + 2f\left(\frac{a+b}{2}\right) + f(b)\right]
\end{aligned}
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx \frac{x}{2 \cdot 2}\,\left[1 + 2 \cdot \frac{2}{1+1+x} + \frac{1}{1+x}\right] \\
&= \frac{x}{4}\,\left[1 + \frac{4}{2+x} + \frac{1}{1+x}\right] \\
&= x - \frac{x^2}{2} + \frac{3 x^3}{8} + \mathcal{O}(x^4) \\
\end{aligned}
\)

Difference

\(
\frac{3}{8} - \frac{1}{3} = \frac{1}{24} \approx 0.0416667
\)

Example

Code:

I = x/4*(1+4/(2+x)+1/(1+x))
I, L, (I-L)/L

(0.1825757575757576, 0.18232155679395462, 0.0013942442477620518)


Program

Code:

00 { 22-Byte Prgm }
01 4
02 2
03 RCL+ ST Z
04 ÷
05 1
06 RCL+ ST Z
07 1/X
08 +
09 1
10 +
11 ×
12 4
13 ÷
14 END

0.182575757576
0.182321556794


Midpoint Rule

\(
\begin{aligned}
\int_{a}^{b}f(x)\,dx \approx \Delta x\left[f\left(a+{\tfrac {\Delta x}{2}}\right)+f\left(a+{\tfrac {3\Delta x}{2}}\right)+\ldots +f\left(b-{\tfrac {\Delta x}{2}}\right)\right]
\end{aligned}
\)

1 Interval

\(
\begin{aligned}
\int_{a}^{b}f(x)\,dx\approx (b-a)\cdot f\left(\frac{a+b}{2}\right)
\end{aligned}
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx x \, \frac{2}{2+x} \\
&= \frac{2x}{2+x} \\
&= x - \frac{x^2}{2} + \frac{x^3}{4} + \mathcal{O}(x^4) \\
\end{aligned}
\)

Difference

\(
\frac{1}{4} - \frac{1}{3} = - \frac{1}{12} \approx -0.0833333
\)

Example

Code:

I = 2*x/(2+x)
I, L, (I-L)/L

(0.18181818181818182, 0.18232155679395462, -0.0027609185914404585)


Program

Code:

00 { 9-Byte Prgm }
01 2
02 RCL+ ST Y
03 ÷
04 2
05 ×
06 END

0.181818181818
0.182321556794


2 Intervals

\(
\begin{aligned}
\int_{a}^{b}f(x)\,dx\approx \frac{b-a}{2}\,\left[f\left(\frac{3a+b}{4}\right) + f\left(\frac{a+3b}{4}\right)\right]
\end{aligned}
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx \frac{x}{2}\,\left[\frac{4}{3+1+x} + \frac{4}{1+3+3x}\right] \\
&= 2x\,\left[\frac{1}{4+x} + \frac{1}{4+3x}\right] \\
&= x - \frac{x^2}{2} + \frac{5 x^3}{16} + \mathcal{O}(x^4) \\
\end{aligned}
\)

Difference

\(
\begin{aligned}
\frac{5}{16} - \frac{1}{3} = - \frac{1}{48} \approx -0.0208333
\end{aligned}
\)

Example

Code:

I = 2*x*(1/(4+x)+1/(4+3*x))
I, L, (I-L)/L

(0.1821946169772257, 0.18232155679395462, -0.0006962414042590937)


Program

Code:

00 { 20-Byte Prgm }
01 4
02 RCL+ ST Y
03 1/X
04 3
05 RCL× ST Z
06 4
07 +
08 1/X
09 +
10 ×
11 2
12 ×
13 END

0.182194616977
0.182321556794


Simpson's Rule

1/3 rule

\(
\begin{aligned}
\int_{a}^{b}f(x)\,dx
&\approx \frac{1}{3} h \,\left[f(a)+4f\left(\frac {a+b}{2}\right)+f(b)\right]\\
&=\frac{b-a}{6}\,\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]
\end{aligned}
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx \frac{x}{6}\,\left[1 + 4 \frac{2}{1+1+x} + \frac{1}{1+x}\right] \\
&= \frac{x}{6}\,\left[1 + \frac{8}{2+x} + \frac{1}{1+x}\right] \\
&= x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \frac{5 x^5}{24} + \mathcal{O}(x^6) \\
\end{aligned}
\)

Difference

\(
\begin{aligned}
\frac{5}{24} - \frac{1}{5} = \frac{1}{120} \approx 0.00833333
\end{aligned}
\)

Example

Code:

I = x/6*(1+8/(2+x)+1/(1+x))
I, L, (I-L)/L

(0.18232323232323233, 0.18232155679395462, 9.189968027780061e-06)


Program

Code:

00 { 22-Byte Prgm }
01 8
02 2
03 RCL+ ST Z
04 ÷
05 1
06 RCL+ ST Z
07 1/X
08 +
09 1
10 +
11 ×
12 6
13 ÷
14 END

0.182323232323
0.182321556794


3/8 rule

\(
\begin{aligned}
\int _{a}^{b}f(x)\,dx
&\approx \frac{3}{8}h\,\left[f(a)+3f\left(\frac{2a+b}{3}\right)+3f\left(\frac{a+2b}{3}\right)+f(b)\right]\\
&= \frac{b-a}{8}\,\left[f(a)+3f\left(\frac{2a+b}{3}\right)+3f\left(\frac{a+2b}{3}\right)+f(b)\right] \\
\end{aligned}
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx \frac{x}{8}\,\left[1 + 3 \frac{3}{2+1+x} + 3 \frac{3}{1+2+2x} + \frac{1}{1+x}\right] \\
&= \frac{x}{8}\,\left[1 + \frac{9}{3+x} + \frac{9}{3+2x} + \frac{1}{1+x}\right] \\
&= x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \frac{11 x^5}{54} + \mathcal{O}(x^6) \\
\end{aligned}
\)

Difference

\(
\frac{11}{54} - \frac{1}{5} = \frac{1}{270} \approx 0.00370370
\)

Example

Code:

I = x/8*(1+9/(3+x)+9/(3+2*x)+1/(1+x))
I, L, (I-L)/L

(0.18232230392156862, 0.18232155679395462, 4.097856705131325e-06)


Program

Code:

00 { 33-Byte Prgm }
01 2
02 RCL× ST Y
03 3
04 +
05 1/X
06 3
07 RCL+ ST Z
08 1/X
09 +
10 9
11 ×
12 1
13 RCL+ ST Z
14 1/X
15 +
16 1
17 +
18 ×
19 8
20 ÷
21 END

0.182322303922
0.182321556794


Gauss–Legendre Quadrature

\(
\begin{aligned}
\int_{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i})
\end{aligned}
\)

2 Point

\(
\begin{array}{|c|c|}
\hline
\text{Point} & \text{Weight} \\
\hline
-\frac{1}{\sqrt{3}} & 1 \\
\hline
+\frac{1}{\sqrt{3}} & 1 \\
\hline
\end{array}
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx \frac{x}{2}\,\left[\frac{1}{1 + \left(\frac{1}{2} - \frac{1}{\sqrt{12}}\right)x} + \frac{1}{1 + \left(\frac{1}{2} + \frac{1}{\sqrt{12}}\right)x}\right] \\
&= x\,\left[\frac{1}{2 + \left(1 - \frac{1}{\sqrt{3}}\right)x} + \frac{1}{2 + \left(1 + \frac{1}{\sqrt{3}}\right)x}\right] \\
&= x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \frac{7 x^5}{36} + \mathcal{O}(x^6) \\
\end{aligned}
\)

Difference

\(
\begin{aligned}
\frac{7}{36} - \frac{1}{5} = - \frac{1}{180} \approx -0.00555556
\end{aligned}
\)

Example

Code:

from math import sqrt
k_1 = 1 - 1/sqrt(3)
k_2 = 1 + 1/sqrt(3)

I = x*(1/(2+k_1*x)+1/(2+k_2*x))
I, L, (I-L)/L

(0.18232044198895028, 0.18232155679395462, -6.114499151613408e-06)


Program

Code:

00 { 29-Byte Prgm }
01 3
02 SQRT
03 1/X
04 1
05 RCL- ST Y
06 RCL× ST Z
07 2
08 +
09 1/X
10 X<>Y
11 1
12 +
13 RCL× ST Z
14 2
15 +
16 1/X
17 +
18 ×
19 END

0.182320441989
0.182321556794

---

This variant uses registers:

1 ENTER 3 SQRT 1/X - STO 00
1 LASTX + STO 01

Code:

00 { 18-Byte Prgm }
01 RCL 00
02 RCL× ST Y
03 2
04 +
05 1/X
06 RCL 01
07 RCL× ST Z
08 2
09 +
10 1/X
11 +
12 ×
13 END

0.182320441989
0.182321556794


3 Point

\(
\begin{array}{|c|c|}
\hline
\text{Point} & \text{Weight} \\
\hline
-\sqrt{\frac{3}{5}} & \frac{5}{9} \\
\hline
0 & \frac{8}{9} \\
\hline
\sqrt{\frac{3}{5}} & \frac{5}{9} \\
\hline
\end{array}
\)

\(
\begin{aligned}
\int_{1}^{1+x}\frac{1}{t}\,dt
&\approx \frac{x}{2 \cdot 9}\,\left[\frac{5}{1 + \left(\frac{1}{2} - \sqrt{\frac{3}{20}}\right)x} + \frac{8 \cdot 2}{1+1+x} + \frac{5}{1 + \left(\frac{1}{2} + \sqrt{\frac{3}{20}}\right)x}\right] \\
&= \frac{x}{9}\,\left[\frac{5}{2 + \left(1 - \sqrt{\frac{3}{5}}\right)x} + \frac{8}{2+x} + \frac{5}{2 + \left(1 + \sqrt{\frac{3}{5}}\right)x}\right] \\
&= x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \frac{x^5}{5} - \frac{x^6}{6} + \frac{57 x^7}{400} + \mathcal{O}(x^8) \\
\end{aligned}
\)

Difference

\(
\begin{aligned}
\frac{57}{400} - \frac{1}{7} = - \frac{1}{2800} \approx -0.000357143
\end{aligned}
\)

Example

Code:

from math import sqrt
k_1 = 1 - sqrt(3/5)
k_2 = 1 + sqrt(3/5)

I = x/9*(5/(2+k_1*x)+8/(2+x)+5/(2+k_2*x))
I, L, (I-L)/L

(0.18232155441457765, 0.18232155679395462, -1.3050442389081796e-08)


Program

Code:

00 { 46-Byte Prgm }
01 0.6
02 SQRT
03 1
04 RCL- ST Y
05 RCL× ST Z
06 2
07 +
08 1/X
09 X<>Y
10 1
11 +
12 RCL× ST Z
13 2
14 +
15 1/X
16 +
17 5
18 ×
19 2
20 RCL+ ST Z
21 8
22 X<>Y
23 ÷
24 +
25 ×
26 9
27 ÷
28 END

0.182321554415
0.182321556794

---

This variant uses registers:

1 ENTER 0.6 SQRT - STO 00
1 LASTX + STO 01

Code:

00 { 34-Byte Prgm }
01 RCL 00
02 RCL× ST Y
03 2
04 +
05 1/X
06 RCL 01
07 RCL× ST Z
08 2
09 +
10 1/X
11 +
12 5
13 ×
14 2
15 RCL+ ST Z
16 8
17 X<>Y
18 ÷
19 +
20 ×
21 9
22 ÷
23 END

0.182321554415
0.182321556794

[Albert Chan]

We should compare sum with trapezoids vs squares, with same function evaluations (not intervals)

To simplify, we let t=s*x+1, to get integration limit from 0 to 1.
RHS integral value is simply averaged height of its integrand.

\(\displaystyle \ln(1+x) = \int_1^{1+x} \frac{1}{t}\,dt= \int_0^1 \frac{x}{x\,s + 1}\,ds \)

For midpoint rule with 3 function evaluation, we trisect intervals.

lua> function f(s) return x/(x*s+1) end
lua> x = 0.2
lua> log1p(x)
0.18232155679395465

Trapezoidal Rule, and squares equivalent:

lua> t1 = (f(0) + f(1))/2
lua> t2 = (t1 + f(1/2))/2
lua> t1, t2
0.18333333333333335      0.18257575757575759
lua> s1 = f(1/2)
lua> s2 = (s1 + f(1/6) + f(5/6))/3
lua> s1, s2
0.18181818181818182      0.1822650467811758

Simpson 1/3 rule, and squares equivalent:

lua> t2 + (t2-t1)/(4-1) -- weight = {1,4,1} / 6
0.18232323232323233
lua> s2 + (s2-s1)/(9-1) -- weight = {3,2,3} / 8
0.18232090490155006

Simpson 3/8 rule, and squares equivalent:

lua> require'fun'()
lua> function acc(s,x,y) return s+x*y end
lua> function dot(x,y) return foldl(acc,0,zip(x,y)) end

lua> dot({f(0/3),f(1/3),f(2/3),f(3/3)}, {1,3,3,1}) / 8
0.18232230392156865
lua> dot({f(1/8),f(3/8),f(5/8),f(7/8)}, {13,11,11,13}) / 48
0.18232121889089584

For ∫(1/x), sum with squares avoided edge bias, converge slightly better.