Leibniz vs. Newton#
Madhava series - Wikipedia (https://en.wikipedia.org/w/index.php?title=Madhava_series&oldid=1141010090)
(c.1340–c.1425) Madhava of Sangamagrama - Wikipedia
!pip install -U lxml
Collecting lxml
Using cached lxml-5.2.2-cp311-cp311-manylinux_2_28_x86_64.whl.metadata (3.4 kB)
Using cached lxml-5.2.2-cp311-cp311-manylinux_2_28_x86_64.whl (5.0 MB)
Installing collected packages: lxml
Successfully installed lxml-5.2.2
import pandas as pd
wiki = pd.read_html('https://en.wikipedia.org/w/index.php?title=Madhava_series&oldid=1141010090')
wiki[0]
| No. | Series | Name | Western discoverers of the series and approximate dates of discovery[6] | |
|---|---|---|---|---|
| 0 | 1 | sin x = x − x3/3! + x5/5! − x7/7! + ... | Madhava's sine series | Isaac Newton (1669) and Wilhelm Leibniz (1676) |
| 1 | 2 | cos x = 1 − x2/2! + x4/4! − x6/6! + ... | Madhava's cosine series | Isaac Newton (1669) and Wilhelm Leibniz (1676) |
| 2 | 3 | arctan x = x − x3/3 + x5/5 − x7/7 + ... | Madhava's series for arctangent | James Gregory (1671) and Wilhelm Leibniz (1673) |
| 3 | 4 | π/4 = 1 − 1/3 + 1/5 − 1/7 + ... | Madhava's formula for π | James Gregory (1671) and Wilhelm Leibniz (1673) |
SymPyで検証する#
from sympy import *
from sympy.abc import x
series(sin(x), x, 0, 10)
\[\displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \frac{x^{9}}{362880} + O\left(x^{10}\right)\]
series(cos(x), x, 0, 10)
\[\displaystyle 1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} - \frac{x^{6}}{720} + \frac{x^{8}}{40320} + O\left(x^{10}\right)\]
series(atan(x), x, 0, 10)
\[\displaystyle x - \frac{x^{3}}{3} + \frac{x^{5}}{5} - \frac{x^{7}}{7} + \frac{x^{9}}{9} + O\left(x^{10}\right)\]
float(series(atan(x), x, 0, 10).removeO().subs(x, 1))*4
3.33968253968254
float(series(atan(x), x, 0, 20).removeO().subs(x, 1))*4
3.0418396189294024
float(series(atan(x), x, 0, 40).removeO().subs(x, 1))*4
3.0916238066678385
float(series(atan(x), x, 0, 80).removeO().subs(x, 1))*4
3.1165965567938323
float(series(atan(x), x, 0, 1000).removeO().subs(x, 1))*4
3.1395926555897833
float(series(atan(x), x, 0, 2000).removeO().subs(x, 1))*4
3.140592653839793
Newtonによる \(\sin(x)\) の級数展開#
\(x = \sin(y)\) の逆関数 \(y = \arcsin(x)\) を級数展開する#
\(\sin(y) = x\) とおく。\(x\) を \(x\) 軸上の点とすると、\(y = \sin^{-1}(x) = \arcsin(x)\) は、次の図の円弧の長さに対応する。
単位円の周長 \(l\) と面積 \(s\) は、\(l : 2 \pi = s : \pi\) より \(\frac{l}{2} = s\) の関係がある。
\[
\frac{1}{2} y = \frac{1}{2}\sin^{-1}(x) = \int_{0}^{x} \sqrt{1-x^2} - \frac{1}{2} x \sqrt{1-x^2}
\]
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import numpy as np
x = np.cos(np.radians(30))
y = np.sin(np.radians(30))
fig, ax = plt.subplots()
arc1 = patches.Arc([0,0], 2,2, theta1=0, theta2=30, linestyle=':')
arc2 = patches.Arc([0,0], 2,2, theta1=30, theta2=90)
wdg2 = patches.Wedge([0,0], 1, theta1=30, theta2=90, fill=True, alpha=0.3)
ax.add_patch(arc1)
ax.add_patch(arc2)
ax.add_patch(wdg2)
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.set_aspect("equal")
plt.plot(x, 0, "ok")
plt.plot(x, y, "ok")
plt.plot([0,x],[y,y], 'k:')
plt.plot([x,x],[0,y], 'k:')
plt.text(x/2, y/2, '1')
plt.text(x+.02, .04, 'x')
plt.text(x-.25, .9, 'y = arcsin(x)')
plt.text(x-.25, .2, 'arccos(x)')
plt.show()
-
“Isaac Newton is generally credited with the generalized binomial theorem, valid for any rational exponent.”
「アイザック・ニュートンは、一般に、任意の有理数の指数に対して有効な一般二項定理を導いたとされている。」
\[
(x + y)^r = \sum_{k=0}^\infty {r \choose k} x^{r-k} y^k
\]
\(\sqrt{1-x^2} = (1-x^2)^\frac{1}{2}\) を二項級数に展開する:
from sympy.abc import x
x
\[\displaystyle x\]
(1-x**2)**(1/2)
\[\displaystyle \left(1 - x^{2}\right)^{0.5}\]
nsimplify(_)
\[\displaystyle \sqrt{1 - x^{2}}\]
series((1-x**2)**(1/2))
\[\displaystyle 1 - 0.5 x^{2} - 0.125 x^{4} + O\left(x^{6}\right)\]
nsimplify(series((1-x**2)**(1/2), x, 0, 10))
\[\displaystyle 1 - \frac{x^{2}}{2} - \frac{x^{4}}{8} - \frac{x^{6}}{16} - \frac{5 x^{8}}{128} + O\left(x^{10}\right)\]
積分する。多項式なので項ごとに積分できる:
nsimplify(integrate(series((1-x**2)**(1/2), x, 0, 10), x))
\[\displaystyle x - \frac{x^{3}}{6} - \frac{x^{5}}{40} - \frac{x^{7}}{112} - \frac{5 x^{9}}{1152} + O\left(x^{11}\right)\]
これらによって、\(\sin(x)\) の逆関数 \(\sin^{-1}(x) = \arcsin(x)\) の級数展開はすぐに求められる:
arcsin_x = 2 * integrate(series((1-x**2)**(1/2), x, 0, 10), x) - expand(x * series((1-x**2)**(1/2), x, 0, 10))
nsimplify(arcsin_x)
\[\displaystyle x + \frac{x^{3}}{6} + \frac{3 x^{5}}{40} + \frac{5 x^{7}}{112} + \frac{35 x^{9}}{1152} + O\left(x^{11}\right)\]
\(0.1\) と \(1\) のときの値を求めて検算する:
arcsin_x.removeO().subs(x,.1), asin(.1)
(0.100167421161334, 0.100167421161560)
arcsin_x.removeO().subs(x,1.), asin(1.)
(1.31669146825397, 1.57079632679490)
import matplotlib.pyplot as plt
import numpy as np
X = np.linspace(-1, 1, 100)
Y = [arcsin_x.removeO().subs(x,xi) for xi in X]
plt.plot(X, Y)
plt.plot(X, np.arcsin(X))
[<matplotlib.lines.Line2D at 0x7f89609d0370>]
SymPyのseries()による検算:
series(asin(x),x, 0, 15)
\[\displaystyle x + \frac{x^{3}}{6} + \frac{3 x^{5}}{40} + \frac{5 x^{7}}{112} + \frac{35 x^{9}}{1152} + \frac{63 x^{11}}{2816} + \frac{231 x^{13}}{13312} + O\left(x^{15}\right)\]
Newtonの力業#
恒等式 \(\sin^{-1}(\sin(x)) = x\) がなりたつ。さらに、\(\sin(x)\) の級数展開を \(\sin(x) = a_{0} + a_{1} x + a_{2} x^2 + a_{3} x^3 + \cdots\) とおいて、上で得られた \(\sin^{-1}(x)\) の級数展開式へ代入しする。変数 \(x\) の係数が、一次のみ \(1\)、それ以外の次数では \(0\) となることから、順次 \(a_{0}, a_{1}, a_{2},...\) を求める。
\[\begin{split}
\begin{align}
\sin(x) &\approx a_{0} \\
\sin(x) &\approx a_{0} + a_{1} x \\
\sin(x) &\approx a_{0} + a_{1} x + a_{2} x^2 \\
\sin(x) &\approx a_{0} + a_{1} x + a_{2} x^2 + a_{3} x^3 \\
\end{align}
\end{split}\]
a=symbols('a:11', Real=True)
sin_x=a[0]+a[1]*x+a[2]*x**2+a[3]*x**3+a[4]*x**4+a[5]*x**5+a[6]*x**6+a[7]*x**7+a[8]*x**8+a[9]*x**9
sol=[[a[0],0],[a[1],0],[a[2],0],[a[3],0],[a[4],0],[a[5],0],[a[6],0],[a[7],0],[a[8],0],[a[9],0]]
sol
[[a0, 0],
[a1, 0],
[a2, 0],
[a3, 0],
[a4, 0],
[a5, 0],
[a6, 0],
[a7, 0],
[a8, 0],
[a9, 0]]
sin_x.subs(sol)
\[\displaystyle 0\]
ax=Symbol('ax', real=True)
sol[1][1]=ax
sol
[[a0, 0],
[a1, ax],
[a2, 0],
[a3, 0],
[a4, 0],
[a5, 0],
[a6, 0],
[a7, 0],
[a8, 0],
[a9, 0]]
sin_x.subs(sol)
\[\displaystyle ax x\]
bx=nsimplify(expand(arcsin_x.removeO().subs(x, sin_x.subs(sol))))
bx
\[\displaystyle \frac{35 ax^{9} x^{9}}{1152} + \frac{5 ax^{7} x^{7}}{112} + \frac{3 ax^{5} x^{5}}{40} + \frac{ax^{3} x^{3}}{6} + ax x\]
恒等式 \(\sin^{-1}(\sin(x)) = x\) より、多項式の係数を比較すると、恒等式が成立する条件は \(a_{x} = 1\) のとき。
cx=poly(bx-x,x).all_coeffs()[::-1]
cx[1]
\[\displaystyle ax - 1\]
解析的に解く#
solve(cx[1],ax)
[1]
このように、\(a_{0}\) から順に係数 \(a_1\), \(a_2\), \(a_3\)… を決めていくプログラムは次のようになる:
sol=[[a[0],0],[a[1],0],[a[2],0],[a[3],0],[a[4],0],[a[5],0],[a[6],0],[a[7],0],[a[8],0],[a[9],0]]
for i in range(10):
sol[i][1]=ax
print("sin(x) =", sin_x.subs(sol))
bx=nsimplify(expand(arcsin_x.removeO().subs(x, sin_x.subs(sol))))
cx=poly(bx-x,x).all_coeffs()[::-1]
s=solve(cx[i],ax)
print("a[{:d}] =".format(i), s[0])
sol[i][1]=s[0]
nsimplify(sin_x.subs(sol))
sin(x) = ax
a[0] = 0
sin(x) = ax*x
a[1] = 1
sin(x) = ax*x**2 + x
a[2] = 0
sin(x) = ax*x**3 + x
a[3] = -1/6
sin(x) = ax*x**4 - x**3/6 + x
a[4] = 0
sin(x) = ax*x**5 - x**3/6 + x
a[5] = 1/120
sin(x) = ax*x**6 + x**5/120 - x**3/6 + x
a[6] = 0
sin(x) = ax*x**7 + x**5/120 - x**3/6 + x
a[7] = -39682539682541/200000000000000000
sin(x) = ax*x**8 - 39682539682541*x**7/200000000000000000 + x**5/120 - x**3/6 + x
a[8] = 0
sin(x) = ax*x**9 - 39682539682541*x**7/200000000000000000 + x**5/120 - x**3/6 + x
a[9] = 137786596120411/50000000000000000000
\[\displaystyle \frac{137786596120411 x^{9}}{50000000000000000000} - \frac{39682539682541 x^{7}}{200000000000000000} + \frac{x^{5}}{120} - \frac{x^{3}}{6} + x\]
1/(39682539682541/200000000000000000)
5039.999999999833
1/(137786596120411/50000000000000000000)
362879.9999987318
import math
list(map(math.factorial, range(11)))
[1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800]
# SymPyのseries()を呼び出して検算
series(sin(x),x, 0, 10)
\[\displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \frac{x^{9}}{362880} + O\left(x^{10}\right)\]
# グラフプロットで可視化
import matplotlib.pyplot as plt
import numpy as np
X = np.linspace(-2*np.pi, 2*np.pi, 100)
Y = [sin_x.subs(sol).subs(x,xi) for xi in X]
plt.plot(X, Y)
plt.plot(X, np.sin(X))
[<matplotlib.lines.Line2D at 0x7f89608eca30>]
力業の続き: 代数的に解く#
以下、力業で解く場合の続き:
定数項は \(0\) なので、\(a_{0} = 0\)
nsimplify(expand(asx.removeO().subs(x, a[1]*x)))
\[\displaystyle \frac{35 a_{1}^{9} x^{9}}{1152} + \frac{5 a_{1}^{7} x^{7}}{112} + \frac{3 a_{1}^{5} x^{5}}{40} + \frac{a_{1}^{3} x^{3}}{6} + a_{1} x\]
poly(nsimplify(expand(asx.removeO().subs(x, a[1]*x))),x).all_coeffs()[::-1][1]
\[\displaystyle a_{1}\]
\(x\) の係数は \(1\) なので、\(a_{1} = 1\)
nsimplify(expand(asx.removeO().subs(x, x+a[2]*x**2)))
\[\displaystyle \frac{35 a_{2}^{9} x^{18}}{1152} + \frac{35 a_{2}^{8} x^{17}}{128} + \frac{35 a_{2}^{7} x^{16}}{32} + \frac{5 a_{2}^{7} x^{14}}{112} + \frac{245 a_{2}^{6} x^{15}}{96} + \frac{5 a_{2}^{6} x^{13}}{16} + \frac{245 a_{2}^{5} x^{14}}{64} + \frac{15 a_{2}^{5} x^{12}}{16} + \frac{3 a_{2}^{5} x^{10}}{40} + \frac{245 a_{2}^{4} x^{13}}{64} + \frac{25 a_{2}^{4} x^{11}}{16} + \frac{3 a_{2}^{4} x^{9}}{8} + \frac{245 a_{2}^{3} x^{12}}{96} + \frac{25 a_{2}^{3} x^{10}}{16} + \frac{3 a_{2}^{3} x^{8}}{4} + \frac{a_{2}^{3} x^{6}}{6} + \frac{35 a_{2}^{2} x^{11}}{32} + \frac{15 a_{2}^{2} x^{9}}{16} + \frac{3 a_{2}^{2} x^{7}}{4} + \frac{a_{2}^{2} x^{5}}{2} + \frac{35 a_{2} x^{10}}{128} + \frac{5 a_{2} x^{8}}{16} + \frac{3 a_{2} x^{6}}{8} + \frac{a_{2} x^{4}}{2} + a_{2} x^{2} + \frac{35 x^{9}}{1152} + \frac{5 x^{7}}{112} + \frac{3 x^{5}}{40} + \frac{x^{3}}{6} + x\]
poly(nsimplify(expand(asx.removeO().subs(x, x+a[2]*x**2))),x).all_coeffs()[::-1][2]
\[\displaystyle a_{2}\]
\(a_{2} = 0\)
nsimplify(expand(asx.removeO().subs(x, x+a[3]*x**3)))
\[\displaystyle \frac{35 a_{3}^{9} x^{27}}{1152} + \frac{35 a_{3}^{8} x^{25}}{128} + \frac{35 a_{3}^{7} x^{23}}{32} + \frac{5 a_{3}^{7} x^{21}}{112} + \frac{245 a_{3}^{6} x^{21}}{96} + \frac{5 a_{3}^{6} x^{19}}{16} + \frac{245 a_{3}^{5} x^{19}}{64} + \frac{15 a_{3}^{5} x^{17}}{16} + \frac{3 a_{3}^{5} x^{15}}{40} + \frac{245 a_{3}^{4} x^{17}}{64} + \frac{25 a_{3}^{4} x^{15}}{16} + \frac{3 a_{3}^{4} x^{13}}{8} + \frac{245 a_{3}^{3} x^{15}}{96} + \frac{25 a_{3}^{3} x^{13}}{16} + \frac{3 a_{3}^{3} x^{11}}{4} + \frac{a_{3}^{3} x^{9}}{6} + \frac{35 a_{3}^{2} x^{13}}{32} + \frac{15 a_{3}^{2} x^{11}}{16} + \frac{3 a_{3}^{2} x^{9}}{4} + \frac{a_{3}^{2} x^{7}}{2} + \frac{35 a_{3} x^{11}}{128} + \frac{5 a_{3} x^{9}}{16} + \frac{3 a_{3} x^{7}}{8} + \frac{a_{3} x^{5}}{2} + a_{3} x^{3} + \frac{35 x^{9}}{1152} + \frac{5 x^{7}}{112} + \frac{3 x^{5}}{40} + \frac{x^{3}}{6} + x\]
poly(nsimplify(expand(asx.removeO().subs(x, x+a[3]*x**3))),x).all_coeffs()[::-1]
[0,
1,
0,
a3 + 1/6,
0,
a3/2 + 3/40,
0,
a3**2/2 + 3*a3/8 + 5/112,
0,
a3**3/6 + 3*a3**2/4 + 5*a3/16 + 35/1152,
0,
3*a3**3/4 + 15*a3**2/16 + 35*a3/128,
0,
3*a3**4/8 + 25*a3**3/16 + 35*a3**2/32,
0,
3*a3**5/40 + 25*a3**4/16 + 245*a3**3/96,
0,
15*a3**5/16 + 245*a3**4/64,
0,
5*a3**6/16 + 245*a3**5/64,
0,
5*a3**7/112 + 245*a3**6/96,
0,
35*a3**7/32,
0,
35*a3**8/128,
0,
35*a3**9/1152]
poly(nsimplify(expand(asx.removeO().subs(x, x+a[3]*x**3))),x).all_coeffs()[::-1][3]
\[\displaystyle a_{3} + \frac{1}{6}\]
\(a_{3} + \frac{1}{6} = 0\) より、\(a_{3} = - \frac{1}{6}\)
poly(nsimplify(expand(asx.removeO().subs(x, x-1/6*x**3+a[5]*x**5))),x).all_coeffs()[::-1][5]
\[\displaystyle a_{5} - \frac{1}{120}\]
\(a_{5} - \frac{1}{120} = 0\) より、\(a_{5} = \frac{1}{120}\)
\(a_7\) は計算精度がたりない・・・
poly(nsimplify(expand(asx.removeO().subs(x, x-1/6*x**3+1/120*x**5+a[7]*x**7))),x).all_coeffs()[::-1][7]
\[\displaystyle a_{7} + \frac{39682539682541}{200000000000000000}\]
200000000000000000/39682539682541
5039.999999999833