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
74
75
76
77
78
79
80
|
#!/usr/bin/python3
""" Plot phi0, and (later) anaylze approximations
"""
import numpy as np
import matplotlib.pyplot as plt
import math
import itertools
import sys
def impl_reference(x):
if (x< 9.08e-5 ): r = 10
else:
z = math.exp(x)
r = math.log( (z+1) / (z-1) )
return(r)
####
# Table T1
t1_tbl = [] # A list of tuples (x, r)
t1_adj = 0.50 # Offset ratio for step of 1, estimate
def t1_loop(upper, lower, step):
for i in range((int)(upper * step), (int)(lower * step), -1):
x = i / step
z = math.exp(x + (t1_adj / step))
r = math.log((z+1)/(z-1))
t1_tbl.append( (x, r) )
print(x, r)
# 9:5 step 1
t1_loop(9, 4.999, 2)
# 5:1 step 1/16
t1_loop(4.999, 0.999, 16)
# 1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
for i in range(1, 28):
e = (2 ** (i/2))
x = 1.0 / e
z = math.exp(x + (0.19 / e))
r = math.log((z+1)/(z-1))
t1_tbl.append( (x, r) )
def impl_t1(x):
if (x > 10): return(0)
else:
for t in t1_tbl:
if (x > t[0]): return(t[1])
return(10)
##########################
# Plot, scanning from 10 to 0
x_vals = np.logspace(1, -4, 2000, endpoint=False)
ref_vals = []
t1_vals = []
for x in x_vals:
ref_vals.append(impl_reference(x))
t1_vals.append(impl_t1(x))
# Sum errors
errsum = errsum2 = 0
for i in range(len(ref_vals)):
error = t1_vals[i] - ref_vals[i]
errsum += error
errsum2 += error * error
print("Net error {}".format(errsum))
print("avg error {}".format(errsum/len(ref_vals)))
print("rms error {}".format(math.sqrt(errsum2/len(ref_vals))))
plt.xscale('log')
plt.plot(x_vals, ref_vals, 'g', x_vals, t1_vals, 'r')
plt.show()
|
