```python
%matplotlib inline
import numpy as np
from matplotlib.pyplot import *

import sympy as sym
import pandas as pd
```

Consider the following one-dimensional PDE:
$$
-u_{xx}(x) = f(x)\quad\mathrm{ in }\ \Omega = (0, \pi)
$$
$$
u(x) = 0, \quad\mathrm{ on }\ \partial\Omega = \{0, \pi\}
$$

Given the following $4^{th}$ order finite difference approximation of the second order derivative:

$$u_{xx}(x_i) = \frac{-u_{i-2}+16u_{i-1}-30u_i+16u_{i+1}-u_{i+2}}{12h^2}$$

Implement a function that given the domain interval, the forcing function, the number of discretization points, the boundary conditions, returns the matrix $A$ and the the right hand side $b$.


```python
def finDif( omega, f, n, bc ):
    '''
    Given the domain interval omega, 
    the forcing function f, 
    the number of discretization points n,
    the boundary conditions bc, 
    finDif returns the matrix A and the the right hand side b
    '''
    assert n>1
    assert omega[1]>omega[0]
    x=np.linspace(omega[0],omega[1],n)
    h = (omega[1]-omega[0])/(n-1)
    a2 = np.ones((n-2,)) # Offdiagonal entries
    a1 = -16 * np.ones((n-1,)) # Offdiagonal entries
    a0 = 30 * np.ones((n,)) # Diagonal entries
    A = (np.diag(a2, -2) + np.diag(a1, -1) + np.diag(a0, 0 )+ np.diag(a1, 1 )+ np.diag(a2, 2))/(12)
    b=f(x)
    # Change first row of the matrix A
    A[0,:] = 0
    A[:,0] = 0
    A[0,0] = 1
    b[0] = bc[0]/(h**2)
    # Change last row of the matrix A
    A[-1,:] = 0
    A[:,-1] = 0
    A[-1,-1] = 1
    b[-1] = bc[1]/(h**2)
    A = A / ( h**2 )
    return A, b
```

Call the function using:


```python
omega = [0, np.pi]
f = lambda x : np.sin(x)
n=100
bc = [0,0]
A, b = finDif(omega, f, n, bc)
```

Implement two functions that compute the LU and the Cholesky factorization of the system matrix $A$


```python
def LU(A):
    '''
    LU returns the LU factorization of the system matrix A
    '''
    A = A.copy()
    N=len(A)
    for k in range(N-1):
        if (abs(A[k,k]) < 1e-15):
            raise RuntimeError("Null pivot")
        A[k+1:N,k] /= A[k,k]
        for j in range(k+1,N):
            A[k+1:N,j] -= A[k+1:N,k]*A[k,j]
    L=np.tril(A)
    for i in range(N):
        L[i,i]=1.0
    U = np.triu(A)
    return L, U

L, U = LU(A)
```


```python
def cholesky(A):
    '''
    cholesky returns the Cholesky factorization of the system matrix A
    '''
    A = A.copy()
    N = len(A)
    for k in range(N-1):
        A[k,k] = np.sqrt(A[k,k])
        A[k+1:N,k] = A[k+1:N,k]/A[k,k]
        for j in range(k+1,N):
            A[j:N,j] = A[j:N,j] - A[j:N,k]*A[j,k]
    A[-1,-1] = np.sqrt(A[-1,-1])
    L=np.tril(A)
    return L, L.transpose()

HT, H = cholesky(A)
```

Implement forward and backward substitution functions to exploit the developed factorization methods to solve the derived linear system of equations.


```python
def L_solve(L,rhs):
    x = np.zeros_like(rhs)
    N = len(L)
    x[0] = rhs[0]/L[0,0]
    for i in range(1,N):
        x[i] = (rhs[i] - np.dot(L[i, 0:i], x[0:i]))/L[i,i]
    return x
```


```python
def U_solve(U,rhs):
    x = np.zeros_like(rhs)
    N=len(U)
    x[-1] = rhs[-1]/U[-1,-1]
    for i in reversed(range(N)):
        x[i] = (rhs[i] - np.dot(U[i, i+1:N], x[i+1:N]))/U[i,i]
    return x
```

Solve the derived linear system using the implemented functions and plot the computed solution:


```python
x=np.linspace(omega[0],omega[1],n)
# My solution
w = L_solve(L,b)
u = U_solve(U,w)
# Numpy solution
u_ex = np.linalg.solve(A, b)
# Lets plot
_ = plot(x,u,'o-b', label="LU")
_ = plot(x,u_ex,'r', label="numpy")
legend()
show()
print(np.linalg.norm(u-u_ex,2)/np.linalg.norm(u,2))
```

Considering the new domain $\Omega = (0,1)$ and the forcing term $f(x) = x(1-x)$ with B.C. $u(x) = 0$, on $\partial \Omega = {0,1}$ produce a plot and a table where you show the decay of the error w.r.t. the number of grid points.
(The analytical solution for the above problems is $u_{an} = \frac{x^4}{12} - \frac{x^3}{6} + \frac{x}{12}$)


```python
def calculate_error(omega, f, bc, Nmin, Nmax, sol_exact):
    '''
    calculate the error between the LU solution and the function sol_exact
    given:
    the interval omega
    the forcing force f
    the boundary condictions bc
    as a function of n integer in between Nmin and Nmax
    '''
    errors = []
    N=[]
    error = 0.
    for n in range(Nmin,Nmax):
        x=np.linspace(omega[0],omega[1],n)
        A, b = finDif(omega, f, n, bc)
        L, U = LU(A)
        w = L_solve(L,b)
        u = U_solve(U,w)
        #u = np.linalg.solve(A, b)
        u_ex =sol_exact(x)
        error=max(abs(u-u_ex))
        errors.append(error)
        N.append(n)
    return N,errors


sol_exact = lambda x :  (x**4)/12. - (x**3)/6. + x/12
omega2 = [0,1]
f2 = lambda x : x*(1-x)
bc2 = [0,0]
m,errors=calculate_error(omega2, f2, bc2, 5, 101, sol_exact)
semilogy(errors, label="error")
title("Decay of the error w.r.t. the number of grid points")
legend()
df = pd.DataFrame({'n' : m, 'error' : errors})  
display(df)
# df.style
```

Exploit the derived LU factorizations to compute the condition number of the system's matrix $A$ using the original problem formulation.


```python
def IPM(A, x0, mu, tol, nmax ):
    '''
    Power method with shift (or inverse power method):  
    given the matrix A and vector x0, together with the shift mu
    IPM returns eigenvalue close to mu and the relative eigenvector
    with tolerance tol or after nmax iterations 
    '''
    M = A - mu * np.identity(len(A))
    L,U = LU(M)
    err = tol + 1.
    it = 0
    q = x0/(np.linalg.norm(x0,2))
    while it < nmax and err > tol :
        y = L_solve(L,q)
        x = U_solve(U,y)
        q = x/(np.linalg.norm(x,2))
        z = np.dot(A,q)
        l = np.dot(q.T,z)
        err =np.linalg.norm(z-l*q,2)
        it += 1
    return l,q

def PM(A, z0, tol, nmax):
    '''
    Power method:
    given A and z0
    PM returns the maximum eigenvalue and the associated eigenvector
    with tolerance tol or after nmax iterations
    '''
    q = z0/(np.linalg.norm(z0,2))
    it = 0
    err = tol + 1.
    while it < nmax and err > tol:
        z = np.dot(A,q)
        l = np.dot(q.T,z)
        err = np.linalg.norm(z-l*q,2)
        q = z/(np.linalg.norm(z,2))
        it += 1
    return l,q
    

def condNumb(A, eps, tol, nmax):
    '''
    condNumb returns the conditional number of the matrix A:
    it uses the power method to find the max eigenvalue
    and the Inverse power method with shift to find the minimum eigenvalue.
    The tolerance is tol and nmax is the max number of iterations.
    eps is a regulator for the Inverse power method.
    Actually, we look for the minimum eigenvalue close to zero (regulated by eps) 
    as our matrix is symmetric positve definite
    '''
    n = np.shape(A)[0]
    z0 = np.random.rand(n)
    lmin, _ = IPM(A, z0, eps, tol, nmax)
    lmax, _ = PM( A, z0, tol, nmax)
    condNu = lmax/lmin
    return condNu


print(condNumb(A, 1e-10, 1e-10, 10000 ))
print(np.linalg.cond(A))
```

    5278.068747699028
    5278.068747922936


Implement a preconditioned Conjugant Gradient method to solve the original linear system of equations using an iterative method:


```python
def conjugate_gradient(A, b, P, nmax, eps):
    '''
    Preconditioned Conjugant Gradient method:
    Given the matrix A, the vector b and the preconditioner matrix P
    conjugate_gradient returns the solution of the original system 
    after nmax iterations or with tolerance eps
    '''
    N=len(A)
    x = np.zeros_like(b)
    tol = eps + 1
    it = 0
    r = b - np.dot(A,x)
    rho_old = 1.
    p_old = np.zeros_like(b)
    while (it < nmax and tol > eps):
        it += 1
        z = np.linalg.solve(P,r)
        rho = np.dot(r,z)
        if (it > 1):
            beta = rho/rho_old
            p = z + beta*p_old
        else:
            p = z
        q = np.dot(A,p)
        alpha = rho/(np.dot(p,q))
        x += p*alpha
        r -= q*alpha
        p_old = p
        rho_old = rho
        tol = np.linalg.norm(r,2)    
    return x


# conjugate gradient solution
sol_conjugate_gradient = conjugate_gradient(A,b, np.identity(len(A)), np.shape(A)[0], 1e-10)
x=np.linspace(omega[0],omega[1],n)
# LU solution
L,U = LU(A)
w = L_solve(L,b)
u = U_solve(U,w)
# Numpy solution
u_ex = np.linalg.solve(A, b)
# Lets plot
_ = plot(x,sol_conjugate_gradient,'ro', label="conjugate gradient")
_ = plot(x,u_ex,'b', label="numpy")
legend()
show()
_ = plot(x,sol_conjugate_gradient,'ro', label="conjugate gradient")
_ = plot(x,u,'g', label="LU")
legend()
show()

print(np.linalg.norm(sol_conjugate_gradient-u_ex,2)/np.linalg.norm(u_ex,2))
print(np.linalg.norm(sol_conjugate_gradient-u,2)/np.linalg.norm(u,2))
```

Consider the following time dependent variation of the PDE starting from the orginal problem formulation:
$$
u'(t)-u_{xx} = \alpha(t)f(x)
$$

for $t\in [0,T]$, with $\alpha(t) = \cos(t)$ and $T = 6\pi$

Use the same finite difference scheme to derive the semi-discrete formulation and solve it using a forward Euler's method.

Plot the time dependent solution solution at $x = \pi/2$, $x=1$, 
$x=\pi$



```python
def fe_pde(omega_x, omega_t, bc, y0, h, n, f, eps ):
    '''
    Given the intervals omega_t and omega_x
    with boundary conditions y0 and bc
    h is the time distance between two consecutive nodes
    n is the number of spatial discretization points
    f is the forcing function
    eps a small parameter to calculate the timestep correctly
    fe_pde returns the solutions of the pde together with the timesteps and xsteps
    '''
    assert n > 0 
    assert h > 0
    A, b = finDif(omega_x, f, n, bc)
    timesteps = np.arange(omega_t[0],omega_t[1]+eps, h) 
    m = len(timesteps)
    sol = np.zeros((m,n))
    Au = np.zeros((m,n))
    xsteps=np.linspace(omega_x[0],omega_x[1],n)
    sol[0] = y0
    for i in range(0,m-1):
        Au[i] = np.dot(A,sol[i])   
        sol[i+1] = sol[i]  + h *(b*np.cos(timesteps[i])- Au[i] )
    return sol, timesteps, xsteps

def sol_x_fixed(sol, omega_x, x):
    assert x>= omega_x[0]
    assert x<= omega_x[1]
    n = np.shape(sol)[1]
    assert n>0
    x_index = int((n-1)*(x-omega_x[0])//(omega_x[1]-omega_x[0]))
    #if x_index<0: 
    #    x_index=0
    #if x_index>=n: 
    #    x_index=n-1    
    sol_x_fixed = sol[:,x_index]
    return sol_x_fixed

def sol_t_fixed(sol, omega_t, t):
    assert t>= omega_t[0]
    assert t<= omega_t[1]
    m = np.shape(sol)[0]
    assert m>0
    t_index = int((m-1)*(t-omega_t[0])//(omega_t[1]-omega_t[0]))
    #if t_index<0: 
    #    t_index=0
    #if t_index>=m: 
    #    t_index=m-1    
    sol_t_fixed = sol[t_index,:]
    return sol_t_fixed


### Boundary condition y0 = 0

f = lambda x : np.sin(x)
omega_x = [0,np.pi]
omega_t = [0,6*np.pi]
n=100
x=np.linspace(omega[0],omega[1],n)
z0 = np.ones_like(x)
bc = [0,0]
A, b = finDif(omega, f, n, bc)
l_max,_ = PM(A,z0,1e-15,15000)
h = 0.0001
assert h < 2/ abs(l_max)
y0 = 0
sol, tsteps, xsteps = fe_pde(omega_x, omega_t, bc, y0, h, n, f, 1e-10)

# Comparison with exact solution
t= tsteps
sol_t = sol_x_fixed(sol, omega_x, 1)
plot(t,sol_t, 'blue')
plot(t,-(1/2)*np.exp(-t)*np.sin(1)+1/2*np.cos(t)*np.sin(1)+1/2*np.sin(t)*np.sin(1),'black',linestyle='dotted', label="analytical")
legend()
show()

sol_t = sol_x_fixed(sol, omega_x,np.pi/2)
plot(t,sol_t,'green')
plot(t,-(1/2)*np.exp(-t)*np.sin(np.pi/2)+1/2*np.cos(t)*np.sin(np.pi/2)+1/2*np.sin(t)*np.sin(np.pi/2),'black',linestyle='dotted', label="analytical")
legend()
show()

sol_t = sol_x_fixed(sol, omega_x,np.pi)
plot(t,sol_t, 'red')
plot(t,-(1/2)*np.exp(-t)*np.sin(np.pi)+1/2*np.cos(t)*np.sin(np.pi)+1/2*np.sin(t)*np.sin(np.pi),'black',linestyle='dotted', label="analytical")
legend()
show()

### Boundary condition y0 = sin(x)

f = lambda x : np.sin(x)

omega_x = [0,np.pi]
omega_t = [0,6*np.pi]
n=100
x=np.linspace(omega[0],omega[1],n)
bc = [0,0]
A, b = finDif(omega, f, n, bc)
l_max,_ = PM(A,z0,1e-15,15000)
h = 0.0003
assert h < 2/ abs(l_max)
y0 = np.sin(x)

sol, tsteps, xsteps = fe_pde(omega_x, omega_t, bc, y0, h, n, f, 1e-10)

# Comparison with exact solution
t= tsteps
sol_t = sol_x_fixed(sol, omega_x, 1)
plot(t,sol_t, 'blue')
plot(t,(1/2)*np.exp(-t)*np.sin(1)+1/2*np.cos(t)*np.sin(1)+1/2*np.sin(t)*np.sin(1),'black',linestyle='dotted', label="analytical")
legend()
show()

sol_t = sol_x_fixed(sol, omega_x, np.pi/2)
plot(t,sol_t,'green')
plot(t,(1/2)*np.exp(-t)*np.sin(np.pi/2)+1/2*np.cos(t)*np.sin(np.pi/2)+1/2*np.sin(t)*np.sin(np.pi/2),'black',linestyle='dotted', label="analytical")
legend()
show()

sol_t = sol_x_fixed(sol, omega_x, np.pi)
plot(t,sol_t, 'red')
plot(t,(1/2)*np.exp(-t)*np.sin(np.pi)+1/2*np.cos(t)*np.sin(np.pi)+1/2*np.sin(t)*np.sin(np.pi),'black',linestyle='dotted', label="analytical")
legend()
show()
## Few differences between analytic solution and numerical solution but of the order of 1e-16
```

Given the original $Au = b$ system, implement an algorithm to compute the eigenvalues and eigenvectors of the matrix $A$. Exploit the computed LU factorization


```python
def lu_eigvalues(A, tol, nmax):
    """
    Given the matrix A lu_eigvalues returns the eigenvalues
    with tolerance tol or after nmax iterations 
    """
    A_new=A.copy()
    n = np.shape(A_new)[0]
    eigvals_old = np.diag(A_new)
    it = 0
    err = tol + 1.
    while(it < nmax and err > tol ):
        L, U = LU(A_new)
        A_new = np.matmul(U,L)
        eigvals_new = np.diag(A_new) 
        err = np.linalg.norm(eigvals_new-eigvals_old,2)
        it = it + 1
        eigvals_old = eigvals_new
    return eigvals_new

def eig(A, eps, tol_eigevalues, nmax_eigevalues, tol_IPM, nmax_IPM):
    """
    eig returns eigenvalues and eigenvectors of A
    eig first finds the approximated eigenvalues using LU decomposition (lueigvals)
    with tolerance tol_eigevalues or using nmax_eigevalues iterations 
    Then we use the IPM function to find the eigenvalues,
    close to the approximated eigenvalue (regulated by eps), 
    and the relative eigenvectors
    with tolerance tol_IPM or using nmax_IPM iterations
    """
    n = np.shape(A)[0]
    z0 = np.random.rand(n)
    eigvals = np.zeros(n)
    eigvecs = np.zeros((n,n))
    lueigvals = lu_eigvalues(A, tol_eigevalues, nmax_eigevalues)
    for i, qre in enumerate(lueigvals):
        eigvals[i], eigvecs[i] = IPM(A, z0, qre + eps, tol_IPM, nmax_IPM )
    return eigvals, eigvecs

l,v = eig(A, 1e-10, 1e-10, 3000, 1e-10, 10000 )
l_numpy,v_numpy = np.linalg.eig(A)
print(max(abs(np.sort(l)-np.sort(l_numpy))))
print()
print(np.sort(v)-np.sort(v_numpy))
# It takes some time...
# The algorith is not that efficient so in order to have good accuracy for the eigenvalues I need many iterations
# Still it is not enough to resolve correctly the eigenvectors 

```

    1.5006662579253316e-11
    
    [[-6.75407314e-25 -6.58353206e-25 -6.33976992e-25 ...  6.88137895e-25
       3.81118185e-01 -7.54736733e-02]
     [ 5.87127558e-04  6.27310410e-04  8.44491526e-04 ... -5.67468240e-05
      -3.55711170e-04 -6.09596892e-04]
     [ 1.32406417e-04 -1.88459072e-04  3.70312990e-04 ... -4.96275872e-05
      -1.34983203e-04 -2.40706189e-04]
     ...
     [-5.88838847e-04 -1.00888656e-03 -1.26785559e-03 ... -2.20569435e-04
      -8.34925547e-05  6.30041547e-05]
     [ 2.29808483e-04 -1.17553961e-03 -1.79866824e-03 ... -1.46720716e-01
      -1.42627391e-01 -1.42649860e-01]
     [-6.75407314e-25 -6.58353206e-25 -6.33976992e-25 ...  6.88137895e-25
       3.81118185e-01 -7.54736733e-02]]


Compute the inverse of the matrix A exploiting the derived LU factorization


```python
def inverse(A):
    """
    inverse returns the inverse of the matrix A
    """
    m = np.shape(A)[0]
    Ainv= np.ones_like(A)
    L,U = LU(A)
    for j in range(m):
        e=np.zeros((m,))
        e[j]=1
        winv = L_solve(L,e)
        uinv = U_solve(U,winv)
        Ainv[:,j]=uinv
    return Ainv    

Ainv = inverse(A)
Ainv_numpy = np.linalg.inv(A)
print(np.linalg.norm(Ainv-Ainv_numpy,2))
```

    6.492862786734044e-14


Consider the following Cauchy problem
$$
\begin{cases}
y'= -ty^2 \quad 0\le t \le 2\\
y(0) = 1
\end{cases}
$$
Implement a Backward Euler's method in a suitable function and solve the resulting non-linear equation using a Newton's method.


```python
def newton(f, f_prime, x0, eps, n_max):
    '''
    Newton's method:
    f: function
    f_prime: derivative of the function
    x0: initial point
    eps: tolerance
    n_max: max number of iterations
    '''
    err = abs(f(x0))
    it = 0
    x = x0
    while (err > eps and it < n_max):
        if abs(f_prime(x)) < 1e-12:
            raise RuntimeError("f_prime(x) is close to zero")
        x_new = x - f(x)/f_prime(x)
        err = abs(f(x_new))
        x = x_new
        it += 1
    return x

def fe(t0, tf, h, y0, eps=1e-10):
    '''
    Forward Euler's method:
    given
    t0:  initial time
    tf:  final time
    h:   time distance between two consecutive nodes
    y0:  boundary condition at y[t0]
    eps: small epsilon to calculate the time steps correctly
    fe returns the solution of the ODE and the timesteps
    '''
    timesteps = np.arange(t0,tf+eps, h)
    sol = np.zeros_like(timesteps)
    sol[0] = y0
    for i in range(1,len(sol)):
        sol[i] = sol[i-1] - h * timesteps[i-1] * sol[i-1]**2   
    return sol, timesteps


def be(t0, tf, h, y0, tol=1e-10, nmax=1000, eps=1e-10):
    """
    Backward Euler's method:
    given
    t0: initial time
    tf: final time
    h: the time distance between two consecutive nodes
    y0: boundary condition at y[t0]
    tol: tolerance for the iterations
    nmax: max number of iterations
    eps: small epsilon to calculate the time steps correctly
    be returns the solution of the ODE and the timesteps
    """
    timesteps = np.arange(t0, tf + eps, h)
    sol = np.zeros_like(timesteps)
    sol[0] = y0
    it = 0
    for i in range(len(sol)-1):
        t = sym.symbols('t')
        f_sym = t - sol[i] + h * timesteps[i+1] * t**2 
        f_prime_sym = sym.diff(f_sym,t)
        f = sym.lambdify(t, f_sym, 'numpy')
        f_prime = sym.lambdify(t, f_prime_sym, 'numpy')
        sol[i+1] = newton(f, f_prime, sol[i], tol, nmax) 
    return sol, timesteps

l = -5.
t0 = 0.
tf = 2.
y0 = 1.
s = np.linspace(t0,tf,100)
h=0.1
#h=0.01
f_exact = lambda x: 2/(2 + x**2)
f_fe, t = fe(t0, tf, h, y0)
f_be, t = be(t0, tf, h, y0)
_ = plot(t, f_fe, 'r', label="forward euler")
_ = plot(t, f_be, 'g', label="backward euler")
_ = plot(t, f_exact(t), 'b', label="exact")
legend()
error1 = max(abs(f_exact(t) - f_fe))
error2 = max(abs(f_exact(t) - f_be))
print(error1, error2)
```


```python

```

ed804951b94c8e069cd9b8981093c04f6714577c
ipynb
Jupyter Notebook
FinalProject_2021_2022.ipynb
marco-celoria/numerical-analysis-2021-2022
28e455aff1ef603276f64be3627a94e8035c2680
CC-BY-4.0
FinalProject_2021_2022.ipynb
marco-celoria/numerical-analysis-2021-2022
28e455aff1ef603276f64be3627a94e8035c2680
CC-BY-4.0
FinalProject_2021_2022.ipynb
marco-celoria/numerical-analysis-2021-2022
28e455aff1ef603276f64be3627a94e8035c2680
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES