Pixel Pedals of Tomakomai

北海道苫小牧市出身の初老の日常

問題3.19の解答(マクローリン展開)

概説微分積分の問題3.19の解答。

exp(0.5)、cos(1)、log(1.1)を小数第二位まで求める問題。ラグランジュの剰余項Rn(x)の絶対値を評価する必要があるのだけど、評価方法を経験則的にやってるので一般化できてない。

import math

def trunc(x):
    if x > 0: return math.floor(x)
    return math.ceil(x)

def solve(x, taylor_term, eval_R):
    n = 0
    while True:
        # calc f(x)
        result = 0.
        for i in range(0, n + 1):
            result += taylor_term(x, i)

        # calc |Rn(x)|
        r = eval_R(x, n + 1)

        # Check the range of f(x)
        min_fx = trunc((result + r) * 100) / 100
        max_fx = trunc((result - r) * 100) / 100

        print "f(x) = %.4f, |Rn(x)| = %.4f" % (result, r)
        if min_fx == max_fx: return min_fx, n

        # run next loop
        n += 1

def factor(n):
    factor_n = 1
    for i in range(2, n + 1): factor_n *= i

    return factor_n

def taylor_term1(x, n):
    return float(x ** n) / factor(n)

def R1(x, n):
    # since c = 1/2 < 2
    return 2. * (x ** n) / factor(n)

def taylor_term2(x, n):
    if n % 2 == 1: return 0
    return (float(x) ** n) * ((-1) ** (n / 2)) / factor(n)

def R2(x, n):
    # since sin(or cos) <= 1
    return (float(x) ** n) / factor(n)

def taylor_term3(x, n):
    if n == 0: return 0
    return float(x ** n) * ((-1) ** (n - 1)) / n

def R3(x, n):
    # since 1 + c > 1
    return float(x ** n) / float(n)

print "f(x) = %.2f, n = %d" % solve(1. / 2, taylor_term1, R1)
print "f(x) = %.2f, n = %d" % solve(1., taylor_term2, R2)
print "f(x) = %.2f, n = %d" % solve(.1, taylor_term3, R3)
f(x) = 1.0000, |Rn(x)| = 1.0000
f(x) = 1.5000, |Rn(x)| = 0.2500
f(x) = 1.6250, |Rn(x)| = 0.0417
f(x) = 1.6458, |Rn(x)| = 0.0052
f(x) = 1.6484, |Rn(x)| = 0.0005
f(x) = 1.64, n = 4
f(x) = 1.0000, |Rn(x)| = 1.0000
f(x) = 1.0000, |Rn(x)| = 0.5000
f(x) = 0.5000, |Rn(x)| = 0.1667
f(x) = 0.5000, |Rn(x)| = 0.0417
f(x) = 0.5417, |Rn(x)| = 0.0083
f(x) = 0.5417, |Rn(x)| = 0.0014
f(x) = 0.54, n = 5
f(x) = 0.0000, |Rn(x)| = 0.1000
f(x) = 0.1000, |Rn(x)| = 0.0050
f(x) = 0.0950, |Rn(x)| = 0.0003
f(x) = 0.09, n = 2