removed redundant output

33c3982a · Benjamin Berkels · 2344f192 · 33c3982a
Commit 33c3982a authored 4 years ago by Benjamin Berkels
--- a/exp-355.ipynb
+++ b/exp-355.ipynb
@@ -83,7 +83,6 @@
      "metadata": {},
      "source": [
        "alpha_index = np.argmin(np.abs(discrepancy-epsilon_a))\n",
-        "display(Math(rf\"\\alpha_d = {alphas[alpha_index]}\"))\n",
        "\n",
        "j = np.arange(1, n+1)\n",
        "plt.title(rf\"Tikhonov, $\\alpha_d$ = {alphas[alpha_index]}\")\n",

 %% Cell type:markdown id: tags:
 ### IPython notebook for Example 3.5.5 from the lecture
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib_inline.backend_inline import set_matplotlib_formats
 from IPython.display import display, Math
 set_matplotlib_formats('svg')
 def k(s, t, d):
    return d / np.power(d**2 + (s-t)**2, 3/2)
 d = 0.1
 n = 80
 h = 1/n
 display(Math(r'\text{Create }A\text{ from Example 1.5.1, }\xi_j=\begin{cases}2&1\leq j\leq40\\1&41\leq j\leq80\end{cases},\ y=A\xi'))
 s = (np.arange(n) + 0.5) * h
 t = (np.arange(n) + 0.5) * h
 A = h * k(s[:, np.newaxis], t[np.newaxis, :], d)
 xi = np.ones(n)
 xi[:n//2] = 2
 y = np.matmul(A, xi)
 ```
 %% Cell type:code id: tags:
 ``` python
 display(Math(r'\tilde y= y+\delta y'))
 np.random.seed(0)
 y_tilde = y + 0.02*(np.random.rand(n)-0.5)
 epsilon_a = np.linalg.norm(y_tilde-y)
 display(Math(rf"\epsilon_a=\|\tilde y-y\|_2 = {epsilon_a:.3f}"))
 ```
 %% Cell type:code id: tags:
 ``` python
 discrepancy = np.zeros(100)
 alphas = np.zeros(100)
 xi_alphas = np.zeros((100, n))
 for i in range(1, 101):
    alpha = 0.005*i
    A_Tikh = np.concatenate((A, alpha*np.eye(n)), axis=0)
    y_tilde_Tikh = np.concatenate((y_tilde, np.zeros(n)), axis=0)
    xi_alpha = np.linalg.lstsq(A_Tikh, y_tilde_Tikh, rcond=0)[0]
    xi_alphas[i-1, :] = xi_alpha
    alphas[i-1] = alpha
    discrepancy[i-1] = np.linalg.norm(np.matmul(A, xi_alpha)-y_tilde)
 plt.title("Tikhonov, discrepancy")
 plt.xlabel(r'$\alpha$')
 plt.ylabel(r'$||A\tilde \xi_{\alpha}- \tilde y||_2$')
 plt.plot(alphas, discrepancy, 'k+')
 plt.hlines(epsilon_a, xmin=alphas[0], xmax=alphas[-1])
 plt.show()
 plt.close()
 ```
 %% Cell type:code id: tags:
 ``` python
 alpha_index = np.argmin(np.abs(discrepancy-epsilon_a))
-display(Math(rf"\alpha_d = {alphas[alpha_index]}"))
 j = np.arange(1, n+1)
 plt.title(rf"Tikhonov, $\alpha_d$ = {alphas[alpha_index]}")
 plt.plot(j, xi_alphas[alpha_index, :], 'kx')
 plt.show()
 plt.close()
 ```
 %% Cell type:code id: tags:
 ``` python
 # CGNE
 xi_tilde = np.zeros(n)
 xi_ks = np.zeros((101, n))
 # xi_tilde is zero, so the initial residuum is the input data.
 r = y_tilde.copy()
 d = np.matmul(A.T, r)
 p = d.copy()
 p_norm = np.linalg.norm(p)
 discrepancy = np.zeros(101)
 discrepancy[0] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)
 for k in range(1, 101):
    q = np.matmul(A, d)
    beta = (np.linalg.norm(p)/np.linalg.norm(q))**2
    xi_tilde += beta * d
    r += -beta*q
    p = np.matmul(A.T, r)
    p_norm_new = np.linalg.norm(p)
    gamma = (p_norm_new/p_norm)**2
    p_norm = p_norm_new
    d = p + gamma*d
    discrepancy[k] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)
    xi_ks[k, :] = xi_tilde
 plt.title("CGNE, discrepancy")
 plt.xlabel(r'$k$')
 plt.ylabel(r'$||A\tilde \xi_k- \xi||_2$')
 plt.yscale('log')
 plt.plot(discrepancy, 'k+')
 plt.hlines(epsilon_a, xmin=0, xmax=discrepancy.shape[0])
 plt.show()
 plt.close()
 ```
 %% Cell type:code id: tags:
 ``` python
 k_d = np.argmin(np.abs(discrepancy-epsilon_a))
 plt.title(rf"CGNE, $k_d$ = {k_d}")
 plt.plot(j, xi_ks[k_d, :], 'kx')
 plt.show()
 plt.close()
 ```
 %% Cell type:code id: tags:
 ``` python
 ```