<jupyter_start><jupyter_text>Iterative methods for solving linear systems Let us consider an elastic cable fixed at the extremes occupying at rest the space $[0,1]$. We apply on it a given force $f(x)$. Its displacement at the point $x$, $u(x)$ is the solution of the following problem:$$-u_h''(x) = f(x)\ \mathrm{ in }\ [0,1]$$$$u_h(0) = 0, \quad u_h(1) = 0$$For the sake of simplicity we consider a uniform discretisation of the interval $[0,1]$ using $n$ equispaced *sample points* to construct a **Finite Difference Approximation**.<jupyter_code>%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
import numpy.linalg

n = 33
h = 1./(n-1)

x=linspace(0,1,n)<jupyter_output><empty_output><jupyter_text>Central finite difference approximationWe consider the approximation$$-u''_i := \frac{-u_{i-1} + 2u_i - u_{i+1}}{h^2}$$and solve the linear system$$\frac{-u_{i-1} + 2u_i - u_{i+1}}{h^2} = f_i \qquad i=1,\ldots,n-2$$This results in a tridiagonal matrix with $n-2$ entries. We assemble all of the entries anyway, and then perform a **line elimination** to make the system solvable. This is typical of problems with **Dirichlet Boundary conditions**.<jupyter_code>a = -ones((n-1,)) # Offdiagonal entries
b = 2*ones((n,)) # Diagonal entries
A = (diag(a, -1) + diag(b, 0) + diag(a, +1))
print(A)

# Now we scale it to 1/h^2
A /= h**2

print(linalg.cond(A))<jupyter_output>[[ 2. -1.  0. ...,  0.  0.  0.]
 [-1.  2. -1. ...,  0.  0.  0.]
 [ 0. -1.  2. ...,  0.  0.  0.]
 ..., 
 [ 0.  0.  0. ...,  2. -1.  0.]
 [ 0.  0.  0. ..., -1.  2. -1.]
 [ 0.  0.  0. ...,  0. -1.  2.]]
467.842628839<jupyter_text>Now let's consider a simple right hand side:$$f_i = x_i(1-x_i)$$which gives as exact solution$$u_{\mathrm{exact}} = \frac{x^4}{12} - \frac{x^3}{6} +\frac{x}{12}$$<jupyter_code>f =x*(1.-x)

# Imposing homogeneous BCs keeping A SPD
A[0,:] = 0
A[:,0] = 0
A[0,0] = 1.

A[:,-1] = 0
A[-1,:] = 0
A[-1,-1] = 1

f[0] = 0
f[-1] = 0

A1 = A.copy()

u = linalg.solve(A1, f)

u_ex = (x**4)/12. - (x**3)/6. + x/12.

_ = plot(x,u,'ro')
_ = plot(x,u_ex)<jupyter_output><empty_output><jupyter_text>Jacobi$$ x_i^{k+1} = \frac{1}{A_{ii}} \times \left(b_i - \sum_{j\neq i} a_{ij}x_j^k\right)$$<jupyter_code>def jacobi(A,b,nmax=10000,eps=1e-10):
    N=len(A)
    x=zeros_like(b)
    x_old=ones_like(b)
    tol=eps+1
    it=0
    while(it<nmax and tol>eps):
        it+=1
        for i in range(N):
            x[i]=1./A[i,i]*(b[i]-dot(A[i,0:i],x_old[0:i]) - dot(A[i,i+1:N],x_old[i+1:N]))
    res= b-dot(A,x)
    tol=numpy.linalg.norm(res,2)
    
    x_old =x.copy()
    
    print it, tol
    return x

sol_jacobi = jacobi(A,f)<jupyter_output>10000 723.996585827<jupyter_text>Gauss-Seidel$$ x_i^{k+1} = \frac{1}{A_{ii}} \times \left(b_i - \sum_{j=0}^{i-1} a_{ij}x_j^{k+1} - \sum_{j=i+1}^{N} a_{ij}x_j^k\right)$$<jupyter_code>def gauss_seidel(A,b,nmax=10000,eps=1e-10):
    N=len(A)
    x=zeros_like(b)
    x_old=ones_like(b)
    tol=eps+1
    it=0
    while(it<nmax and tol>eps):
        it+=1
        for i in range(N):
            x[i]=1./A[i,i]*(b[i]-dot(A[i,0:i],x[0:i]) - dot(A[i,i+1:N],x_old[i+1:N]))
        
        res= b-dot(A,x)
        tol=numpy.linalg.norm(res,2)
    
        x_old =x.copy()
    
    print it, tol
    return x

GS = gauss_seidel(A,f)<jupyter_output>2790 9.99401683743e-11<jupyter_text>Gradient method   $$   {\bf r}^k = {\bf b} - A {\bf x}^k   $$      $$   \alpha^k = \frac{{\bf r}^{k^{T}} {\bf r}^k}{{\bf r}^{k^{T}} A{\bf r}^k}   $$      $$   {\bf x}^{k+1} = {\bf x}^k + \alpha^k {\bf r}^k   $$       preconditioned      $$   {\bf z}^k = P^{-1} {\bf r}^k   $$      $$   \alpha^k = \frac{{\bf z}^{k^{T}} {\bf r}^k}{{\bf z}^{k^{T}} A{\bf z}^k}   $$      $$   \bf x^{k+1} =\bf x^{k}+\alpha^{k}\bf z^{k}   $$  $$  $$<jupyter_code>def gradient(A,b,P,nmax=8000,eps=1e-10):
    n=len(A)
    z=zeros_like(b)
    tol=eps+1
    r=b-dot(A,x)
    it=0
    while(it<nmax and tol>eps):
        z=linalg.solve(P,r)
        alpha=dot(r,z)/(dot(z,dot(A,z)))
        x+=z*alpha
        r-=dot(A,z)*alpha
        tol=numpy.linalg.norm(r,2)
        it+=1
    print it, tol

gradient(A,f,identity(len(A)))
print ""
gradient(A,f,A)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2264331974094507, 'fraction_numerical': 0.05037179179659391, 'mean_word_length': 3.028985507246377, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '662513', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2032, 'n_tokens_neox': 1908, 'n_words': 390, 'path': 'python-lectures/10_linear_systems_iterative.ipynb', 'repo_name': 'denocris/Introduction-to-Numerical-Analysis', 'size': '4169'}