{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1. Exercise on Interpolation\n",
"\n",
"We aim to implement a numerical method to approximately solve Laplace's equation with homogeneous Dirichlet boundary conditions:\n",
"\n",
"$$ u\\in C^2([0,1]),\\quad\\left\\{\\begin{aligned} -u''(x) & = \\varphi(x) & \\forall\\, x\\in(0,1),\\\\ u(0) & = u(1) = 0. \\end{aligned}\\right.$$\n",
"\n",
"To do this, we approximate $\\varphi$ with an interpolating polynomial $\\widehat \\varphi$,\n",
"and then we exactly solve Laplace's equation with the right-hand side $\\widehat \\varphi$ instead of $\\varphi$. Write a function `approx(n)` implementing this approach.\n",
" The function should return a polynomial approximation of the solution based on an interpolation of **degree** $n$ of the right-hand side at equidistant points between 0 and 1, inclusive.\n",
" We will take the right-hand side to be the function\n",
" $$\\varphi(x) = \\exp\\Bigl(\\sin(2\\pi x)\\Bigr) \\Bigl(\\sin(2\\pi x)-\\cos(2\\pi x)^2 \\Bigr),$$\n",
" in which case the analytical solution is given by\n",
" $$ u(x)=(2\\pi)^{-2}\\Bigl(\\exp\\bigl(\\sin(2\\pi x)\\bigr)-1\\Bigr).$$\n",
"\n",
"*Hints:*\n",
"- You can use the `fit` function from the `Polynomials.jl` library to obtain the interpolating polynomial:\n",
"\n",
" ```julia\n",
" p = fit(x, y)\n",
" ```\n",
"\n",
" where `x` are the interpolation nodes, and `y` are the values of the function to interpolate.\n",
" \n",
"- To calculate the analytical solution with a polynomial right-hand side, you can notice that all solutions are polynomials, and without boundary conditions, the solution is unique modulo $\\mathbf{P}_1$.\n",
" \n",
"- You can use the `integrate` function from the `Polynomials.jl` library, which calculates an antiderivative of a polynomial:\n",
"\n",
" ```julia\n",
" P = integrate(p)\n",
" ```\n",
"\n",
"- Use the `BigFloat` format to limit rounding errors."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. Numerical Integration\n",
"\n",
"1. Create a function that takes as input a function `u`, bounds `a` and `b`, and a number `N` of integration nodes and returns the numerical integration by Gauss-Legendre of the function `u` over the interval `[a, b]`.\n",
"\n",
"1. Create a function that takes as input a function `u`, bounds `a` and `b`, and a step size `h`, and returns the integration of `u` between `a` and `b` using the composite trapezoidal method with step size `h`. Start by converting the step size `h` into the number of intervals that best approximates `h` with the following formula `N = round(Int,(b-a)/h)`.\n",
"\n",
"1. Do the same with Simpson's method.\n",
"\n",
"1. Integrate the Runge function $u(x)=1/(1+25x^2)$ on $[-1,1]$ analytically, then plot the relative errors for the different numerical integration methods as a function of the number of integration intervals. You can use a logarithmic progression of the number of intervals between $2^3$ and $2^8$ (create an array with `2 .^(3:8)`) and plot the errors on logarithmic scales for both x and y axes (use options `xaxis=:log10` and `yaxis=:log10` in a `plot`).\n",
"\n",
"1. Apply the Romberg method to 1 and 2 iterations using the trapezoidal integration formula. In other words, implement $$J_1(h)=\\frac{4J(h)-J(2h)}{3}, \\qquad J_2(h)=\\frac{16J_1(h)-J_1(2h)}{15}$$\n",
"if $J(h)$ represents the result of the trapezoidal integration, and add the error graphs obtained to the previous ones.\n",
"\n",
"*Hints:*\n",
"\n",
"- If `U` and `V` are column vectors, `[U ; V]` returns the column vector obtained by concatenating `U` and `V`.\n",
"\n",
"- `sum([U ; V])` returns the sum of all elements in `U` and `V`, which is equivalent to `sum(U) + sum(V)`.\n",
"\n",
"- `U[1]` or `U[begin]` returns the first element of vector `U`.\n",
"\n",
"- `U[end]` returns the last element of vector `U`.\n",
"\n",
"- `U[end-1]` returns the second-to-last element of vector `U`, and so on.\n",
"\n",
"- `U[2:2:end-1]` returns the vector extracted from `U` containing only even indices between the second and second-to-last elements.\n",
"\n",
"- `U[3:2:end-2]` returns the vector extracted from `U` containing only odd indices between the third and third-to-last elements."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. Exercise on Cholesky Factorization\n",
"\n",
"The goal of this exercise is to develop an algorithm for performing the Cholesky factorization of a real positive definite matrix $\\mathsf{A}\\in\\mathbb{R}^{n×n}$, i.e., to find the real lower triangular matrix $\\mathsf{C}$ such that:\n",
"\n",
"$$\n",
"\\tag{1}\n",
"\\mathsf{A}=\\mathsf{C}\\mathsf{C}^T\n",
"$$\n",
"\n",
"1. Create a function that takes a matrix `A` as an argument and returns the matrix `C`.\n",
" To compute $\\mathsf{C}$, use a sequential identification of the columns of `C`, starting with the diagonal element and then the rest of the column, by comparing both sides of (1).\n",
"\n",
"2. Create a function to generate a random positive definite matrix in $\\mathbb{R}^{n×n}$. You can generate a random matrix `B` using `rand(n, n)` and then calculate `B'B+I`, where the identity `I` is the `UniformScaling` object from the `LinearAlgebra` library. (*Adding the identity helps avoid potential eigenvalues too close to `0`.*)\n",
"\n",
"3. Test the above algorithms with code like the following.\n",
" ```julia\n",
" n = 1000\n",
" A = generate_defpos_matrix(n)\n",
" C = cholesky(A)\n",
" norm(C*C' - A, Inf)\n",
" ```\n",
"\n",
"4. Now, suppose the real positive definite matrix `A` has a bandwidth `b` that is much smaller than `n`. Rewrite the Cholesky decomposition function to take advantage of the bandwidth.\n",
"\n",
"5. Create a function to generate a random positive definite matrix with a given bandwidth `b`. You can start by generating a random lower triangular matrix `B` with bandwidth `b`, and then return `B'B+I`.\n",
"\n",
"6. Test the algorithms from the previous two questions. For example, use `n=1000` and `b=4`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function cholesky(A)\n",
" # Your code comes here\n",
" return C\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function generate_defpos_matrix(n)\n",
" B = rand(n,n)\n",
" return B'B+I\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n = 1000\n",
"A = generate_defpos_matrix(n)\n",
"C = cholesky(A)\n",
"norm(C*C' - A, Inf)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function cholesky_banded(A, b)\n",
" # Your code comes here\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function generate_banded_matrix(n, b)\n",
" C = [j≤i≤j+b ? rand() : zero(Float64) for i in 1:n, j in 1:n]\n",
" return C*C'+I\n",
"end\n",
"\n",
"n = 10 ; b = 4\n",
"A = generate_banded_matrix(n,b)\n",
"C = cholesky_banded(A,b)\n",
"norm(C*C' - A, Inf)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.9.3",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}