```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve_triangular
import warnings
from functools import partial
```

## Finite differences

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,second_order=False):
    """
    Compute the Finite Difference matrix for the differential 
    equation :math:`-u'' = f` in the given domain `omega`, using
    a grid with `n` elements and the given boundary conditions
    `bc`.
    
    :return: The matrix for the FD problem and the RHS vector.
    """
    non_null_component = np.array([-1, 16, -30, 16, -1])
    f = np.vectorize(f)
    
    h = (omega[1] - omega[0]) / n
    
    A = np.zeros((n,n))
    rhs = f(np.linspace(*omega, n))
    
    if not second_order:
        # fourth order
        for i in range(2,n-2):
            A[i,i-2:i+3] = np.array(non_null_component) / 12
        # second order if the fourth order is not applicable
        A[1,0:3] = A[-2,-3:] = np.array([1,-2,1])
    else:
        for i in range(1,n-1):
            A[i,i-1:i+2] = np.array([1,-2,1])
    
    A /= h*h
    
    # boundary conditions
    A[0,0] = A[-1,-1] = 1
    rhs[0] = bc[0]
    rhs[-1] = bc[1]
    
    return A, rhs
```

Call the function using:


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

## LU and Cholesky factorization 

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


```python
# row-partial pivoting
def gauss_factorization(A, row_pivoting=True):
    """
    Operate Gauss factorization on the matrix `A` (with optional
    row-pivoting depending on the value of the parameter 
    `row_pivoting`).
    
    :return: The upper triangular matrix obtained at the end of 
                the algorithm, the lower triangular Gauss coefficient
                matrix, and (if `row_pivoting` is `True`) the 
                permutation matrix.
    """ 
    U = np.array(A)
    gauss_coeff_matrix = np.zeros_like(U)
    exchanged_rows = []
    
    for i in range(U.shape[0]-1):
        if row_pivoting:
            # we look for the biggest element in the column i. we need 
            # to add i to row_max since the index given by argmax is 
            # relative to the submatrix A[i:,i].
            row_max = i + np.argmax(np.abs(U[i:,i])) 
            exchanged_rows.append((i, row_max))
            U[[i,row_max]] = U[[row_max,i]]
            gauss_coeff_matrix[[i,row_max]] = gauss_coeff_matrix[[row_max,i]]
        
        gauss_coeff_matrix[i+1:,i] = U[i+1:,i] / U[i,i]
        U[i+1:] -= np.multiply(gauss_coeff_matrix[i+1:,i][:,None], np.repeat(U[i][None], U.shape[0] - (i+1), axis=0))
    
    P = np.identity(A.shape[0])
    for a,b in exchanged_rows[::-1]:
        Pi = np.identity(U.shape[0])
        Pi[[a,b]] = Pi[[b,a]]
        P = P.dot(Pi)
            
    if row_pivoting:
        return U, gauss_coeff_matrix, P
    else:
        return U, gauss_coeff_matrix

def LU(A, row_pivoting=False):
    tp = gauss_factorization(A, row_pivoting)
    return (tp[1] + np.diag(np.ones(A.shape[0])), tp[0], *tp[2:])

# test gaussian elimination
temp_matrix = np.array([[4,2,3],[2,4,5],[7,8,9]])
np.testing.assert_almost_equal(np.dot(*LU(temp_matrix)), temp_matrix)

# test gaussian elimination with pivoting
L,U,P = LU(temp_matrix, row_pivoting=True)
np.testing.assert_almost_equal(L.dot(U), P.dot(temp_matrix))
```


```python
def cholesky(A):
    """
    Attemp Cholesky factorization on the matrix `A`. If the output
    contains `np.nan`, the factorization failed. No checks are
    performed.
    """
    warnings.filterwarnings('ignore')
    
    H = np.zeros_like(A)
    H[0,0] = np.sqrt(A[0,0])
    for j in range(1, A.shape[0]):
        for i in range(j):
            H[i,j] = (A[i,j] - H[:i,i].dot(H[:i,j])) / H[i,i]
        H[j,j] = np.sqrt(A[j,j] - np.einsum('i,i', H[:j,j], H[:j,j]))
        
    warnings.resetwarnings()
    return H.T, H

# test cholesky factorization on a positive definite matrix
temp_matrix = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]], dtype=float)
np.testing.assert_almost_equal(np.dot(*cholesky(temp_matrix)), temp_matrix)
```


```python
def triangular_factorization(A):
    """
    Attempt a triangular factorization of the given matrix using
    Cholesky factorization, fallbacks to LU factorization with
    Gauss elimination method.
    """
    L,U = cholesky(A)
    if np.count_nonzero(np.isnan(L)) > 0:
        L,U = LU(A)
    return L,U

def check_triangular(T, tolerance=1.e-14, verbose=False):
    """
    Check that the given matrix is triangular (elements whose module is
    smaller than `tolerance` are assumed to be zero).
    
    ..warning
    
        Diagonal matrices are recognized as upper triangular matrices.
    
    :return: `'upper'` if `T` is upper triangular, `'lower'` if it is
                lower triangular, `None` otherwise.
    """
    if T.shape[0] != T.shape[1]:
        return False
    if T.shape[0] == 1:
        return True
    
    T = np.abs(T)
    
    upper = lower = True
    for i in range(T.shape[0]):
        uvec = T[i+1:,i]
        if verbose:
            print(i,uvec)
        if upper and np.count_nonzero(uvec > tolerance) != 0:
            upper = False
            if verbose:
                print('non upper (i={}): {}'.format(i,uvec))
            
        lvec = T[i,i+1:]
        if lower and np.count_nonzero(lvec > tolerance) != 0:
            lower = False
            if verbose:
                print('non lower (i={}): {}'.format(i,lvec))
    
    if upper:
        return 'upper'
    elif lower:
        return 'lower'
    return None
    
# test factorization on a non-(positive definite) matrix
temp_matrix = np.array([[-2,1,0],[-1,2,-1],[0,-1,2]], dtype=float)
np.testing.assert_almost_equal(np.dot(*triangular_factorization(temp_matrix)), temp_matrix)

# test check triangular
assert check_triangular(np.array([[1,1,1],[0,1,1],[0,0,0]])) == 'upper'
assert check_triangular(np.array([[1,0,0],[-1,1,0],[1,1,0]])) == 'lower'
assert check_triangular(np.array([[1,1,1],[0,1,1],[1,0,0]])) == None
```

## Direct methods for linear systems

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


```python
def solve_triangular_system(M,rhs,lower):
    """
    Substitution method for triangular matrices (upper or lower depending
    on the parameter `lower`).
    """
    x = np.zeros_like(M[0])
    
    rng = range(M.shape[0]) if lower else range(M.shape[0]-1,-1,-1)
    for i in rng:
        if lower:
            mi = M[i,:i]
            xi = x[:i]
        else:
            mi = M[i,i+1:]
            xi = x[i+1:]
        x[i] = (rhs[i] - np.dot(xi, mi)) / M[i,i]
        
    np.testing.assert_almost_equal(x, solve_triangular(M,rhs,lower=lower))
    return x
```


```python
L_solve = partial(solve_triangular_system, lower=True)
U_solve = partial(solve_triangular_system, lower=False)
```


```python
def solve_direct(A,rhs,*,assertions=True,L=None,U=None,P=None):
    """
    Solve the given linear system using a direct 
    method (i.e. LU or Cholesky).
    """
    if L is None or U is None:
        L,U = triangular_factorization(A)
    if P is None:
        P = np.identity(len(A))
    
    if assertions:
        np.testing.assert_almost_equal(P.dot(A), L.dot(U))
        assert check_triangular(L, tolerance=5.e-10) == 'lower'
        assert check_triangular(U, tolerance=5.e-10) == 'upper'
    
    y = L_solve(L,P.dot(rhs))
    if assertions:
        np.testing.assert_almost_equal(L.dot(y), P.dot(rhs))

    x = U_solve(U,y)
    if assertions:
        np.testing.assert_almost_equal(U.dot(x), y)
        
    return x
```

## Solution of the FD problem

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


```python
def solve_fin_diff_problem(omega, f, n, bc, solver=solve_direct, second_order=False):
    """
    Return the function `x` which satisfies the ODE 
    :math:`-x'' = f`, and the discrete grid used to 
    find the solution.
    
    By default the function uses a direct solver, but different
    solvers may be used by changing the parameter `solver`.
    
    :return: A tuple whose first item is the discrete grid used 
                to compute the solution, the second item is the
                solution computed on the grid.
    """
    A, rhs = finDif(omega, f, n, bc, second_order=second_order)
    x = -solver(A, rhs)
    np.testing.assert_almost_equal(-A.dot(x), rhs)
    return np.linspace(*omega, n), x
```


```python
plt.xlabel('$\\Omega$')
plt.xticks(omega, ['0','$\\pi$'])
plt.plot(*solve_fin_diff_problem(omega, f, n, bc), label='IV order')
plt.plot(*solve_fin_diff_problem(omega, f, n, bc, second_order=True), label='II order')
plt.legend()
plt.show()
```

## Error plot

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
omega2 = [0,1]
f2 = lambda x : x*(1-x)
bc2 = [0,0]
ns = np.arange(10,1011,100)

exact = lambda x: np.power(x,4) / 12 - np.power(x,3) / 6 + x / 12
```


```python
def compute_error(n, error_function):
    xs, x = solve_fin_diff_problem(omega2, f2, n, bc2)
    return error_function(np.abs(exact(xs) - x))

error_infty_norm = partial(compute_error, error_function=np.max)
error_mean = partial(compute_error, error_function=np.mean)
err_funcs = [error_infty_norm, error_mean]

labels = ['$\\Vert \\cdot \\Vert_\\infty$', 'mean']
    
plt.title('Errors')
for err_func, label in zip(err_funcs, labels):
    err = list(map(lambda n: err_func(n), ns))
    plt.plot(ns, err, label=label)
plt.legend()

plt.xlabel('n')
plt.ylabel('Error')

plt.yscale('log')

plt.show()
```

## Inverse

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


```python
def inverse(A):
    """
    Compute the inverse matrix of `A` via triangular factorization.
    """
    L,U = triangular_factorization(A)
    A_inv = np.zeros_like(A)
    identity = np.identity(A.shape[0])
    
    for j in range(A.shape[0]):
        A_inv[:,j] = solve_direct(A, identity[:,j])
    return A_inv

# a quick test
plt.pcolormesh(inverse(A).dot(A));
```

## Condition number

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


```python
def infty_norm(A):
    """
    Returns :math:`\\Vert A \\Vert_\\infty`.
    """
    return np.max(np.einsum('ij->i', np.abs(A)))

# infty norm
def condNumb(A):
    """
    Compute the condition number of `A` in the 
    :math:`\\Vert \\cdot \\Vert_\\infty` norm."""
    return infty_norm(A) * infty_norm(inverse(A))

condNumb(A)
```




    11937.129794259363



## Preconditioned conjugate gradient method

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=None, eps=1e-10):
    """
    Preconditioned gradient method with tolerance `eps` and
    preconditioner `P`. If `P` is `None`, the preconditioner
    is the diagonal matrix of `A`.
    """
    if P is None:
        P = np.identity(A.shape[0])
    
    x = np.zeros_like(A[0])
    r = b - A.dot(x)
    z = np.linalg.solve(P, r)
    p = z
    
    for i in range(len(A)):
        if np.linalg.norm(r) < eps:
            break
            
        a = p.dot(r) / A.dot(p).dot(p)
        x += a*p
        r -= a*A.dot(p)
        z = np.linalg.solve(P, r)
        beta = A.dot(p).dot(z) / A.dot(p).dot(p)
        p = z - beta*p
        
    if np.linalg.norm(A.dot(x) - b) > eps:
        warnings.warn("The given matrix may not be symmetric "
                      "positive definite, PCG did not converge.")
        
    return x
```


```python
xs, sol_conj = solve_fin_diff_problem(omega, f, 100, bc, solver=conjugate_gradient, second_order=True)
_, sol_direct = solve_fin_diff_problem(omega, f, 100, bc)

np.linalg.norm(sol_conj - sol_direct) / np.linalg.norm(sol_direct)
```




    8.390564146200464e-05




```python
plt.figure(figsize=(20,4))

plt.subplot(1,3,1)
plt.xlabel('$\\Omega$')
plt.plot(xs, sol_direct)
plt.xticks(omega, ['0','$\\pi$'])
plt.title('Direct')

plt.subplot(1,3,2)
plt.xlabel('$\\Omega$')
plt.plot(xs, sol_conj)
plt.xticks(omega, ['0','$\\pi$'])
plt.title('PCG')

plt.subplot(1,3,3)
plt.xlabel('$\\Omega$')
plt.plot(xs, np.abs(sol_direct - sol_conj))
plt.xticks(omega, ['0','$\\pi$'])
plt.title('Difference between PCG and direct solution')

plt.show()
```

## Forward euler

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$

### Forward Euler method
We discretize the time using the Forward Euler approximation for the first derivative with respect to $t$:
$$u_{xx} + \alpha(t)f(x) = \frac{\partial u}{\partial t}(x,t) \approx \frac{\partial u}{\partial t}(x,t_n) \approx \frac{u(x,t_{n+1}) - u(x,t_n)}{h_{\text{t}}}$$
$$\implies u(x,t_{n+1}) = u(x,t_n) + h_{\text{t}} u_{xx}(x,t_n) + h_{\text{t}} \alpha(t_n)f(x)$$

We discretize also the space:
$$u_{xx}(x,t) \approx u_{xx}(x_i,t) = \frac{-u(x_{i-2},t)+16u(x_{i-1},t)-30u(x_i,t)+16u(x_{i+1},t)-u(x_{i+2},t)}{12h_x^2}$$

And we obtain the following numerical scheme:
$$\begin{align} u(x_m,t_{n+1}) &= u(x_m,t_n) + h_{\text{t}} \frac{-u(x_{m-2},t_n)+16u(x_{m-1},t_n)-30u(x_m,t_n)+16u(x_{m+1},t_n)-u(x_{m+2},t_n)}{12h_x^2} + h_{\text{t}} \alpha(t_n)f(x_m)\\
&= h_{\text{t}} \frac{-u(x_{m-2},t_n)+16u(x_{m-1},t_n)-(30 - 12h_x^2/h_t)u(x_m,t_n)+16u(x_{m+1},t_n)-u(x_{m+2},t_n)}{12h_x^2} + h_{\text{t}} \alpha(t_n)f(x_m)\\
&= h_{\text{t}} \left(\frac{-u(x_{m-2},t_n)+16u(x_{m-1},t_n)-(30 - 12h_x^2/h_t)u(x_m,t_n)+16u(x_{m+1},t_n)-u(x_{m+2},t_n)}{12h_x^2} + \alpha(t_n)f(x_m)\right)
\end{align}$$

We can now write the problem in the following form for some matrix $A \in \mathbb{R}^{N_x \times N_x}$:
$$\begin{bmatrix} u(x_0)\\ \vdots \\ u(x_{N_x})\end{bmatrix}(t_{n+1}) = A \begin{bmatrix} u(x_0)\\ \vdots \\ u(x_{N_x})\end{bmatrix}(t_n) + \begin{bmatrix} f(x_0)\\ \vdots \\ f(x_{N_x})\end{bmatrix} h_t \alpha(t_n)$$
where $N_x$ is an integer such that $N_x h_x$ is equal to the right boundary of the space domain.


```python
def forward_euler_A(ht, hx, Nx):
    A = np.zeros((Nx,Nx), dtype=float)
    
    non_null_component = np.array([-1, 16, -30+12*hx**2/ht, 16, -1])

    # fourth order
    for i in range(2,Nx-2):
        A[i,i-2:i+3] = np.array(non_null_component) / 12
    # second order if the fourth order is not applicable
    A[1,0:3] = A[-2,-3:] = np.array([1,-2+hx**2/ht,1])
    
    A /= hx**2
    A *= ht
    
    return A

def forward_euler(u0, f, alpha, omega_x, omega_t, Nx, Nt):
    if not isinstance(u0, np.ndarray):
        u0 = np.repeat(u0, Nx)
    
    ht = (omega_t[1] - omega_t[0]) / Nt
    hx = (omega_x[1] - omega_x[0]) / Nx
    
    xs = np.linspace(*omega_x, Nx)
    ts = np.linspace(*omega_t, Nt)
    
    A = forward_euler_A(ht, hx, Nx)
    f = np.vectorize(f)(xs)
    
    u = np.zeros((Nx, Nt), dtype=float)
    u[:,0] = u0
    
    # boundary conditions
    u[[0,-1],0] = 0
    f[[0,-1]] = 0
    
    for i in range(1,Nt):
        u[:,i] = A.dot(u[:,i-1]) + ht * alpha(i*ht) * f
        
    return xs,ts,u
```


```python
def x_idx(xs, x):
    """
    Find the best value in xs to approximate the scalar x.
    """
    return np.argmin(np.abs(xs - x))

def show_solution(Nx, Nt, *, xs=None):
    xs, ts, u = forward_euler(0, np.sin, np.cos, (0,np.pi), (0,6*np.pi), Nx=Nx, Nt=Nt)

    pi_2_idx = x_idx(xs, np.pi/2)
    one_idx = x_idx(xs, 1)
    pi_idx = -1
    
    my_pi2 = np.round(xs[pi_2_idx],4)
    my_one = np.round(xs[one_idx],4)
    my_pi = np.round(xs[pi_idx], 4)

    time_slice = slice(None,None,100)

    plt.plot(ts[time_slice], u[pi_2_idx,time_slice], '-', label='x={}'.format(my_pi2))
    plt.plot(ts[time_slice], u[one_idx,time_slice], '-.', label='x={}'.format(my_one))
    plt.plot(ts[time_slice], u[pi_idx,time_slice], '--', label='x={}'.format(my_pi))

    plt.xticks([i*np.pi for i in range(7)], ['0', *('{}$\\pi$'.format(i) for i in range(1,7))])

    plt.xlabel('t')
    plt.ylabel('u(x,t)')

    plt.title('Time dependent solution (Nx={}, Nt={})'.format(Nx,Nt))
    plt.legend()
    
plt.figure(figsize=(20,4))

plt.subplot(1,3,1)
show_solution(50,10000)
plt.subplot(1,3,2)
show_solution(50,100000)
plt.subplot(1,3,3)
show_solution(100,100000)
plt.show()
```

## Eigenvalues and eigenvectors

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 power_method(A, max_iter=1000, max_residual=1.e-10, inverse=False, verbose=False):
    x = np.ones_like(A[0])
    
    if inverse:
        L,U,P = LU(A, row_pivoting=True)
        if np.count_nonzero(np.isnan(L)) > 0 or np.count_nonzero(np.isnan(U)) > 0:
            return None,None
        np.testing.assert_almost_equal(L.dot(U), P.dot(A))
    
    for i in range(max_iter):
        y = x / np.linalg.norm(x)
        
        # update x
        if inverse:
            x = solve_direct(A, y, assertions=False, L=L, U=U, P=P)
        else:
            x = np.dot(A,y)
        
        # rayleigh equation
        eig = complex(y.dot(x))
        
        if np.linalg.norm(A.dot(y) - y * eig) < max_residual:
            break
    if verbose:
        print('Stopped after {} iterations'.format(i+1))
        
    return 1/eig if inverse else eig, y

def compute_some_eigenvalues(A, equal_tolerance=1.e-10):
    eigs = []
    eigenvalues = []
    
    # we look for the eigenvalues near the center of the gershgorin circles (i.e.
    # diagonal elements of A).
    for diag_element in [*np.diag(A)]:
        eig, vec = power_method(A - diag_element * np.identity(len(A)), inverse=True)
        if eig != None:
            eig += eig + diag_element
            eigs.append(eig)
            eigenvalues.append(vec)
            
    # add the smallest eigenvalue of A
    eig, vec = power_method(A, inverse=True)
    if eig != None:
        eigs.append(eig)
        eigenvalues.append(vec)
    # and the biggest
    eig, vec = power_method(A)
    eigs.append(eig)
    eigenvalues.append(vec)
        
    srt_idx = np.argsort(np.abs(eigs))
    eigs = np.array(eigs)[srt_idx]
    eigenvalues = np.array(eigenvalues)[srt_idx]
    
    # indexes of unique eigenvalues, since we may find the same eigenvalue
    # multiple times
    non_repeated_eigs_idxes = [0]
    for i in range(1,len(eigs)):
        if np.abs(eigs[non_repeated_eigs_idxes[-1]] - eigs[i]) > equal_tolerance:
            non_repeated_eigs_idxes.append(i)
            
    eigs = eigs[non_repeated_eigs_idxes]
    eigenvalues = eigenvalues[non_repeated_eigs_idxes]
    
    if len(eigs) > A.shape[0]:
        warnings.warn("The number of unique eigenvalues found is greater than the dimension "
                      "of the matrix, there are probably some spurious eigenvalues that "
                      "should be discarded.")
    return eigs, eigenvalues
```


```python
A, _ = finDif(omega, f, n, bc)
eigs, eigenvecs = compute_some_eigenvalues(A)
```

    /u/f/fandreuz/.local/lib/python3.6/site-packages/ipykernel_launcher.py:27: RuntimeWarning: invalid value encountered in true_divide



```python
true_eigs = np.linalg.eig(A)[0]
for eig, eigenvec in zip(eigs,eigenvecs):
    idx = np.argmin(np.abs(true_eigs - eig))
    print('Found {}, relative error={}, residual={}'.format(eig, 
                                                            np.abs(true_eigs - eig)[idx] / np.abs(true_eigs[idx]),
                                                            np.linalg.norm(A.dot(eigenvec) - eig * eigenvec)))
```

    Found (0.9999999980238907+0j), relative error=1.9761092762138333e-09, residual=4.412018730890187e-09
    Found (-2083.933510491916+0j), relative error=0.01399144528808244, residual=28.75491882258004
    Found (-2589.5416917742764+0j), relative error=0.011031979370446372, residual=28.25605035791613
    Found (-5396.812957782051+0j), relative error=6.997473027363431e-06, residual=0.8857371160328106



```python
# a more refined approach may use this function to locate the area
# of the complex plane where the eigenvalues are located
def gershgorin_radius(A, row=True):
    # we ignore diagonal entries
    A = np.array(A) - np.diag(np.diag(A))
    return np.einsum('ij->{}'.format('i' if row else 'j'), np.abs(A))
```

## Cauchy problem with Backward Euler method

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.

### Backward Euler method:
$$-ty_n^2 = y' \approx \frac{y_{n} - y_{n-1}}{h}$$
$$\implies thy_n^2 + y_n - y_{n-1} =: \phi(y_n) = 0$$
$$\phi'(x) = 2th + 1$$
### Newton:
$$\begin{align}x_{n+1} &= x_n - \frac{\phi'(x_{n})}{\phi(x_{n})}\\
&= x_n - \frac{2th + 1}{thx_n^2 + x_n - y_{n-1}}
\end{align}$$


```python
def newton(y, y_prev, n, h, phi, phi_prime, iterations=0, nmax=100, tolerance=1.e-12):
    """
    Apply Netwon method to the problem above to compute the n-th 
    value of the solution y in the discretized time grid 
    (i.e. y_n).
    
    y_prev is the value of the function computed in the time instant
    n-1 (i.e. y_{n-1}).
    """
    
    t = n*h
    y -= phi(y,y_prev,t) / phi_prime(t)
    
    if phi(y,y_prev,t) < tolerance or iterations >= nmax:
        return y
    else:
        return newton(y, y_prev, n, h, phi, phi_prime, iterations+1, nmax, tolerance)

def backward_euler(omega,n,bc):
    h = (omega[1] - omega[0]) / n
    
    def phi(x, y, t):
        return t*h*x**2 + x - y

    def phi_prime(t):
        return 2*t*h + 1
    
    y = np.zeros(n, dtype=float)
    y[0] = bc
    for i in range(1,n):
        y[i] = newton(y[i-1], y[i-1], i, h, phi, phi_prime)
    return y

def exact(x):
    return 1/(x**2/2+1)

n = 50
omega3 = (0,2)
plt.plot(np.linspace(*omega3,n), backward_euler(omega3, n, 1), '.')
plt.plot(np.linspace(*omega3,n), exact(np.linspace(*omega3,n)));
```

7258a89a13ad9ab4eca231da30b93a6eacca91df
ipynb
Jupyter Notebook
FinalProject_2021_2022.ipynb
fAndreuzzi/numerical-analysis-2021-2022
49b7366727a7b4343fabafe7252fd0828dcb644b
CC-BY-4.0
FinalProject_2021_2022.ipynb
fAndreuzzi/numerical-analysis-2021-2022
49b7366727a7b4343fabafe7252fd0828dcb644b
CC-BY-4.0
FinalProject_2021_2022.ipynb
fAndreuzzi/numerical-analysis-2021-2022
49b7366727a7b4343fabafe7252fd0828dcb644b
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES