Finish first question

_code/pset_10/py/pade_approximants.py

+import matplotlib.pyplot as plt
+import numpy as np
+import sympy as sp
+
+x = sp.symbols("x", real=True)
+
+def pade_approximant(function, N, M):
+    # Compute relevant taylor series terms
+    derivative = function
+    taylor_coefficients = [ function.subs(x, 0) ]
+    for i in range(1, N + M + 1):
+        derivative = sp.diff(derivative, x)
+        taylor_coefficients.append(derivative.subs(x, 0) / sp.factorial(i))
+
+    # Build matrix
+    one_one = [1] + [0] * M
+    matrix_rows = [one_one]
+    for i in range(1, M + 1):
+        new_row = []
+        for j in range(0, min(N + i, M) + 1):
+            new_row.append(taylor_coefficients[N + i - j])
+        for j in range(min(N + i, M) + 1, M + 1):
+            new_row.append(0)
+        matrix_rows.append(new_row)
+    A = sp.Matrix(matrix_rows)
+    b = sp.Matrix(one_one)
+    sp.pprint(A)
+    sp.pprint(b)
+
+    # Solve
+    answer = sp.linsolve((A, b))
+    b_coeffs = list(answer.args[0])
+
+    a_coeffs = [
+        sum(b_coeffs[m] * taylor_coefficients[n - m] for m in range(0, min(n, M) + 1))
+        for n in range(0, N + 1)
+    ]
+
+    print(a_coeffs)
+    print(b_coeffs)
+    print("")
+
+    return lambda x_val: sum(a_coeffs[n] * x_val ** n for n in range(0, N + 1)) \
+        / sum(b_coeffs[m] * x_val ** m for m in range(0, M + 1))
+
+# print pade approximant values
+function = sp.exp(x)
+pade_approximations = [ pade_approximant(function, i, i) for i in range(1, 6) ]
+for approx in pade_approximations:
+    print(approx(1))
+print("")
+for approx in pade_approximations:
+    print(approx(1.0))
+print("")
+pade_errors = [ abs(approx(1.0) - function.subs(x, 1.0)) for approx in pade_approximations ]
+
+# Create polynomial approximations and print their values
+poly_approximations = [
+    #lambda x_val: sum(taylor_coefficients[n] * x**n for n in range(0, order)) for order in [3, 5, 7, 9, 11]
+    sum(x**n / sp.factorial(n) for n in range(0, order)) for order in [3, 5, 7, 9, 11]
+]
+for approx in poly_approximations:
+    print(approx)
+    print(approx.subs(x, 1))
+    print(approx.subs(x, 1.0))
+    print("")
+poly_errors = [ abs(approx.subs(x, 1.0) - function.subs(x, 1.0)) for approx in poly_approximations ]
+
+# Graph errors
+fig1 = plt.figure()
+left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
+ax1 = fig1.add_axes([left, bottom, width, height])
+x_vals = [3, 5, 7, 9, 11] # number of free parameters
+ax1.set_yscale("log")
+ax1.plot(x_vals, pade_errors, label="padé")
+ax1.plot(x_vals, poly_errors, label="poly")
+ax1.set_xlabel("free parameters")
+ax1.set_ylabel("absolute error")
+ax1.legend()
+ax1.set_title("Absolute errors of Padé and polynomial approximations")
+fig1.savefig("../../../assets/img/10_errors.png", transparent=True)

_psets/10.md

@@ -96,3 +96,46 @@ equations, and pull new ones from higher order terms (i.e. consider  $$L > N + M
 possible that all additional rows will be degenerate; in this case, the approximant is legitimately
 underdetermined. This means it can match the function exactly, and you can reduce $$M$$ until your
 system is nonsingular.
+
+Ok, so now to actually answer the question. I wrote a SymPy
+[script](https://gitlab.cba.mit.edu/erik/nmm_2020_site/-/tree/master/_code/pset_10/py/pade_approximants.py)
+that uses the strategy we just derived to build approximants for me. For $$f(x) = e^x$$, I get the
+following approximations:
+
+$$
+\begin{aligned}
+[1/1]_ f(x) &= \frac{1 + x/2}{1 - x/2} \\
+[2/2]_ f(x) &= \frac{1 + x/2 + x^2/12}{1 - x/2 + x^2/12} \\
+[3/3]_ f(x) &= \frac{1 + x/2 + x^2/10 + x^3/120}{1 - x/2 + x^2/10 - x^3/120} \\
+[4/4]_ f(x) &= \frac{1 + x/2 + 3x^2/28 + x^3/84 + x^4/1,680}{1 - x/2 + 3x^2/28 - x^3/84 + x^4/1,680} \\
+[5/5]_ f(x) &= \frac{1 + x/2 + x^2/9 + x^3/72 + x^4/1,008 + x^5/30,240}{1 - x/2 + x^2/9 - x^3/72 + x^4/1,008 - x^5/30,240} \\
+\end{aligned}
+$$
+
+These give the corresponding approximations for $$e$$:
+
+$$
+\begin{aligned}
+[1/1]_ f(1) &= 3 &= 3.00000000000000 \\
+[2/2]_ f(1) &= \frac{19}{7} &\approx 2.71428571428571 \\
+[3/3]_ f(1) &= \frac{193}{71} &\approx 2.71830985915493 \\
+[4/4]_ f(1) &= \frac{2,721}{1,001} &\approx 2.71828171828172 \\
+[5/5]_ f(1) &= \frac{49,171}{18,089} &\approx 2.71828182873569
+\end{aligned}
+$$
+
+Polynomial approximations, on the other hand (to equivalent orders), give
+
+$$
+\begin{aligned}
+\sum_{n = 0}^2 \frac{x^n}{n!} &= \frac{5}{2} &= 2.50000000000000 \\
+\sum_{n = 0}^4 \frac{x^n}{n!} &= \frac{65}{24} &\approx 2.70833333333333 \\
+\sum_{n = 0}^6 \frac{x^n}{n!} &= \frac{1,957}{720} &\approx 2.71805555555556 \\
+\sum_{n = 0}^8 \frac{x^n}{n!} &= \frac{109,601}{40,320} &\approx 2.71827876984127 \\
+\sum_{n = 0}^{10} \frac{x^n}{n!} &= \frac{9,864,101}{3,628,800} &\approx 2.71828180114638
+\end{aligned}
+$$
+
+Here are the different errors.
+
+![errors](../../assets/img/10_errors.png)