Stability (hard to prove for nontrivial equations): $\|A_h^{-1}\| \leq C$ as $h\to 0$.
Then $\|u_h - u\| \leq C \|A_h u_h - A_h u\| = C \|f - A_h u\| = C \epsilon_h~~~$ ($=O(h^2)$ in our case)
Exact problem: $ u'' = f, \quad u(0)=u(1) = 0 $
Approximate problem: $$ \begin{align*} h^{-2} (u _h(x+h)-2u _h(x)+u _h(x-h)) &= f(x), \\ u _h(0) = u _h(1) &= 0 \end{align*} $$ on a grid with $h=1/N$
Approximation (we need to substitute $u$ into the approximate problem): $$ \begin{align*} & h^{-2} (u(x+h)-2u(x)+u(x-h)) - f(x) \\&= u''(x) - f(x) + O(h^2) = O(h^2) \end{align*} $$ (easy!)
Stability... gotta work hard...
$\|A_h^{-1}\| \leq C$ $\Leftrightarrow$ $|\lambda_\min(A_h)| \geq C^{-1}$, where $\lambda_\min$ is the eigenvalue with the smallest absolute value
Exercise: Do some math:
In a tabular form we write it as $$ k \Delta_h \sim \begin{pmatrix} 0 & k & 0 \\ k & -4k & k \\ 0 & k & 0 \end{pmatrix} $$
If $k\ne{\rm const}$ then we shall be more careful: $$\scriptstyle {\rm div}(k\, \nabla u)(x,y) = \left(\begin{array}{ccc} 0 & k(x,y-h/2) & 0 \\ k(x-h/2,y) & \text{guess} & k(x+h/2,y) \\ 0 & k(x,y+h/2) & 0 \end{array}\right) $$
Let's see how it works in 1D:
This is an explicit discretization: we don't need to invert $\Delta_h$ (also called the "forward difference scheme")
Exercise: Order of approximation: $O(\tau + h^2)$
This is called the "Crank-Nicolson" scheme
Here we do need to invert $\Delta_h$: $$ u_h(t+\tau) = (I - \frac{1}{2} \tau \Delta_h)^{-1} (I + \frac{1}{2} \tau \Delta_h) u_h(t) $$
Exercise: approximation error is $O(\tau^2 + h^2)$
Exercise: stability is unconditional (i.e., always stable!)
If we can efficiently invert $(I - \frac{1}{2} \tau \Delta_h)^{-1}$ then it is better than the explicit!
In [1]:
from IPython.html.widgets import interact
from PIL import Image
from numba import autojit
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
img=Image.open('fig/photo.png')
imarray_orig=np.array(img)
img=Image.open('fig/photo_noisy.png')
imarray=np.array(img)
imarray_noisy=imarray;
import PIL
from IPython.display import display, Image
from io import BytesIO
def display_img_array(ima):
im = PIL.Image.fromarray(ima)
bio = BytesIO()
im.save(bio, format='png')
display(Image(bio.getvalue(), format='png', retina=True))
import pyamg
A = sp.sparse.csr_matrix(pyamg.gallery.laplacian.poisson(imarray.shape))
imarray = imarray.astype(float)
vec = imarray.reshape(-1,1);
for iter in range(1, 20):
vec += -0.1*A*vec;
imarray = vec.reshape(imarray.shape);
imarray_blur = imarray.astype(dtype="uint8")
def set_cursor(i):
if i==0:
display_img_array(imarray_orig)
elif i==1:
display_img_array(imarray_noisy.astype(dtype="uint8"))
else:
display_img_array(imarray_blur)
interact(set_cursor, i=(0, 2, 1))
#display_img_array(imarray_orig)
#display_img_array(imarray_noisy.astype(dtype="uint8"))
#display_img_array(imarray_blur)
In [2]:
imarray.shape
Out[2]:
Let's solve $$ \begin{align*} u' &= {\rm div} (k \nabla u) \qquad \text{for }(x,y)\in\Omega \\ u|_{t=0} &= u_0 \\ u|_{\Gamma} &= u_0|_{\Gamma} \end{align*}, $$ where $k=k(\nabla u) = 1/\sqrt{1+|\nabla u|^2}$ ($|v|$ is the length of $v$)
We use the explicit finite difference method:
We discretize it with $\tau=0.25$ below:
In [2]:
imarray = imarray_noisy.astype(float)
def time_step(imarray):
new_imarray = imarray;
for i in range(1,imarray.shape[0]-2):
for j in range(1,imarray.shape[1]-2):
dx1 = (imarray[i+1,j]-imarray[i,j]);
dx0 = (imarray[i,j]-imarray[i-1,j]);
dx = (dx1+dx0)/2;
dy1 = (imarray[i,j+1]-imarray[i,j]);
dy0 = (imarray[i,j]-imarray[i,j-1]);
dy = (dy0+dy1)/2;
new_imarray[i,j] += 0.25*(dx1/np.sqrt(1+dx1*dx1+dy*dy)-dx0/np.sqrt(1+dx0*dx0+dy*dy));
new_imarray[i,j] += 0.25*(dy1/np.sqrt(1+dy1*dy1+dx*dx)-dy0/np.sqrt(1+dy0*dy0+dx*dx));
imarray = new_imarray;
time_step = autojit(time_step)
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.animation as animation
num_steps = 20
num_images = 10
imarray_denoise = []
for outer_iter in xrange(num_images): # how many images will be produced
for iter in xrange(num_steps): # how many time steps to do for each image
time_step(imarray)
imarray_denoise.append(imarray.astype(dtype="uint8"))
#plt.imshow(imarray_denoise)
#display_img_array(imarray_denoise)
def set_cursor(i):
display_img_array(imarray_denoise[i])
#display_img_array(imarray_orig)
interact(set_cursor, i=(0, num_images-1, 1))
Out[2]:
In order to view animation of noise removal you need to install JSAnimation library. The easiset way to do this is to run from command line:
conda install --channel https://conda.binstar.org/IOOS jsanimation
Win-64 users:
conda install --channel https://conda.binstar.org/melund JSAnimation
You can explore all available packages of that library by running:
binstar search -t conda jsanimation
In [3]:
from JSAnimation import IPython_display
fig = plt.figure(figsize=(12, 9))
plt.axis('off')
im = plt.imshow(imarray_denoise[0], cmap='gray')
def animate(i):
im.set_data(imarray_denoise[i])
return [im]
def init():
im.set_data(imarray_denoise[0])
return im
animation.FuncAnimation(fig, animate, frames=num_images,
init_func=init,interval=300, blit=False)
Out[3]:
Idea:
In [4]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/alex.css", "r").read()
return HTML(styles)
css_styling()
Out[4]:
In [ ]: