1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
"""
Interval arithmetic demo: estimating error of numerical Taylor series.
This module can be run interactively with
python taylor.py
"""
from __future__ import print_function
from mpmath import *
def taylor(x, n):
print("-"*75)
t = x = mpi(x)
s = 1
print("adding 1")
print(s, "\n")
s += t
print("adding x")
print(s, "\n")
for k in range(2, n+1):
t = (t * x) / k
s += t
print("adding x^%i / %i! ~= %s" % (k, k, t.mid))
print(s, "\n")
print("-"*75)
return s
# Note: this should really be computed using interval arithmetic too!
def remainder(x, n):
xi = max(0, x)
r = exp(xi) / factorial(n+1)
r = r * x**(n+1)
return abs(r)
def exponential(x, n):
"""
Compute exp(x) using n terms of the Taylor series for exp using
intervals, and print detailed error analysis.
"""
t = taylor(x, n)
r = remainder(x, n)
expx = exp(x)
print("Correct value of exp(x): ", expx)
print()
print("Computed interval: ")
print(t)
print()
print("Computed value (midpoint): ", t.mid)
print()
print("Estimated rounding error: ", t.delta)
print("Estimated truncation error: ", r)
print("Estimated total error: ", t.delta + r)
print("Actual error ", abs(expx - t.mid))
print()
u = t + mpi(-r, r)
print("Interval with est. truncation error added:")
print(u)
print()
print("Correct value contained in computed interval:", t.a <= expx <= t.b)
print("When accounting for truncation error:", u.a <= expx <= u.b)
if __name__ == "__main__":
print("Interval arithmetic demo")
print()
print("This script sums the Taylor series for exp(x) using interval arithmetic,")
print("and then compares the numerical errors due to rounding and truncation.")
print()
x = mpf(raw_input("Enter the value of x (e.g. 3.5): "))
n = int(raw_input("Enter the number of terms n (e.g. 10): "))
print()
exponential(x, n)
|
