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