5.3 用 SymPy 展開級數並估計餘項#

import sympy as sp
import inspect

print("sp.series 簽名:", inspect.signature(sp.series))
print("Expr.nseries 簽名:", inspect.signature(sp.exp(sp.Symbol('x')).nseries))

sp.series 簽名: (expr, x=None, x0=0, n=6, dir='+')
Expr.nseries 簽名: (x=None, x0=0, n=6, dir='+', logx=None, cdir=0)

import sympy as sp

x = sp.symbols('x', real=True)

# sp.series: 符號展開（回傳 Add 型別，內含 O 物件）
s = sp.series(sp.exp(x), x, 0, 6)
print("exp(x) 的 series =", s)
print("型別 =", type(s))
print("O 物件 =", s.getO())
print("removeO 後 =", s.removeO())
print("係數字典 =", dict(s.as_coefficients_dict()))

# Expr.nseries: 快速數值展開（方法的 dir 參數）
print()
print("exp(x).nseries(x, 0, 6) =", sp.exp(x).nseries(x, 0, 6))
print("數值代入 x=1 =", sp.N(sp.exp(x).nseries(x, 0, 6).removeO().subs(x, 1), 15))
print("exp(1) 真值 =", sp.N(sp.exp(1), 15))

exp(x) 的 series = 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
型別 = <class 'sympy.core.add.Add'>
O 物件 = O(x**6)
removeO 後 = x**5/120 + x**4/24 + x**3/6 + x**2/2 + x + 1
係數字典 = {1: 1, x: 1, x**2: 1/2, x**3: 1/6, x**4: 1/24, x**5: 1/120, O(x**6): 1}

exp(x).nseries(x, 0, 6) = 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
數值代入 x=1 = 2.71666666666667
exp(1) 真值 = 2.71828182845905

import sympy as sp

x = sp.symbols('x', real=True)

# Order 物件：SymPy 如何表達「O(x^n)」
O6 = sp.Order(x**6)
print("O 物件 =", O6, "  型別 =", type(O6))
print("O(x^6) + O(x^6) =", O6 + O6)
print("O(x^6) + O(x^3) =", O6 + sp.Order(x**3))
print("O(x^6) * x^2 =", O6 * x**2)
print("O(x^6) / x^2 =", O6 / x**2)
print("as_leading_term =", O6.as_leading_term(x))

# series 物件裡如何取出 O
print()
s = sp.series(sp.exp(x), x, 0, 4)
print("series =", s)
print("getO() =", s.getO())
print("removeO() =", s.removeO())
print("O 物件的 args =", s.getO().args)

# 比較兩個 O 的階數：取「較低階」者
print()
print("O(x^6) - O(x^3) =", sp.Order(x**6) - sp.Order(x**3))
print("O(x^3) - O(x**6) =", sp.Order(x**3) - sp.Order(x**6))
print("O(x^3).contains(x^3)? ", sp.Order(x**3).contains(x**3))
print("O(x^6).contains(x^3)? ", sp.Order(x**6).contains(x**3))

# dir 參數：log(x) 在 0 處的左右展開不同
print()
print("log(x) dir='+' =", sp.series(sp.log(x), x, 0, 3, dir='+'))
print("log(x) dir='-' =", sp.series(sp.log(x), x, 0, 3, dir='-'))

# 在非零展開點也可以展開
print()
print("sin(x) 在 x=pi 展開 =", sp.series(sp.sin(x), x, sp.pi, 4))

# nseries 的 dir 參數
print()
print("exp(x).nseries(x, 0, 5, dir='-') =", sp.exp(x).nseries(x, 0, 5, dir='-'))

O 物件 = O(x**6)   型別 = <class 'sympy.series.order.Order'>
O(x^6) + O(x^6) = O(x**6)
O(x^6) + O(x^3) = O(x**3)
O(x^6) * x^2 = O(x**8)
O(x^6) / x^2 = O(x**4)
as_leading_term = O(x**6)

series = 1 + x + x**2/2 + x**3/6 + O(x**4)
getO() = O(x**4)
removeO() = x**3/6 + x**2/2 + x + 1
O 物件的 args = (x**4, (x, 0))

O(x^6) - O(x^3) = O(x**3)
O(x^3) - O(x**6) = O(x**3)
O(x^3).contains(x^3)?  True
O(x^6).contains(x^3)?  False

log(x) dir='+' = log(x)
log(x) dir='-' = log(-x) + I*pi

sin(x) 在 x=pi 展開 = pi + (x - pi)**3/6 - x + O((x - pi)**4, (x, pi))

exp(x).nseries(x, 0, 5, dir='-') = 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

import sympy as sp

x, xi = sp.symbols('x xi', real=True)

# 驗證 Lagrange 餘項：R_n = f^(n+1)(xi)/(n+1)! * (x-a)^(n+1)
f = sp.exp(x)
a, b = sp.Integer(0), sp.Integer(1)

for n in [1, 2, 3]:
    T = sum(sp.diff(f, x, k).subs(x, a) * (x - a)**k / sp.factorial(k) for k in range(n+1))
    R = sp.simplify(f.subs(x, b) - T.subs(x, b))
    expr = sp.diff(f, x, n+1).subs(x, xi) / sp.factorial(n+1) * (b - a)**(n+1)
    sol = sp.solve(sp.Eq(expr, R), xi)
    print(f"n = {n}: T_n = {sp.expand(T)}")
    print(f"      R_n = {R}")
    for s_ in sol:
        ok = sp.simplify(a < s_ < b)
        print(f"      Lagrange 的 xi = {s_}  (在 (a,b) 內? {ok}, 數值 {float(s_):.6f})")

print()
print("=== 數值核對：sin 在 pi/2 的餘項 ===")
g = sp.sin(x)
a2 = sp.Integer(0)
b2 = sp.pi/2
n2 = 3
T2 = sum(sp.diff(g, x, k).subs(x, a2) * (x - a2)**k / sp.factorial(k) for k in range(n2+1))
R2 = sp.simplify(g.subs(x, b2) - T2.subs(x, b2))
print("T_3 =", sp.expand(T2), " R_3 =", R2)
expr2 = sp.diff(g, x, n2+1).subs(x, xi) / sp.factorial(n2+1) * (b2 - a2)**(n2+1)
print("xi 的方程 =", sp.Eq(expr2, R2))
print("數值解 xi ≈", sp.nsolve(expr2 - R2, xi, 0.5))

n = 1: T_n = x + 1
      R_n = -2 + E
      Lagrange 的 xi = log(-4 + 2*E)  (在 (a,b) 內? True, 數值 0.362254)
n = 2: T_n = x**2/2 + x + 1
      R_n = -5/2 + E
      Lagrange 的 xi = log(-15 + 6*E)  (在 (a,b) 內? True, 數值 0.269791)
n = 3: T_n = x**3/6 + x**2/2 + x + 1
      R_n = -8/3 + E
      Lagrange 的 xi = log(-64 + 24*E)  (在 (a,b) 內? True, 數值 0.214114)

=== 數值核對：sin 在 pi/2 的餘項 ===
T_3 = -x**3/6 + x  R_3 = -pi/2 + pi**3/48 + 1
xi 的方程 = Eq(pi**4*sin(xi)/384, -pi/2 + pi**3/48 + 1)
數值解 xi ≈ 0.300839035833927

import sympy as sp
import numpy as np

x = sp.symbols('x', real=True)

# 通用工具：對 f 求「在 x=a、要誤差 < eps、需要幾項」
def needed_terms(f, a, h, eps, samples=201):
    """回傳使 Lagrange 餘項界 < eps 的最小 n（M 由區間取樣估計）"""
    grid = np.linspace(float(a), float(a + h), samples)
    for n in range(0, 40):
        d = sp.diff(f, x, n + 1)
        M = max(abs(complex(sp.N(d.subs(x, t)))) for t in grid)
        bound = M / sp.factorial(n + 1) * abs(float(h)) ** (n + 1)
        if bound < eps:
            return n, bound
    return None, None

for (name, f) in [("e^1", sp.exp(x)), ("sin(1)", sp.sin(x)), ("cos(1)", sp.cos(x)),
                  ("ln(2)", sp.log(1 + x)), ("arctan(1)", sp.atan(x))]:
    for eps in [1e-6, 1e-10, 1e-15]:
        n, bound = needed_terms(f, 0, 1, eps)
        if n is None:
            print(f"{name:10s} 目標 {eps:.0e}: Lagrange 界在 40 階內都不夠（改用交錯級數界）")
        else:
            print(f"{name:10s} 目標 {eps:.0e}: 需要 n = {n} (項數 {n+1}), 餘項界 = {bound:.3e}")
    print()

# ln(2) 的 Maclaurin 級數 1 - 1/2 + 1/3 - ... 收斂極慢：
# 交錯級數界是 1/(n+1)，要達到 eps 需要約 1/eps 項
print("ln(2) = 1 - 1/2 + 1/3 - 1/4 + ...（交錯級數）")
for eps in [1e-6, 1e-10, 1e-15]:
    n = int(1/eps) + 1
    print(f"  目標 {eps:.0e}: 交錯級數界要求 n >= 1/eps ≈ {n} 項")

# 改寫：ln 2 = 2 arctanh(1/3) = 2 * sum (1/3)^(2k+1)/(2k+1)，收斂快得多
print("\n改寫 ln 2 = 2 arctanh(1/3) = 2 * sum_{k>=0} (1/3)^(2k+1)/(2k+1)")
for eps in [1e-6, 1e-10, 1e-15]:
    k = 0
    while float(2 * sp.Rational(1, 3)**(2*k+3)/(2*k+3)) >= eps:
        k += 1
    partial = 2 * sum(sp.Rational(1, 3)**(2*m+1)/(2*m+1) for m in range(k+1))
    print(f"  目標 {eps:.0e}: 只要 {k+1} 項，和 = {float(partial):.15f}   誤差 = {abs(float(partial)-float(sp.log(2))):.3e}")

e^1        目標 1e-06: 需要 n = 9 (項數 10), 餘項界 = 7.491e-7
e^1        目標 1e-10: 需要 n = 13 (項數 14), 餘項界 = 3.118e-11
e^1        目標 1e-15: 需要 n = 17 (項數 18), 餘項界 = 4.246e-16

sin(1)     目標 1e-06: 需要 n = 9 (項數 10), 餘項界 = 2.319e-7
sin(1)     目標 1e-10: 需要 n = 13 (項數 14), 餘項界 = 9.652e-12
sin(1)     目標 1e-15: 需要 n = 17 (項數 18), 餘項界 = 1.314e-16

cos(1)     目標 1e-06: 需要 n = 9 (項數 10), 餘項界 = 2.756e-7
cos(1)     目標 1e-10: 需要 n = 13 (項數 14), 餘項界 = 1.147e-11
cos(1)     目標 1e-15: 需要 n = 17 (項數 18), 餘項界 = 1.562e-16

ln(2)      目標 1e-06: Lagrange 界在 40 階內都不夠（改用交錯級數界）
ln(2)      目標 1e-10: Lagrange 界在 40 階內都不夠（改用交錯級數界）
ln(2)      目標 1e-15: Lagrange 界在 40 階內都不夠（改用交錯級數界）

arctan(1)  目標 1e-06: Lagrange 界在 40 階內都不夠（改用交錯級數界）
arctan(1)  目標 1e-10: Lagrange 界在 40 階內都不夠（改用交錯級數界）
arctan(1)  目標 1e-15: Lagrange 界在 40 階內都不夠（改用交錯級數界）

ln(2) = 1 - 1/2 + 1/3 - 1/4 + ...（交錯級數）
  目標 1e-06: 交錯級數界要求 n >= 1/eps ≈ 1000001 項
  目標 1e-10: 交錯級數界要求 n >= 1/eps ≈ 10000000001 項
  目標 1e-15: 交錯級數界要求 n >= 1/eps ≈ 1000000000000000 項

改寫 ln 2 = 2 arctanh(1/3) = 2 * sum_{k>=0} (1/3)^(2k+1)/(2k+1)
  目標 1e-06: 只要 6 項，和 = 0.693147073759785   誤差 = 1.068e-07
  目標 1e-10: 只要 9 項，和 = 0.693147180459244   誤差 = 1.007e-10
  目標 1e-15: 只要 15 項，和 = 0.693147180559945   誤差 = 1.110e-16

import sympy as sp

x = sp.symbols('x', real=True)
n = sp.symbols('n', integer=True, nonnegative=True)
k = sp.symbols('k', integer=True, positive=True)

# 比較 1/(1-x) 與 1/(1+x)：同樣的 series 結構，收斂行為卻不同
f1 = 1/(1-x)
f2 = 1/(1+x)
print("1/(1-x) 的 series =", sp.series(f1, x, 0, 6))
print("1/(1+x) 的 series =", sp.series(f2, x, 0, 6))
print("係數逐項比較（奇數項差一個正負號）:")
c1 = sp.series(f1, x, 0, 6).removeO().as_coefficients_dict()
c2 = sp.series(f2, x, 0, 6).removeO().as_coefficients_dict()
for key in sorted(c1, key=lambda t: sp.degree(t, x)):
    print(f"  {key}: a_k = {c1[key]}, b_k = {c2[key]}")

# 在 x = 0.9 與 x = -0.9 處分別計算部分和
for xv in [0.9, -0.9]:
    s1 = sum(sp.Rational(xv)**k for k in range(300))
    s2 = sum((-sp.Rational(xv))**k for k in range(300))
    print(f"\nx = {xv}:")
    print(f"  1/(1-x) 的 300 項部分和 = {float(s1):.6f}   真值 = {float(1/(1-xv)):.6f}")
    print(f"  1/(1+x) 的 300 項部分和 = {float(s2):.6f}   真值 = {float(1/(1+xv)):.6f}")

# 收斂半徑的符號判據：|a_k|^(1/k) 的上極限
print()
print("1/(1-x): a_k = 1 => limsup |a_k|^(1/k) =", sp.limit(1**(1/k), k, sp.oo))
print("  故 R = 1 / 1 = 1")
print("1/(1+x): a_k = (-1)^k => limsup |a_k|^(1/k) =", sp.limit(abs((-1)**k)**(1/k), k, sp.oo))
print("  故 R = 1 / 1 = 1")

# SymPy 的符號求和直接給出帶收斂條件的封閉形式（Piecewise）
print()
print("sum x^n (n=0..oo) =", sp.summation(x**n, (n, 0, sp.oo)))
print("sum (-x)^n (n=0..oo) =", sp.summation((-x)**n, (n, 0, sp.oo)))

# 收斂半徑之外會怎樣？
print()
for xv in [0.5, 0.9, 0.99, 1.0, 1.5, -1.5]:
    s = sum(sp.Rational(xv)**k for k in range(200))
    exact = float(1/(1-xv)) if xv != 1 else float('inf')
    print(f"x = {xv:>5}: 200 項部分和 = {float(s):>14.6f}   1/(1-x) = {exact:>12.6f}")

1/(1-x) 的 series = 1 + x + x**2 + x**3 + x**4 + x**5 + O(x**6)
1/(1+x) 的 series = 1 - x + x**2 - x**3 + x**4 - x**5 + O(x**6)
係數逐項比較（奇數項差一個正負號）:
  1: a_k = 1, b_k = 1
  x: a_k = 1, b_k = -1
  x**2: a_k = 1, b_k = 1
  x**3: a_k = 1, b_k = -1
  x**4: a_k = 1, b_k = 1
  x**5: a_k = 1, b_k = -1

x = 0.9:
  1/(1-x) 的 300 項部分和 = 10.000000   真值 = 10.000000
  1/(1+x) 的 300 項部分和 = 0.526316   真值 = 0.526316

x = -0.9:
  1/(1-x) 的 300 項部分和 = 0.526316   真值 = 0.526316
  1/(1+x) 的 300 項部分和 = 10.000000   真值 = 10.000000

1/(1-x): a_k = 1 => limsup |a_k|^(1/k) = 1
  故 R = 1 / 1 = 1
1/(1+x): a_k = (-1)^k => limsup |a_k|^(1/k) = 1
  故 R = 1 / 1 = 1

sum x^n (n=0..oo) = Piecewise((1/(1 - x), Abs(x) < 1), (Sum(x**n, (n, 0, oo)), True))
sum (-x)^n (n=0..oo) = Piecewise((1/(x + 1), Abs(x) < 1), (Sum((-1)**n*x**n, (n, 0, oo)), True))

x =   0.5: 200 項部分和 =       2.000000   1/(1-x) =     2.000000
x =   0.9: 200 項部分和 =      10.000000   1/(1-x) =    10.000000
x =  0.99: 200 項部分和 =      86.602033   1/(1-x) =   100.000000
x =   1.0: 200 項部分和 =     200.000000   1/(1-x) =          inf
x =   1.5: 200 項部分和 = 330583982157641620045218217747021824.000000   1/(1-x) =    -2.000000
x =  -1.5: 200 項部分和 = -66116796431528318475020421436538880.000000   1/(1-x) =     0.400000

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 5))
xs = np.linspace(-2.5, 2.5, 800)

for N, color in [(5, "#c0392b"), (10, "#2980b9"), (30, "#27ae60"), (200, "#8e44ad")]:
    ys = np.array([sum(x**k for k in range(N+1)) for x in xs])
    ax.plot(xs, ys, color=color, label=f"$\\sum_{{k=0}}^{{{N}}} x^k$")

exact = np.where(np.abs(xs) < 1, 1/(1-xs), np.nan)
ax.plot(xs, exact, "k--", linewidth=2, label="$1/(1-x)$")
for v in (1, -1):
    ax.axvline(v, color="gray", linestyle=":", linewidth=1.5)
ax.set_ylim(-5, 5); ax.set_xlabel("x"); ax.set_ylabel("y")
ax.legend(loc="upper left"); ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.savefig("/tmp/geometric_series.png", dpi=120)

print("收斂半徑 R = 1；|x| < 1 時部分和 -> 1/(1-x)，|x| > 1 時發散")
print("saved /tmp/geometric_series.png")

收斂半徑 R = 1；|x| < 1 時部分和 -> 1/(1-x)，|x| > 1 時發散
saved /tmp/geometric_series.png

import sympy as sp
import numpy as np

x = sp.symbols('x', real=True)

# 數值微分：中心差分 = 二階 Taylor 展開
h = sp.symbols('h', positive=True)
f = sp.sin(x)
center = sp.series(f.subs(x, x + h) - f.subs(x, x - h), h, 0, 5)
print("sin 的中心差分展開 =", center)
print("=> (sin(x+h)-sin(x-h))/(2h) = cos(x) + O(h^2)")

# 符號驗證：中心差分與真導數的誤差是 O(h^2)
cd = sp.series((f.subs(x, x+h) - f.subs(x, x-h))/(2*h), h, 0, 3)
print("中心差分 =", cd)

# 插值：等距節點的 Lagrange 插值多項式 = 差點表
print()
pts = [(0, sp.Rational(1,1)), (1, sp.Rational(3,2)), (2, sp.Rational(14,3))]
xs = sp.symbols('xs')
lag = sum(y * sp.prod([(xs - xj)/(xi - xj) for (xj, _) in pts if xj != xi])
           for (xi, y) in pts)
print("Lagrange 插值多項式 =", sp.expand(lag))
print("代入 x=1 得", sp.simplify(lag.subs(xs, 1)), " 應為 3/2")

# 插值多項式 = 該點的 Taylor 多項式（余項解釋誤差）
print()
f2 = sp.exp(xs)
T2 = sp.series(f2, xs, 0, 3).removeO()
print("e^x 的 T_2 =", T2, "  在 x=1 的值 =", sp.N(T2.subs(xs, 1)))

# 插值與微分的關係：差商 -> 導數
print()
for hv in [0.1, 0.01, 0.001]:
    num = float((sp.sin(1+hv) - sp.sin(1-hv))/(2*hv))
    print(f"h = {hv}: 中心差分 = {num:.12f}   真值 cos(1) = {float(sp.cos(1)):.12f}   誤差 = {abs(num-float(sp.cos(1))):.3e}")

sin 的中心差分展開 = 2*h*cos(x) - h**3*cos(x)/3 + O(h**5)
=> (sin(x+h)-sin(x-h))/(2h) = cos(x) + O(h^2)
中心差分 = cos(x) - h**2*cos(x)/6 + O(h**3)

Lagrange 插值多項式 = 4*xs**2/3 - 5*xs/6 + 1
代入 x=1 得 3/2  應為 3/2

e^x 的 T_2 = xs**2/2 + xs + 1   在 x=1 的值 = 2.50000000000000

h = 0.1: 中心差分 = 0.539402252170   真值 cos(1) = 0.540302305868   誤差 = 9.001e-04
h = 0.01: 中心差分 = 0.540293300875   真值 cos(1) = 0.540302305868   誤差 = 9.005e-06
h = 0.001: 中心差分 = 0.540302215818   真值 cos(1) = 0.540302305868   誤差 = 9.005e-08