Chebyshev coefficients for $e^{-x}$ for $n=1$ from Cody, Meinardus, and Varga (table III)
n = 1
---------------------------------
i p_i q_i
---------------------------------
0 1.0669 ( 00) 1.0000 ( 00)
1 -1.1535 (-01) 1.7275 ( 00)
In [4]:
from sympy import init_session
init_session()
t = symbols('t', real=True)
%matplotlib inline
(From Maria_Pusa2.pdf)
Remez algorithm
Assume $\left\{t_i \right\}_{i=1}^{2k+2} \subset [−1, 1]$ and find real polynomials $p_k$ and $q_k$ and a parameter $\epsilon > 0$ such that $$ \begin{align} e^{\phi(t_i)} − \frac{p(\phi(t_i))}{q(\phi(t_i))} +(−1)^i\epsilon=0,&& i=1,...,2k+2,\\ \end{align} $$ where $q_{k+1} = 1$.
Assume $r_{k,k} \in \pi_{k,k}$ and $\epsilon > 0$ and find the $2k + 2$ local extreme points of the function $$ E(t) = e^{\phi(t)} - r_{k,k}(\phi(t)) $$ in the interval $[−1, 1]$.
In [2]:
[chebyshevt_root(5, i) for i in range(5)]
Out[2]:
In [3]:
[chebyshevt_root(5, i).evalf() for i in range(5)]
Out[3]:
In [4]:
chebyshevt(5, x)
Out[4]:
In [5]:
epsilon = symbols("epsilon")
p0, p1, q0, q1 = symbols("p0, p1, q0, q1")
i = symbols("i")
In [6]:
expr = exp(-t) - (p0 + p1*t)/(q0 + q1*t) + (-1)**i*epsilon
expr = expr*(q0 + q1*t)
expr = simplify(expr)
expr
Out[6]:
In [7]:
system = Tuple(*[expr.subs({i: j, t: chebyshevt_root(5, 4-j).evalf()}) for j in range(1,5)])
In [8]:
system = system.subs(q0, 1)
system
Out[8]:
In [9]:
sols = solve(system, [p0, p1, q1, epsilon], dict=True)
sols
Out[9]:
In [10]:
E = exp(-t) - (p0 + p1*t)/(1 + q1*t)
In [11]:
plot(E.subs(sols[0]), (t, -1, 1))
Out[11]:
In [12]:
extreme_x2 = [-1, nsolve(diff(E.subs(sols[0]), t), -1), nsolve(diff(E.subs(sols[0]), t), 0), 1]
extreme_x2
Out[12]:
In [13]:
system2 = Tuple(*[expr.subs({i: j, t: extreme_x2[j]}) for j in range(4)]).subs(q0, 1)
system2
Out[13]:
In [14]:
sols2 = solve(system2, [p0, p1, q1, epsilon], dict=True)
sols2
Out[14]:
In [15]:
plot(E.subs(sols2[1]), (t, -1, 1))
Out[15]:
In [16]:
extreme_x3 = [-1, nsolve(diff(E.subs(sols2[1]), t), -1), nsolve(diff(E.subs(sols2[1]), t), 1), 1]
extreme_x3
Out[16]:
In [17]:
system3 = Tuple(*[expr.subs({i: j, t: extreme_x3[j]}) for j in range(4)]).subs(q0, 1)
system3
Out[17]:
In [18]:
sols3 = solve(system3, [p0, p1, q1, epsilon], dict=True)
sols3
Out[18]:
In [19]:
plot(E.subs(sols3[1]), (t, -1, 1))
Out[19]:
In [20]:
[E.subs(sols3[1]).subs(t, i).evalf() for i in extreme_x3]
Out[20]:
In [21]:
extreme_x4 = [-1, nsolve(diff(E.subs(sols3[1]), t), -1), nsolve(diff(E.subs(sols3[1]), t), 1), 1]
extreme_x4
Out[21]:
In [22]:
system4 = Tuple(*[expr.subs({i: j, t: extreme_x4[j]}) for j in range(4)]).subs(q0, 1)
system4
Out[22]:
In [23]:
sols4 = solve(system4, [p0, p1, q1, epsilon], dict=True)
sols4
Out[23]:
In [24]:
plot(E.subs(sols4[1]), (t, -1, 1))
Out[24]:
In [25]:
[E.subs(sols3[1]).subs(t, i).evalf() for i in extreme_x3]
Out[25]:
In [ ]:
In [ ]:
Another useful resource: http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/backgrounders/remez.html
The method usually adopted to solve these equations is an iterative one: we guess the value of E, solve the equations to obtain a new value for E (as well as the polynomial coefficients), then use the new value of E as the next guess. The method is repeated until E converges on a stable value.
These complications extend the running time required for the development of rational approximations quite considerably. It is often desirable to obtain a rational rather than polynomial approximation none the less: rational approximations will often match more difficult to approximate functions, to greater accuracy, and with greater efficiency, than their polynomial alternatives. For example, if we takes our previous example of an approximation to ex, we obtained 5x10-4 accuracy with an order 4 polynomial. If we move two of the unknowns into the denominator to give a pair of order 2 polynomials, and re-minimise, then the peak relative error drops to 8.7x10-5. That's a 5 fold increase in accuracy, for the same number of terms overall.
In [26]:
expr = exp(-(-t - 1)/(2*t - 2)) - (p0 + p1*t)/(q0 + q1*t) + (-1)**i*epsilon
expr = expr*(q0 + q1*t)
expr = simplify(expr)
expr
Out[26]:
In [27]:
system = Tuple(*[expr.subs({i: j, t: chebyshevt_root(5, 4-j).evalf()}) for j in range(1,5)])
In [28]:
system = system.subs(q0, 1)
system
Out[28]:
In [29]:
sols = solve(system, [p0, p1, q1, epsilon], dict=True)
sols
Out[29]:
In [30]:
E = exp(-(-t - 1)/(2*t - 2)) - (p0 + p1*t)/(1 + q1*t)
In [31]:
plot(E.subs(sols[0]), (t, -1, 1))
Out[31]:
In [ ]:
In [ ]:
In [32]:
plot(E.subs(sols[0]).diff(t), (t, -1, 1))
Out[32]:
In [33]:
print(E.subs(sols[0]).diff(t))
In [34]:
nsolve(diff(E.subs(sols[0]), t), (.5, 0.9), solver='bisect')
Out[34]:
In [35]:
nsolve(diff(E.subs(sols[0]), t), 0.9)
Out[35]:
In [36]:
E.subs(sols[0]).diff(t).subs(t, 1)
Out[36]:
In [37]:
E.subs(sols[0]).diff(t).subs(t, 0.99996577349047597)
Out[37]:
This is the wrong answer. See https://github.com/sympy/sympy/issues/11768
In [38]:
E.subs(sols[0]).diff(t).subs(t, 0.7)
Out[38]:
In [ ]:
In [ ]:
We have to use the weird wrong solution here, because expr(1) = infinity ???
In [39]:
extreme_x2 = [-1, nsolve(diff(E.subs(sols[0]), t), 0), nsolve(diff(E.subs(sols[0]), t), (0.5, .9), solver='bisect'), nsolve(diff(E.subs(sols[0]), t), 0.9)]
extreme_x2
Out[39]:
In [40]:
system2 = Tuple(*[expr.subs({i: j, t: extreme_x2[j]}) for j in range(4)]).subs(q0, 1)
system2
Out[40]:
In [41]:
sols2 = solve(system2, [p0, p1, q1, epsilon], dict=True)
sols2
Out[41]:
In [42]:
plot(E.subs(sols2[1]), (t, -1, 1))
Out[42]:
In [43]:
extreme_x3 = [-1, nsolve(diff(E.subs(sols2[1]), t), -1), nsolve(diff(E.subs(sols2[1]), t), (0.5, 0.9), solver='bisect'), nsolve(diff(E.subs(sols2[1]), t), .9)]
extreme_x3
Out[43]:
In [44]:
plot(exp(-(-t - 1)/(2*t - 2)), (t, -1, 1))
Out[44]:
In [ ]:
In [45]:
system3 = Tuple(*[expr.subs({i: j, t: extreme_x3[j]}) for j in range(4)]).subs(q0, 1)
system3
Out[45]:
In [46]:
sols3 = solve(system3, [p0, p1, q1, epsilon], dict=True)
sols3
Out[46]:
In [47]:
plot(E.subs(sols3[1]), (t, -1, 1))
Out[47]:
In [48]:
[E.subs(sols3[1]).subs(t, i).evalf() for i in extreme_x3]
Out[48]:
In [49]:
r = (p0 + p1*t)/(1 + q1*t)
In [50]:
r.subs(sols3[1])
Out[50]:
In [51]:
n, d = simplify(r.subs(sols3[1]).subs(t, (2*t - 1)/(2*t + 1))).as_numer_denom()
rat_func = (n/Poly(d).TC())/(d/Poly(d).TC())
correct_rat_func = (1.0669 + -1.1535e-1*t)/(1 + 1.7275*t)
rat_func, correct_rat_func
Out[51]:
In [52]:
plot((rat_func, (t, 0, 100)), (exp(-t), (t, 0, 100)), ylim=(-1, 10))
Out[52]:
In [53]:
plot((correct_rat_func - exp(-t), (t, 0, 100)))
Out[53]:
In [54]:
plot((rat_func - exp(-t), (t, 0, 100)))
Out[54]:
That's it!!!
For reference, these are the Möbius transforms for $[-1, 1] \leftrightarrow [0, \infty]$
In [55]:
simplify(1/(1 - x) - S(1)/2)
Out[55]:
In [56]:
solve(1/(1 - x) - S(1)/2 - y, x)
Out[56]:
Now let's try $n = 2$
Chebyshev coefficients for $e^{-x}$ for $n=1$ from Cody, Meinardus, and Varga (table III)
n = 2
---------------------------------
i p_i q_i
---------------------------------
0 9.92641 (-01) 1.00000 ( 00)
1 -1.88332 (-01) 6.69295 (-01)
2 4.21096 (-03) 5.72258 (-01)
In [57]:
[chebyshevt_root(6, i) for i in range(6)]
Out[57]:
In [58]:
epsilon = symbols("epsilon")
p0, p1, p2, q1, q2 = symbols("p0, p1, p2, q1, q2")
i = symbols("i")
In [59]:
expr = exp(-(-t - 1)/(2*t - 2)) - (p0 + p1*t + p2*t**2)/(1 + q1*t + q2*t**2) + (-1)**i*epsilon
expr = expr*(1 + q1*t + q2*t**2)
expr = simplify(expr)
expr
Out[59]:
In [60]:
system = Tuple(*[expr.subs({i: j, t: chebyshevt_root(7, 6-j)}) for j in range(1,7)])
In [61]:
system = system.subs(q0, 1)
system
Out[61]:
In [62]:
sols = [dict(zip([p0, p1, p2, q1, q2, epsilon], nsolve(system, [p0, p1, p2, q1, q2, epsilon], [1, 1, 1, 1, 1, 0])))]
sols
Out[62]:
In [63]:
E = exp(-(-t - 1)/(2*t - 2)) - (p0 + p1*t + p2*t**2)/(1 + q1*t + q2*t**2)
In [64]:
plot(E.subs(sols[0]), (t, -1, 1))
Out[64]:
In [65]:
E.subs(sols[0]).subs(t, 1)
Out[65]:
In [66]:
extreme_x2 = [-1, nsolve(diff(E.subs(sols[0]), t), (-1, 0), solver='bisect'), nsolve(diff(E.subs(sols[0]), t), (0, 0.5), solver='bisect'), nsolve(diff(E.subs(sols[0]), t), (0.5, .7), solver='bisect'), nsolve(diff(E.subs(sols[0]), t), 0.9), 0.999]
extreme_x2
Out[66]:
In [67]:
system2 = Tuple(*[expr.subs({i: j, t: extreme_x2[j]}) for j in range(6)]).subs(q0, 1)
system2
Out[67]:
In [68]:
sols2 = [dict(zip([p0, p1, p2, q1, q2, epsilon], nsolve(system2, [p0, p1, p2, q1, q2, epsilon], [1, 1, 1, 1, 1, 0])))]
sols2
Out[68]:
In [69]:
plot(E.subs(sols2[0]), (t, -1, 1))
Out[69]:
In [70]:
extreme_x3 = [-1, nsolve(diff(E.subs(sols2[0]), t), (-1, 0), solver='bisect'), nsolve(diff(E.subs(sols2[0]), t), (0, 0.5), solver='bisect'), nsolve(diff(E.subs(sols2[0]), t), (0.5, .8), solver='bisect'), nsolve(diff(E.subs(sols2[0]), t), (0.8, 0.99), solver='bisect'), 0.999]
extreme_x3
Out[70]:
In [71]:
system3 = Tuple(*[expr.subs({i: j, t: extreme_x3[j]}) for j in range(6)]).subs(q0, 1)
system3
Out[71]:
In [72]:
sols3 = [dict(zip([p0, p1, p2, q1, q2, epsilon], nsolve(system3, [p0, p1, p2, q1, q2, epsilon], [1, 1, 1, 1, 1, 0])))]
sols3
Out[72]:
In [73]:
plot(E.subs(sols3[0]), (t, -1, 1))
Out[73]:
In [74]:
[E.subs(sols3[0]).subs(t, i).evalf() for i in extreme_x3]
Out[74]:
In [75]:
r = (p0 + p1*t + p2*t**2)/(1 + q1*t + q2*t**2)
In [76]:
r.subs(sols3[0])
Out[76]:
In [77]:
n, d = together(r.subs(sols3[0]).subs(t, (2*t - 1)/(2*t + 1))).as_numer_denom() # simplify/cancel here will add degree to the numerator and denominator
rat_func = (Poly(n)/Poly(d).TC())/(Poly(d)/Poly(d).TC())
correct_rat_func = (9.92641e-1 + -1.88332e-1*t + 4.21096e-3*t**2)/(1 + 6.69295e-1*t + 5.72258e-1*t**2)
rat_func, correct_rat_func
Out[77]:
In [78]:
plot((rat_func, (t, 0, 100)), (exp(-t), (t, 0, 100)), ylim=(-1, 10))
Out[78]:
In [79]:
plot((correct_rat_func - exp(-t), (t, 0, 100)))
Out[79]:
In [80]:
plot((rat_func - exp(-t), (t, 0, 100)))
Out[80]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [81]:
N = symbols("N")
x = symbols("x", real=True)
phi = N*(0.1309 - 0.1149*x**2 + I*0.2500*x)
In [82]:
phi.subs(N, 4).as_real_imag()
Out[82]:
In [83]:
from sympy.plotting import *
plot_parametric((*phi.subs(N, 4).as_real_imag(), (x, -10, 10)))
Out[83]:
In [84]:
correct_rat_func = (1.0669 + -1.1535e-1*t)/(1 + 1.7275*t)
correct_rat_func2 = (9.92641e-1 + -1.88332e-1*t + 4.21096e-3*t**2)/(1 + 6.69295e-1*t + 5.72258e-1*t**2)
In [85]:
apart(correct_rat_func)
Out[85]:
In [86]:
apart(correct_rat_func2, full=True).doit()
Out[86]:
$p_i$
0 9.9999999999998168e-01
1 -2.7495604296300043e-01
2 3.4346984175671475e-02
3 -2.5674439819028618e-03
4 1.2734070715233181e-04
5 -4.3932808492511236e-06
6 1.0753202054485227e-07
7 -1.8710255961089453e-09
8 2.2849515765300155e-11
9 -1.9038640942835345e-13
10 1.0310151365350495e-15
11 -3.3490280333667533e-18
12 5.6743825539523501e-21
13 -3.7794523874503295e-24
14 4.1609826642376613e-28$q_i$
1 1.0000000000000000e+00
2 7.2504395703488666e-01
3 2.5939094125018012e-01
4 6.0968185127283595e-02
5 1.0574049161691156e-02
6 1.4405527154587316e-03
7 1.6019211440612190e-04
8 1.4908234724800242e-05
9 1.1820291576355772e-06
10 8.0163997982357503e-08
11 4.8361800878648281e-09
12 2.2472030400428529e-10
13 1.2922802705792779e-11
14 1.8921838854022449e-13
15 2.2710680218891295e-14
In [87]:
n14pis = """
0 9.9999999999998168e-01
1 -2.7495604296300043e-01
2 3.4346984175671475e-02
3 -2.5674439819028618e-03
4 1.2734070715233181e-04
5 -4.3932808492511236e-06
6 1.0753202054485227e-07
7 -1.8710255961089453e-09
8 2.2849515765300155e-11
9 -1.9038640942835345e-13
10 1.0310151365350495e-15
11 -3.3490280333667533e-18
12 5.6743825539523501e-21
13 -3.7794523874503295e-24
14 4.1609826642376613e-28
"""
n14qis = """
1 1.0000000000000000e+00
2 7.2504395703488666e-01
3 2.5939094125018012e-01
4 6.0968185127283595e-02
5 1.0574049161691156e-02
6 1.4405527154587316e-03
7 1.6019211440612190e-04
8 1.4908234724800242e-05
9 1.1820291576355772e-06
10 8.0163997982357503e-08
11 4.8361800878648281e-09
12 2.2472030400428529e-10
13 1.2922802705792779e-11
14 1.8921838854022449e-13
15 2.2710680218891295e-14
"""
n14pi = [Float(i.split()[-1]) for i in n14pis.strip().split('\n')]
n14qi = [Float(i.split()[-1]) for i in n14qis.strip().split('\n')]
n14qi[0] = 1 # Exactly
n14pi, n14qi
Out[87]:
In [88]:
correct_rat_func14 = Poly(reversed(n14pi), t)/Poly(reversed(n14qi), t)
correct_rat_func14
Out[88]:
In [89]:
plot(exp(-t) - correct_rat_func14, (t, 0, 10))
Out[89]:
In [90]:
plot(exp(-t) - correct_rat_func14, (t, 1e-4, 1e2), xscale='log', adaptive=False, nb_of_points=10000)
Out[90]:
In [91]:
# Use x, which has no assumptions, to work around a bug
#apart_correct_rat_func14 = apart(correct_rat_func14.subs(t, x), x, domain=QQ, full=True).doit().evalf().expand().subs(x, t)
#apart_correct_rat_func14
In [92]:
#plot(exp(-t) - re(apart_correct_rat_func14), (t, 1e-4, 1e2), xscale='log', adaptive=False, nb_of_points=10000)
In [93]:
#re(apart_correct_rat_func14)
In [94]:
#cancel(apart_correct_rat_func14)
In [95]:
correct_rat_func14
Out[95]:
In [96]:
srepr(correct_rat_func14)
Out[96]:
In [1]:
real_part14 = """
"-8.897 773 186468 8888199 X 100
“3.703 275 049423 448 0603 X 100
-0.208 758 638 250 130 125 l X 100
+3.993 3697105785685194 X 100
+5.089 345 060580 6245066 X 100
+5.623 142 572745 977 1248 X 100
+2.269 783 8292311127097 X 100
“7.154288 0635890672853 X 10""5
+9.439 025 310736168 877 9 X 10“3
“3.763 600 387 822 696 871 7 X 10“1
-2.349 823 2091082701191 X 1001
+4.693 327 448 8831293047 X 101
~2.787516194014564646 8 X 101
+4.807 112098 832508 8907 X 100
+1.832 174 378 254041 275 1 X 10“14
""".strip().replace(' ', '').replace('"', '-').replace('“', '-').replace('X10', 'e').replace('--', '-').replace('l', '1').replace('~', '-')
print(real_part14)
In [5]:
theta1r, theta2r, theta3r, theta4r, theta5r, theta6r, theta7r, alpha1r, alpha2r, alpha3r, alpha4r, alpha5r, alpha6r, alpha7r, alpha0r = map(Float, real_part14.strip().split())
theta1r, theta2r, theta3r, theta4r, theta5r, theta6r, theta7r, alpha1r, alpha2r, alpha3r, alpha4r, alpha5r, alpha6r, alpha7r, alpha0r
Out[5]:
In [6]:
imaginary_part14 = """
+ 1.663 098 261 990208 5304 X 101
+ 1.365 637 187 148 3268171 X 101
+1.0991260561901260913 X 101
+6.004 831 642235 0373178 X 100
+3.588 8240290270065102 X 100
+1.194069046 343 966 9766 X 100
+8.461 737 973 0402214019 X 100
+1.436104334954130011 1 X 10““4
--1.718479195 848 301 751 1 X 10-2
+3.3518347029450104214 X 10""1
-5.808 359 129714207 4004 X 100
+4.564 364976 882 776 0791 X 101
- 1.021 473 399 905 645 143 4 X 102
- 1.320 979 383742 872 3881 X 100
+ 0.000 000 000 000 000 000 0 X 100
""".strip().replace(' ', '').replace('"', '-').replace('“', '-').replace('X10', 'e').replace('--', '-').replace('l', '1').replace('~', '-')
print(imaginary_part14)
In [7]:
theta1j, theta2j, theta3j, theta4j, theta5j, theta6j, theta7j, alpha1j, alpha2j, alpha3j, alpha4j, alpha5j, alpha6j, alpha7j, alpha0j = [Float(i)*I for i in imaginary_part14.strip().split()]
theta1j, theta2j, theta3j, theta4j, theta5j, theta6j, theta7j, alpha1j, alpha2j, alpha3j, alpha4j, alpha5j, alpha6j, alpha7j, alpha0j
Out[7]:
In [8]:
theta1, theta2, theta3, theta4, theta5, theta6, theta7, alpha1, alpha2, alpha3, alpha4, alpha5, alpha6, alpha7, alpha0 =\
(r + j for r, j in zip(
(theta1r, theta2r, theta3r, theta4r, theta5r, theta6r, theta7r, alpha1r, alpha2r, alpha3r, alpha4r, alpha5r, alpha6r, alpha7r, alpha0r),
(theta1j, theta2j, theta3j, theta4j, theta5j, theta6j, theta7j, alpha1j, alpha2j, alpha3j, alpha4j, alpha5j, alpha6j, alpha7j, alpha0j)))
In [9]:
thetas = theta1, theta2, theta3, theta4, theta5, theta6, theta7
alphas = alpha1, alpha2, alpha3, alpha4, alpha5, alpha6, alpha7
In [10]:
correct_partial_frac14 = alpha0 + Add(*(alpha/(t - theta) + alpha/(t - conjugate(theta)) for alpha, theta in zip(alphas, thetas)))
In [12]:
correct_re_partial_frac14 = alpha0 + 2*re(Add(*(alpha/(t - theta) for alpha, theta in zip(alphas, thetas))))
print(correct_re_partial_frac14)
In [300]:
correct_partial_frac14
Out[300]:
In [105]:
a, d = map(lambda i: Poly(i, t), cancel(correct_partial_frac14).as_numer_denom())
collect((a*expand_complex(1/d.TC())), t)/collect((d*expand_complex(1/d.TC())), t)
Out[105]:
In [106]:
correct_rat_func14
Out[106]:
There are a lot of errors here. I think there must be a lot of loss of precision. Also, I think one expression is representing $e^t$ and one $e^{-t}$, which would explain the wrong signs on the odd-powered $t$ terms. See https://github.com/sympy/sympy/issues/11795.
In [107]:
together(correct_partial_frac14)
Out[107]:
In [108]:
# var('_w _a')
# issue_answer = S("-1.0*RootSum(2.2710680218891296e-14*_w**14 + 1.8921838854022448e-13*_w**13 + 1.292280270579278e-11*_w**12 + 2.2472030400428528e-10*_w**11 + 4.8361800878648284e-9*_w**10 + 8.0163997982357503e-8*_w**9 + 1.1820291576355772e-6*_w**8 + 1.4908234724800242e-5*_w**7 + 0.00016019211440612189*_w**6 + 0.0014405527154587316*_w**5 + 0.010574049161691156*_w**4 + 0.060968185127283595*_w**3 + 0.25939094125018014*_w**2 + 0.72504395703488667*_w + 1.0, Lambda(_a, (3.7997502361189504e-11*_a**13 - 9.6301335874438056e-11*_a**12 + 2.2157750448607173e-8*_a**11 + 1.4101105165363766e-7*_a**10 + 6.2233711263991178e-6*_a**9 + 6.7342488958142767e-5*_a**8 + 0.0011644621826545941*_a**7 + 0.012030053971086767*_a**6 + 0.12615589337787736*_a**5 + 0.97182363139568362*_a**4 + 6.2658937818543281*_a**3 + 28.994072695533735*_a**2 + 90.256649069691846*_a + 138.65741216748256)/(-_a + x))) + 1.8321699852813983e-14")
# issue_answer
In [109]:
correct_rat_func14
Out[109]:
In [110]:
num = correct_rat_func14.as_numer_denom()[0]
In [111]:
Poly(num, domain=QQ)
Out[111]:
In [112]:
Poly(num, domain=QQ).as_expr().evalf(17)
Out[112]:
In [113]:
num.replace(lambda i: i.is_Float, Rational)
Out[113]:
In [114]:
Float("1.01").replace(lambda i: i.is_Float, lambda i: Rational(i))
Out[114]:
In [115]:
correct_rat_func14_rational = Poly(reversed(list(map(Rational, n14pi))), t)/Poly(reversed(list(map(Rational, n14qi))), t)
correct_rat_func14_rational
Out[115]:
In [116]:
#apart_correct_rat_func14_rational = apart(correct_rat_func14_rational.subs(t, x), x, full=True).doit().evalf(17).expand().subs(x, t)
#apart_correct_rat_func14_rational
In [117]:
#apart_correct_rat_func14_rational*2
In [118]:
#correct_partial_frac14
In [119]:
#apart(correct_rat_func14_rational.subs(t, x), x, full=True).doit()
In [120]:
correct_partial_frac14_rational = Rational(re(alpha0)) + Rational(im(alpha0))*I + Add(*((Rational(re(alpha)) + Rational(im(alpha))*I)/(t - Rational(re(theta)) - Rational(im(theta))*I) + (Rational(re(alpha)) + Rational(im(alpha))*I)/(t - Rational(re(theta)) + Rational(im(theta))*I) for alpha, theta in zip(alphas, thetas)))
correct_partial_frac14_rational
Out[120]:
In [121]:
a, d = map(lambda i: Poly(i, t), cancel(correct_partial_frac14_rational).as_numer_denom())
expr = collect((a*expand_complex(1/d.TC())), t)/collect((d*expand_complex(1/d.TC())), t)
expr
Out[121]:
In [122]:
expr.evalf()
Out[122]:
In [123]:
def canonicalize(expr):
n, d = expr.as_numer_denom()
return (n/Poly(d).TC())/(d/Poly(d).TC())
In [124]:
expr1 = canonicalize(cancel(Rational(1, 3) + Rational(5, 2)/(x - 2 + I) + Rational(5, 2)/(x - 2 - I)))
expr1.evalf()
Out[124]:
In [125]:
expr2 = canonicalize(cancel(Rational(1/3) + Rational(5, 2)/(x - 2 + I) + Rational(5, 2)/(x - 2 - I)))
expr2.evalf()
Out[125]:
In [126]:
expr2 = canonicalize(cancel(Rational(1/3) + Rational(5/2)/(x - Rational(2 + 0.000001) + I) + Rational(5, 2)/(x - 2 - I)))
expr2.evalf()
Out[126]:
In [127]:
expr1 = canonicalize(cancel(Rational(1, 300000000000000) + Rational(5, 2)/(x - 2 + I) + Rational(5, 2)/(x - 2 - I)))
expr1.evalf()
Out[127]:
In [128]:
expr2 = canonicalize(cancel(Rational(1/300000000000000) + Rational(5, 2)/(x - 2 + I) + Rational(5, 2)/(x - 2 - I)))
expr2.evalf()
Out[128]:
In [129]:
correct_rat_func14.as_numer_denom()[1]
Out[129]:
In [130]:
roots = intervals(correct_rat_func14_rational.subs(t, -t).as_numer_denom()[1], all=True)
roots
Out[130]:
In [255]:
list(map(lambda i: i.evalf(), roots[1][0][0]))
Out[255]:
In [256]:
list(map(lambda i: i.evalf(), roots[1][1][0]))
Out[256]:
In [133]:
intervals(x**3 + 3, all=True, eps=1e-100)
Out[133]:
In [134]:
solve(x**3 + 3.)
Out[134]:
In [135]:
minpoly(sqrt(2) + sqrt(3))
Out[135]:
In [136]:
solve(x**4 - 10*x**2 + 1)
Out[136]:
In [137]:
minpoly(sqrt(-2*sqrt(6) + 5))
Out[137]:
In [138]:
roots10 = intervals(correct_rat_func14_rational.subs(t, -t).as_numer_denom()[1], all=True, eps=1e-14)
roots10
Out[138]:
In [258]:
list(map(lambda i: i.evalf(17), roots10[1][0][0]))
Out[258]:
In [140]:
plot(re(correct_partial_frac14.subs(t, -t)), (t, 1e-4, 1e-2), xscale='log', adaptive=False, nb_of_points=10000)
Out[140]:
In [141]:
plot(exp(-t), (t, 1e-4, 1e-2), xscale='log', adaptive=False, nb_of_points=10000)
Out[141]:
In [142]:
plot(exp(-t) - 2*re(correct_partial_frac14.subs(t, -t)), (t, 1e-4, 1e-2), xscale='log', adaptive=False, nb_of_points=10000)
Out[142]:
In [143]:
alpha1
Out[143]:
In [144]:
alpha1*2
Out[144]:
In [145]:
srepr(correct_rat_func14)
Out[145]:
In [146]:
incorrect_rat_func14 = Mul(Add(Mul(Float('4.1609826642376610641e-28', prec=17), Pow(Symbol('t', real=True), Integer(14))), Mul(Integer(-1), Float('3.7794523874503297816e-24', prec=17), Pow(Symbol('t', real=True), Integer(13))), Mul(Float('5.6743825539523504126e-21', prec=17), Pow(Symbol('t', real=True), Integer(12))), Mul(Integer(-1), Float('3.3490280333667534845e-18', prec=17), Pow(Symbol('t', real=True), Integer(11))), Mul(Float('1.0310151365350495889e-15', prec=17), Pow(Symbol('t', real=True), Integer(10))), Mul(Integer(-1), Float('1.9038640942835345932e-13', prec=17), Pow(Symbol('t', real=True), Integer(9))), Mul(Float('2.2849515765300153532e-11', prec=17), Pow(Symbol('t', real=True), Integer(8))), Mul(Integer(-1), Float('1.8710255961089454423e-9', prec=17), Pow(Symbol('t', real=True), Integer(7))), Mul(Float('1.0753202054485226852e-7', prec=17), Pow(Symbol('t', real=True), Integer(6))), Mul(Integer(-1), Float('4.393280849251123483e-6', prec=17), Pow(Symbol('t', real=True), Integer(5))), Mul(Float('0.00012734070715233181928', prec=17), Pow(Symbol('t', real=True), Integer(4))), Mul(Integer(-1), Float('0.0025674439819028619519', prec=17), Pow(Symbol('t', real=True), Integer(3))), Mul(Float('0.034346984175671474437', prec=17), Pow(Symbol('t', real=True), Integer(2))), Mul(Integer(-1), Float('0.27495604296300041325', prec=17), Symbol('t', real=True)), Float('0.99999999999998168132', prec=17)), Pow(Add(Mul(Float('2.2710680218891296266e-14', prec=17), Pow(Symbol('t', real=True), Integer(14))), Mul(Float('1.8921838854022448175e-13', prec=17), Pow(Symbol('t', real=True), Integer(13))), Mul(Float('1.2922802705792779525e-11', prec=17), Pow(Symbol('t', real=True), Integer(12))), Mul(Float('2.2472030400428528176e-10', prec=17), Pow(Symbol('t', real=True), Integer(11))), Mul(Float('4.836180087864828391e-9', prec=17), Pow(Symbol('t', real=True), Integer(10))), Mul(Float('8.0163997982357503177e-8', prec=17), Pow(Symbol('t', real=True), Integer(9))), Mul(Float('1.1820291576355771705e-6', prec=17), Pow(Symbol('t', real=True), Integer(8))), Mul(Float('0.000014908234724800242013', prec=17), Pow(Symbol('t', real=True), Integer(7))), Mul(Float('0.00016019211440612189201', prec=17), Pow(Symbol('t', real=True), Integer(6))), Mul(Float('0.0014405527154587316474', prec=17), Pow(Symbol('t', real=True), Integer(5))), Mul(Float('0.010574049161691155552', prec=17), Pow(Symbol('t', real=True), Integer(4))), Mul(Float('0.060968185127283594515', prec=17), Pow(Symbol('t', real=True), Integer(3))), Mul(Float('0.25939094125018014036', prec=17), Pow(Symbol('t', real=True), Integer(2))), Mul(Float('0.725043', prec=17), Symbol('t', real=True)), Float('1.0', prec=17)), Integer(-1)))
In [147]:
plot(exp(-t) - correct_rat_func14, (t, 1e-4, 1e2), xscale='log', adaptive=False, nb_of_points=10000)
Out[147]:
From "Efficient Solution of Parabolic Equations by Krylov Approximation Methods" by E. Gallopoulos and Y. Saad Appendix B Table 5
In [148]:
real_part_alternate14 = """
-0.562314417475317895E+01
-0.508934679728216110E+01
-0.399337136365302569E+01
-0.226978543095856366E+01
0.208756929753827868E+00
0.370327340957595652E+01
0.889777151877331107E+01
0.557503973136501826E+02
-0.938666838877006739E+02
0.469965415550370835E+02
-0.961424200626061065E+01
0.752722063978321642E+00
-0.188781253158648576E-01
0.143086431411801849E-03
0.183216998528140087E-11
"""
In [149]:
theta1ra, theta2ra, theta3ra, theta4ra, theta5ra, theta6ra, theta7ra, alpha1ra, alpha2ra, alpha3ra, alpha4ra, alpha5ra, alpha6ra, alpha7ra, alpha0ra = map(Float, real_part_alternate14.strip().split())
theta1ra, theta2ra, theta3ra, theta4ra, theta5ra, theta6ra, theta7ra, alpha1ra, alpha2ra, alpha3ra, alpha4ra, alpha5ra, alpha6ra, alpha7ra, alpha0ra
Out[149]:
In [150]:
imaginary_part_alternate14 = """
0.119406921611247440E+01
0.358882439228376881E+01
0.600483209099604664E+01
0.846173881758693369E+01
0.109912615662209418E+02
0.136563731924991884E+02
0.166309842834712071E+02
-0.204295038779771857E+03
0.912874896775456363E+02
-0.116167609985818103E+02
-0.264195613880262669E+01
0.670367365566377770E+00
-0.343696176445802414E-01
0.287221133228814096E-03
0
"""
In [151]:
theta1ja, theta2ja, theta3ja, theta4ja, theta5ja, theta6ja, theta7ja, alpha1ja, alpha2ja, alpha3ja, alpha4ja, alpha5ja, alpha6ja, alpha7ja, alpha0ja = [Float(i)*I for i in imaginary_part_alternate14.strip().split()]
theta1ja, theta2ja, theta3ja, theta4ja, theta5ja, theta6ja, theta7ja, alpha1ja, alpha2ja, alpha3ja, alpha4ja, alpha5ja, alpha6ja, alpha7ja, alpha0ja
Out[151]:
In [152]:
alpha1ja
Out[152]:
In [153]:
theta1a, theta2a, theta3a, theta4a, theta5a, theta6a, theta7a, alpha1a, alpha2a, alpha3a, alpha4a, alpha5a, alpha6a, alpha7a, alpha0a =\
(r + j for r, j in zip(
(theta1ra, theta2ra, theta3ra, theta4ra, theta5ra, theta6ra, theta7ra, alpha1ra, alpha2ra, alpha3ra, alpha4ra, alpha5ra, alpha6ra, alpha7ra, alpha0ra),
(theta1ja, theta2ja, theta3ja, theta4ja, theta5ja, theta6ja, theta7ja, alpha1ja, alpha2ja, alpha3ja, alpha4ja, alpha5ja, alpha6ja, alpha7ja, alpha0ja)))
In [154]:
thetasa = theta1a, theta2a, theta3a, theta4a, theta5a, theta6a, theta7a
alphasa = alpha1a, alpha2a, alpha3a, alpha4a, alpha5a, alpha6a, alpha7a
In [155]:
correct_partial_frac_alternate14 = alpha0a + Add(*(alpha/(t - theta) + alpha/(t - conjugate(theta)) for alpha, theta in zip(alphasa, thetasa)))
correct_partial_frac_alternate14
Out[155]:
In [156]:
plot((2*re(correct_partial_frac_alternate14) - alpha0a, (t, 1e2, 1e4)), (exp(-t), (t, 1e-2, 1e4)), xscale='log', adaptive=False, nb_of_points=10000)
Out[156]:
In [157]:
plot((exp(-t) - re(correct_partial_frac_alternate14) + 0*alpha0a, (t, 1e-2, 1e4)), xscale='log', adaptive=False, nb_of_points=10000)
Out[157]:
In [158]:
N = symbols('N')
plot_parametric((*phi.subs(N, 4).as_real_imag(), (x, pi/14, pi*13/14)))
Out[158]:
In [159]:
phi
Out[159]:
In [160]:
!say done
In [161]:
correct_rat_func
Out[161]:
In [162]:
apart(correct_rat_func)
Out[162]:
In [163]:
together(apart(correct_rat_func))
Out[163]:
In [164]:
correct_rat_func2
Out[164]:
In [165]:
apart(correct_rat_func2, full=True).doit()
Out[165]:
In [166]:
canonicalize(cancel(apart(correct_rat_func2, full=True).doit()))
Out[166]:
In [167]:
def quad_coeffs(k):
import numpy
xvals = numpy.linspace(1, 2, k - 1, dtype=complex)
xquad = pi*x/k
lxquad = lambdify(x, xquad, 'numpy')
theta = phi.subs(N, k)
ltheta = lambdify(x, theta, 'numpy')
w_j = phi.subs(N, k).diff(x)
alpha = I/k*exp(theta)*w_j
lalpha = lambdify(x, alpha, 'numpy')
return (ltheta(lxquad(xvals)), lalpha(lxquad(xvals)))
In [168]:
quad_coeffs(2)
Out[168]:
In [169]:
quad_coeffs(14)
Out[169]:
In [170]:
thetas
Out[170]:
In [171]:
import numpy as np
import matplotlib.pyplot as plt
In [175]:
plt.plot(*zip(*[i.as_real_imag() for i in thetas]))
Out[175]:
In [180]:
plt.plot(*zip(*[i.as_real_imag() for i in thetasa]))
Out[180]:
In [179]:
plot_parametric((*phi.subs(N, 14).as_real_imag(), (x, -pi, pi)))
Out[179]:
In [181]:
thetas
Out[181]:
In [182]:
thetasa
Out[182]:
In [184]:
quad_thetas, quad_alphas = quad_coeffs(14)
quad_thetas, quad_alphas
Out[184]:
In [192]:
quad_rational_func14 = 2*re(Add(*(alpha/(t - theta) for alpha, theta in zip(quad_alphas, quad_thetas))))
quad_rational_func14
Out[192]:
In [198]:
plot(quad_rational_func14)
Out[198]:
In [199]:
plot(exp(t) - quad_rational_func14, (t, -100, 0))
Out[199]:
In [200]:
plot((exp(t) - quad_rational_func14), (t, -1000, 0))
Out[200]:
In [241]:
re_correct_partial_frac14 = alpha0 + 2*re(Add(*(alpha/(t - theta) for alpha, theta in zip(alphas, thetas))))
plot(re_correct_partial_frac14)
Out[241]:
In [242]:
plot(exp(t) - re_correct_partial_frac14, (t, -100, 0))
Out[242]:
In [243]:
plot(exp(t) - re_correct_partial_frac14, (t, -1000, 0))
Out[243]:
In [244]:
print(np.pi*np.linspace(1, 2, 14 - 1, dtype=complex)/14)
plot_parametric((*phi.subs(N, 4).as_real_imag(), (x, pi/14, pi*13/14)))
Out[244]:
In [323]:
def quad_coeffs2(k):
import numpy
xvals = numpy.linspace(1, k - 1, k//2, dtype=complex)
xquad = pi*x/k
lxquad = lambdify(x, xquad, 'numpy')
theta = phi.subs(N, k)
ltheta = lambdify(x, theta, 'numpy')
w_j = phi.subs(N, k).diff(x)
alpha = I/k*exp(theta)*w_j
lalpha = lambdify(x, alpha, 'numpy')
return (ltheta(lxquad(xvals)), lalpha(lxquad(xvals)))
quad2_thetas, quad2_alphas = quad_coeffs2(32)
print(quad2_thetas)
print(quad2_alphas)
quad2_rational_func14 = 2*re(Add(*(alpha/(t - theta) for alpha, theta in zip(quad2_alphas, quad2_thetas))))
plot(quad2_rational_func14)
plot(exp(t) - quad2_rational_func14, (t, -100, -10))
Out[323]:
In [216]:
plot((re(exp(pi/14*phi.subs(N, 14))), (x, -5, 5)), (im(exp(pi/14*phi.subs(N, 14))), (x, -5, 5)))
Out[216]:
In [226]:
plot(exp(t) - re_correct_partial_frac14, (t, -100, 0))
Out[226]:
In [228]:
plot(exp(t) - correct_rat_func14.subs(t, -t), (t, -100, 0))
Out[228]:
In [245]:
re_correct_partial_frac_alternate14 = alpha0a + 2*re(Add(*(alpha/(t - theta) for alpha, theta in zip(alphasa, thetasa))))
re_correct_partial_frac_alternate14 = re_correct_partial_frac_alternate14.subs(t, -t)
plot(re_correct_partial_frac_alternate14)
Out[245]:
In [246]:
plot(re_correct_partial_frac_alternate14, re_correct_partial_frac14)
Out[246]:
In [247]:
plot((exp(t) - re_correct_partial_frac14, (t, -100, 0)), (exp(t) - re_correct_partial_frac_alternate14, (t, -100, 0)))
Out[247]:
In [248]:
alpha0a
Out[248]:
In [249]:
alphasa
Out[249]:
In [250]:
thetasa
Out[250]:
In [251]:
real_part_alternate14.strip().split()
Out[251]:
In [253]:
print(real_part_alternate14)
In [296]:
def my_theta_alphas14():
num, den = correct_rat_func14_rational.subs(t, -t).as_numer_denom()
# Note, eps is NOT the precision. It's the length of the interval.
# If a root is small (say, on the order of 10**-N), then eps will need to be 10**(-N - d)
# to get d digits of precision. For the order 14 approximation, the roots (thetas) are all
# within order 10**-1...10**1, so eps is *roughly* the precision.
roots = intervals(den, all=True, eps=1e-14)[1]
# eps ought to be small enough that either side of the interval is the precision we want,
# but for the sake of with smaller eps, take the average (center of the rectangle)
# XXX: Make sure to change the evalf precision if eps is lowered.
thetas = [((i + j)/2).evalf(17) for ((i, j), _) in roots]
error = [(j - i).evalf(17) for ((i, j), _) in roots]
alphas = []
for theta in thetas:
q, r = div(den, t - theta)
alpha = (num/q).evalf(17, subs={t: theta})
alphas.append(alpha)
alpha0 = LC(num)/LC(den)
return thetas, alphas, alpha0
In [ ]:
my_thetas, my_alphas, my_alpha0 = my_theta_alphas14()
my_thetas, my_alphas, my_alpha0
In [269]:
my_thetas[0]
Out[269]:
In [272]:
num, den = correct_rat_func14.subs(t, -t).as_numer_denom()
(num.subs(t, my_thetas[0])/cancel(den/(t - my_thetas[0])).subs(t, my_thetas[0])).evalf(17)
Out[272]:
In [278]:
(den/(t - my_thetas[0])).evalf(17, subs={t: my_thetas[0]})
Out[278]:
In [279]:
num.subs(t, my_thetas[0]).evalf(17)
Out[279]:
In [280]:
(num/(den/(t - my_thetas[0]))).evalf(17, subs={t: my_thetas[0]})
Out[280]:
In [284]:
q, r = div(den, t - my_thetas[0])
q.evalf(subs={t: my_thetas[0]}), r.evalf(subs={t: my_thetas[0]})
Out[284]:
In [285]:
(num/q).evalf(17, subs={t: my_thetas[0]})
Out[285]:
In [291]:
[2*i for i in my_alphas]
Out[291]:
In [293]:
n14pi[-1]/n14qi[-1]
Out[293]:
In [295]:
cancel(correct_rat_func14.subs(t, 1/t))
Out[295]:
In [303]:
canonicalize(cancel(correct_partial_frac14))
Out[303]:
In [307]:
canonicalize(cancel(correct_re_partial_frac14))
Out[307]:
In [ ]: