Skip to content
Snippets Groups Projects
Commit 90ae0ca7 authored by Benjamin Berkels's avatar Benjamin Berkels
Browse files

improved the IPython notebook for Example 2.2.4

parent 4e9d6f32
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 34 additions and 24 deletions
exp-224.ipynb
+ 34
− 24
View file @ 90ae0ca7
%% Cell type:markdown id: tags:
### IPython notebook for Example 2.2.4 from the lecture
%% Cell type:code id: tags:
``` python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats, display, Math
set_matplotlib_formats('svg')
from matplotlib_inline.backend_inline import set_matplotlib_formats
from IPython.display import display, Math
n = 40
h = 1/n
set_matplotlib_formats("svg")
display(Math(r'''A = h \begin{pmatrix} 1 & & & \\ 1 & 1 & &0 \\ \vdots &\vdots & \ddots & \\ 1 & 1& \dots & 1 \end{pmatrix}'''))
A = h*np.tril(np.ones((n, n), dtype=int), 0)
n = 40
h = 1 / n
display(Math(r'''y= \frac{1}{n}(1,1,\ldots,1)^T'''))
y = np.ones(n)/n
display(Math(r'''\delta y = \frac{1}{n}(1,-1,1,-1, \ldots)^T'''))
display(Math(r"\text{Setting from Section 1.2:}"))
display(Math(r"A = h \begin{pmatrix} 1 & & & \\ 1 & 1 & &0 \\ \vdots &\vdots & \ddots & \\ 1 & 1& \dots & 1 \end{pmatrix}\text{ (discrete integration matrix)}"))
A = h * np.tril(np.ones((n, n), dtype=int), 0)
display(Math(r"y= \frac{1}{n}(1,1,\ldots,1)^T\text{ (low frequency input)}"))
y = np.ones(n) / n
display(Math(r"\delta y = \frac{1}{n}(1,-1,1,-1, \ldots)^T\text{ (high frequency pertubation)}"))
delta_y = y.copy()
delta_y[1::2] = -1/n
display(Math(r'''\Rightarrow \|y\|_2= \|\delta y\|_2= \frac{1}{\sqrt{n}}'''))
delta_y[1::2] = -1 / n
display(Math(r"\Rightarrow \|y\|_2= \|\delta y\|_2= \frac{1}{\sqrt{n}}"))
```
%% Cell type:code id: tags:
``` python
display(Math(r'U^TAV=\Sigma'))
display(Math(r"U^TAV=\Sigma"))
U, sigma, VT = np.linalg.svd(A)
display(Math(r'\hat y=U^Ty'))
display(Math(r"\hat y=U^Ty"))
y_hat = np.matmul(U.T, y)
display(Math(r'\hat{\delta y}:= U^T \delta y'))
display(Math(r"\hat{\delta y}:= U^T \delta y"))
delta_y_hat = np.matmul(U.T, delta_y)
sigma_inv = 1/sigma
sigma_inv = 1 / sigma
j = np.arange(1, 41)
```
%% Cell type:code id: tags:
``` python
plt.yscale('log')
plt.plot(j, np.abs(y_hat), '.', label=r'$|(U^Ty)_j|$')
plt.plot(j, np.abs(delta_y_hat), '+', label=r'$|(U^T \delta y)_j|$')
plt.yscale("log")
plt.plot(j, np.abs(y_hat), ".", label=r"$|(U^Ty)_j|$")
plt.plot(j, np.abs(delta_y_hat), "+", label=r"$|(U^T \delta y)_j|$")
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
plt.close()
plt.yscale('log')
plt.plot(j, sigma_inv, '*', label=r'$\sigma_j^{-1}$')
plt.yscale("log")
plt.plot(j, sigma_inv, "*", label=r"$\sigma_j^{-1}$")
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
plt.close()
plt.yscale('log')
plt.plot(j, sigma_inv*np.abs(y_hat), '.', label=r'$\sigma_j^{-1}|(U^Ty)_j|$')
plt.plot(j, sigma_inv*np.abs(delta_y_hat), '+', label=r'$\sigma_j^{-1}|(U^T \delta y)_j|$')
plt.yscale("log")
plt.plot(j, sigma_inv * np.abs(y_hat), ".", label=r"$\sigma_j^{-1}|(U^Ty)_j|$")
plt.plot(j, sigma_inv * np.abs(delta_y_hat), "+", label=r"$\sigma_j^{-1}|(U^T \delta y)_j|$")
plt.legend()
plt.show()
display(Math(r'''\Rightarrow\ \|A^{-1} y\|_2 = \|\Sigma^{-1} \hat{y}\|_2 \ll \|A^{-1} \delta y\|_2= \|\Sigma^{-1} \hat{\delta y}\|_2\ \Rightarrow\ Q(y,\delta y) \gg 1'''))
display(Math(r"\Rightarrow\ \|A^{-1} y\|_2 = \|\Sigma^{-1} \hat{y}\|_2 \ll \|A^{-1} \delta y\|_2= \|\Sigma^{-1} \hat{\delta y}\|_2\ \Rightarrow\ Q(y,\delta y) \gg 1"))
```
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment