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

improved the IPython notebook for Examples 2.3.2 and 2.3.3

parent 9dc0d8d6
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
exp-232.ipynb
+ 28
− 21
View file @ bdd7dc99
%% Cell type:markdown id: tags:
### IPython notebook for Example 2.3.2 from the lecture
%% Cell type:code id: tags:
``` python
import numpy as np
from IPython.display import display, Math
display(Math(r'''A=\begin{pmatrix} 1 & 1 \\ 0 & 0 \\ 0 & 1 \end{pmatrix}'''))
A = np.array([[1., 1.], [0., 0.], [0., 1.]])
display(Math(r'''y= \begin{pmatrix} 0.01 \\ 1 \\ 0 \end{pmatrix}'''))
y = np.array([0.01, 1., 0.])
display(Math(r'''\delta y= \begin{pmatrix} 0 \\ 0 \\ 0.01 \end{pmatrix}'''))
delta_y = np.array([0., 0., 0.01])
display(Math(r"A=\begin{pmatrix} 1 & 1 \\ 0 & 0 \\ 0 & 1 \end{pmatrix}"))
A = np.array([[1.0, 1.0], [0.0, 0.0], [0.0, 1.0]])
display(Math(r"y= \begin{pmatrix} 0.01 \\ 1 \\ 0 \end{pmatrix}"))
y = np.array([0.01, 1.0, 0.0])
display(Math(r"\delta y= \begin{pmatrix} 0 \\ 0 \\ 0.01 \end{pmatrix}"))
delta_y = np.array([0.0, 0.0, 0.01])
```
%% Cell type:code id: tags:
``` python
print('rank(A) =', np.linalg.matrix_rank(A))
print('cond(A) = {:.2f}'.format(np.linalg.cond(A)))
display(Math(r"\operatorname{rank}(A) = " + f"{np.linalg.matrix_rank(A)}"))
display(Math(r"\operatorname{cond(A)} = " + f"{np.linalg.cond(A):.2f}"))
```
%% 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'P_Ay=(u_1^Ty)u_1+(u_2^Ty)u_2'))
PAy = np.dot(U[:, 0], y)*U[:, 0] + np.dot(U[:, 1], y)*U[:, 1]
display(Math(r'\text{We have }P_Ay=\begin{pmatrix}0.01 \\ 0 \\ 0 \end{pmatrix}\text{, since }\text{range}(A)= \text{span} \{e_1,e_3\}'))
display(Math(r"P_Ay=(u_1^Ty)u_1+(u_2^Ty)u_2"))
PAy = np.dot(U[:, 0], y) * U[:, 0] + np.dot(U[:, 1], y) * U[:, 1]
display(Math(r"\text{We have }P_Ay=\begin{pmatrix}0.01 \\ 0 \\ 0 \end{pmatrix}\text{, since }\text{range}(A)= \text{span} \{e_1,e_3\}"))
```
%% Cell type:code id: tags:
``` python
print('||P_Ay||/||y|| = {:.2f}'.format(np.linalg.norm(PAy)/np.linalg.norm(y)))
data_error = np.linalg.norm(delta_y)/np.linalg.norm(y)
print('||delta y||/||y|| = {:.2f}'.format(data_error))
display(Math(r"\frac{\|P_Ay\|_2}{\|y\|_2} = " + f"{np.linalg.norm(PAy) / np.linalg.norm(y):.2f}"))
data_error = np.linalg.norm(delta_y) / np.linalg.norm(y)
display(Math(r"\frac{\|\delta y\|_2}{\|y\|_2} = " + f"{data_error:.2f}"))
```
%% Cell type:code id: tags:
``` python
A_dagger = np.linalg.pinv(A)
display(Math(r'''\xi = A^\dagger y'''))
display(Math(r"\xi = A^\dagger y"))
xi = np.matmul(A_dagger, y)
display(Math(r'''\tilde\xi = A^\dagger (y+\tilde y)'''))
display(Math(r"\tilde\xi = A^\dagger (y+\delta y)"))
xi_tilde = np.matmul(A_dagger, y + delta_y)
```
%% Cell type:code id: tags:
``` python
solution_error = np.linalg.norm(xi_tilde-xi)/np.linalg.norm(xi)
print('||xi_tilde - xi||/||xi|| = {:.2f}'.format(solution_error))
print('Q(y,delta_y) = {:.2f}'.format(solution_error/data_error))
solution_error = np.linalg.norm(xi_tilde - xi) / np.linalg.norm(xi)
display(Math(r"\frac{\|\tilde\xi-\xi\|_2}{\|\xi\|_2} = " + f"{solution_error:.2f}"))
display(Math(rf"Q(y,\delta y) = {solution_error / data_error:.2f}"))
```
%% Cell type:code id: tags:
``` python
```
......
exp-233.ipynb
+ 30
− 24
View file @ bdd7dc99
%% Cell type:markdown id: tags:
### IPython notebook for Example 2.3.3 from the lecture
%% Cell type:code id: tags:
``` python
import numpy as np
from IPython.display import display, Math
from sympy import Float, Matrix, latex
display(Math(r'A \in \mathbb{R}^{6 \times 3}, A_1 \in \mathbb{R}^{6 \times 4}\text{ from Example 1.1.1}'))
def my_latex(expr):
return latex(expr, mat_delim='(').replace('.0 &', '&').replace('.0\\', '\\')
display(Math(r"A \in \mathbb{R}^{6 \times 3}, A_1 \in \mathbb{R}^{6 \times 4}\text{ from Example 1.1.1}"))
delta_t = 0.15
n = 6
t = np.arange(1, n+1) * delta_t
A_1 = np.stack((t, np.exp(t), t**3, np.sin(t)), axis=1)
t = np.arange(1, n + 1) * delta_t
A_1 = np.stack((t, np.exp(t), t ** 3, np.sin(t)), axis=1)
A = A_1[:, :-1]
print("A_1 = ")
print(A_1)
print("A = ")
print(A)
display(Math(rf'A_1 = {my_latex(Matrix(np.round(A_1, 8)))}'))
display(Math(rf'A = {my_latex(Matrix(np.round(A, 8)))}'))
```
%% 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)
np.set_printoptions(precision=4, suppress=True)
print("sigma(A) =", sigma)
print("kappa(A) = {:.1f}".format(sigma[0]/sigma[-1]))
display(Math(rf'\Sigma(A) = {my_latex(Matrix(np.round(sigma, 4)))}'))
display(Math(rf'\kappa(A) = {sigma[0] / sigma[-1]:.1f}'))
```
%% Cell type:code id: tags:
``` python
display(Math(r'U_1^TA_1V_1=\Sigma_1'))
display(Math(r"U_1^TA_1V_1=\Sigma_1"))
U_1, sigma_1, VT_1 = np.linalg.svd(A_1)
np.set_printoptions(precision=5)
print("sigma(A_1) =", sigma_1)
print("kappa(A_1) = {:.1e}".format(sigma_1[0]/sigma_1[-1]))
display(Math(rf'\Sigma(A_1) = {my_latex(Matrix(np.round(sigma_1, 5)))}'))
display(Math(r'\kappa(A_1) = ' + latex(Float(sigma_1[0] / sigma_1[-1], 2))))
```
%% Cell type:code id: tags:
``` python
display(Math(r"\delta y^1=10^{-2} \begin{pmatrix} 1 & 0 & -1 & -1 & -0.5 & 1 \end{pmatrix}^T"))
delta_y1 = 1e-2 * np.array([1, 0, -1, -1, -0.5, 1])
display(Math(r"\delta y^2=10^{-2}\begin{pmatrix} -1 & 1 & 1 & -0.5 & -2 & 1 \end{pmatrix}^T"))
delta_y2 = 1e-2 * np.array([-1, 1, 1, -0.5, -2, 1])
```
%% Cell type:code id: tags:
``` python
display(Math(r'\delta y^1=10^{-2} \begin{pmatrix} 1 & 0 & -1 & -1 & -0.5 & 1 \end{pmatrix}^T'))
delta_y1 = 1e-2*np.array([1, 0, -1, -1, -0.5, 1])
display(Math(r'\delta y^2=10^{-2}\begin{pmatrix} -1 & 1 & 1 & -0.5 & -2 & 1 \end{pmatrix}^T'))
delta_y2 = 1e-2*np.array([-1, 1, 1, -0.5, -2, 1])
display(Math(rf'U_1^T\delta y^1 = {my_latex(Matrix(np.round(np.matmul(U_1.T, delta_y1), 4)))}'))
display(Math(rf'U_1^T\delta y^2 = {my_latex(Matrix(np.round(np.matmul(U_1.T, delta_y2), 4)))}'))
display(Math(r"\Rightarrow\ \delta y^1\text{ is dominated by }u_3,\ \delta y^2\text{ dominated by }u_4"))
```
%% Cell type:code id: tags:
``` python
np.set_printoptions(precision=4)
print("U_1^T\delta y^1 =", np.matmul(U_1.T, delta_y1))
print("U_1^T\delta y^2 =", np.matmul(U_1.T, delta_y2))
display(Math(r'\Rightarrow\ \delta y^1\text{ is dominated by }u_3, \delta y^2\text{ dominated by }u_4'))
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment