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

improved the IPython notebook for Example 3.4.5

parent 8b29d296
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
### IPython notebook for Example 3.4.5 from the lecture
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats, display, Math
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.7.1, }\xi_j=\begin{cases}2&1\leq j\leq40\\1&41\leq j\leq80\end{cases},\ y=A\xi'))
A = np.zeros((n, n))
for i in range(n):
for j in range(n):
A[i, j] = h * k((i+0.5)*h, (j+0.5)*h, d)
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)
j = np.arange(1, n+1)
plt.title(r"$xi$")
plt.title(r"$\xi$")
plt.plot(j, xi, 'kx')
plt.show()
plt.close()
plt.title(r"$y=A\xi$")
plt.plot(j, y, 'kx')
plt.show()
plt.close()
```
%% 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)
print('||y_tilde-y||/||y|| = {:.4f}'.format(np.linalg.norm(y_tilde-y)/np.linalg.norm(y)))
display(Math(r"\frac{\|\tilde y-y\|_2}{\|y\|_2} = " + f"{np.linalg.norm(y_tilde-y)/np.linalg.norm(y):.4f}"))
plt.title(r"$\tilde \xi = $numpy.linalg.solve$(A, \tilde y)$")
plt.plot(j, np.linalg.solve(A, y_tilde), 'k+')
plt.show()
plt.close()
```
%% Cell type:code id: tags:
``` python
error = np.zeros(100)
alphas = np.zeros(100)
for i in range(1, 101):
alpha = 0.005*i
A_Tikh = np.concatenate((A, np.diag(alpha*np.ones(n))), axis=0)
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]
alphas[i-1] = alpha
error[i-1] = np.linalg.norm(xi_alpha-xi)
if i == 1 or i == 8 or i == 100:
plt.title(r"Tikhonov, $\alpha$ = {}".format(alpha))
plt.title(rf"Tikhonov, $\alpha$ = {alpha}")
plt.plot(j, xi_alpha, 'kx')
plt.show()
plt.close()
plt.title("Tikhonov, error")
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$||\tilde \xi_{\alpha}- \xi||_2$')
plt.plot(alphas, error, 'k+')
plt.show()
plt.close()
print('alpha_opt = {}'.format(alphas[np.argmin(error)]))
display(Math(r'\alpha_\text{opt} = ' + f'{alphas[np.argmin(error)]}'))
```
%% Cell type:code id: tags:
``` python
# CGNE
xi_tilde = np.zeros(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)
error = np.zeros(201)
error[0] = np.linalg.norm(xi_tilde-xi)
for k in range(1, 201):
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
if k == 5 or k == 40 or k == 200:
plt.title(r"CGNE, $k$ = {}".format(k))
plt.title(rf"CGNE, $k$ = {k}")
plt.plot(j, xi_tilde, 'kx')
plt.show()
plt.close()
error[k] = np.linalg.norm(xi_tilde-xi)
plt.title("CGNE, error")
plt.xlabel(r'$k$')
plt.ylabel(r'$||\tilde \xi_k- \xi||_2$')
plt.plot(error, 'k+')
plt.show()
plt.close()
print('k_opt = {}'.format(np.argmin(error)))
display(Math(r'k_\text{opt} = ' + f'{np.argmin(error)}'))
```
%% Cell type:code id: tags:
``` python
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment