{ "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 }