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