psets/pset6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "01ad88a2-a4dc-457b-90d5-c7e14eff71df",
   "metadata": {},
   "source": [
    "# 18.S190/6.S090 Problem Set 6\n",
    "\n",
    "Due Friday 3/21 at **11:00am**; note that this is **earlier** than the usual midnight due date because of MIT regulations.  Email us if you need a *brief* extension.\n",
    "\n",
    "20% penalty if it is turned by **11:59pm 3/22**, and after that late psets will not be accepted.   Submit in PDF format: a decent-quality scan/image of any handwritten solutions (e.g. get a scanner app on your phone or use a tablet), combined with a PDF printout of your Jupyter notebook showing your code and (clearly labeled) results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32a513d9-1e1b-41dc-8e06-83555b60926b",
   "metadata": {},
   "source": [
    "## Problem 1 (2+5+5+5 points)\n",
    "\n",
    "Simpson’s method for integration of ODEs $\\frac{dv}{dt} = f(v, t)$ is given by:\n",
    "\n",
    "$$\n",
    "\\frac{v^{n+1} - v^{n-1}}{2\\Delta t} = \\frac{1}{6} \\left[ f(v^{n+1}, t^{n+1}) + 4f(v^{n}, t^{n}) + f(v^{n-1}, t^{n-1}) \\right].\n",
    "$$\n",
    "\n",
    "1. Is the method explicit or implicit?\n",
    "2. Using a Taylor series analysis, calculate the local truncation error for the method. Specifically, what is the order of accuracy?\n",
    "3. Is the method convergent? Why or why not?\n",
    "4. Plot the linear-stability region ($f(v,t) = \\lambda v$) for this integration method as a function of $\\lambda \\Delta t$ in the complex plane. For purely real, negative $\\lambda < 0$, what is the largest timestep $\\Delta t$ that can be taken while remaining stable for this method?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab979e06-22c6-4832-8ba9-f70bf5b5142b",
   "metadata": {},
   "source": [
    "## Problem 2 (5+5+5 points)\n",
    "\n",
    "\n",
    "Consider the integral\n",
    "$$\n",
    "I = \\int_{-1}^{1} \\sin\\left( x - 0.5 \\right) \\, dx,\n",
    "$$\n",
    "which can be recast as the initial value problem (IVP)\n",
    "$$\n",
    "\\frac{dI}{dx} = \\sin\\left( x - 0.5 \\right), \\quad I(-1) = 0.\n",
    "$$\n",
    "\n",
    "1. **Forward Euler:**  Discretize the ODE using the Forward Euler method. Compute the integral with $n = 10, 100, 1000$ steps, and analyze the convergence rate by comparing with the exact value.\n",
    "    \n",
    "2. **RK4 Method:** Implement the 4th-order Runge-Kutta (RK4) method for the same ODE and compute the integral for $n = 10, 100, 1000$ steps. Compare its accuracy and computational cost with Forward Euler, and determine the observed order of convergence.\n",
    "\n",
    "3.  Explain why ODE solvers like Forward Euler and RK4 are generally not used for numerical integration, addressing issues of stability, accuracy, and computational efficiency."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cbc1b8d-56bc-4c20-bce2-260f716b2e29",
   "metadata": {},
   "source": [
    "## Problem 3 (5+5+5 points)\n",
    "\n",
    "Consider the $5 \\times 5$ matrix:\n",
    "$$\n",
    "A = \\begin{pmatrix}\n",
    "1 & -1 &  &  &  \\\\\n",
    "2 & -1 & 1 &  &  \\\\\n",
    " & 3 & 4 & -2 &  \\\\\n",
    " &  & -2 & 5 & -2 \\\\\n",
    " &  &  & -1 & 3\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "**(a)**\n",
    "\n",
    "Compute (by hand) its LU factorization $A = LU$ via Gaussian elimination (i.e. give both $L$ and $U$).\n",
    "\n",
    "Notice any special pattern to the nonzero entries?\n",
    "\n",
    "**(b)** Explain why you will *always* get this pattern of nonzero entries for LU factorization of an $m \\times m$ tridiagonal matrix, assuming for simplicity that no row swaps are required.   Show that LU factorization of an $m \\times m$ tridiagonal matrix therefore requires $O(m)$ work.\n",
    "\n",
    "**(c)** Given such an LU factorization of a tridiagonal matrix $A$, explain why solving $Ax =b$ then requires $O(m)$ work (instead of $O(m^2)$ for ordinary LU factors) whereas computing $A^{-1}$ requires $O(m^2)$ work (instead of $O(m^3)$)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a640287e-c46c-4ce9-8322-73ebf0d0cae9",
   "metadata": {},
   "source": [
    "## Problem 4 (5+5 points)\n",
    "\n",
    "Consider the matrix\n",
    "$$\n",
    "A = \\begin{pmatrix} 10^{-100} & 1 \\\\ 1 & 1 \\end{pmatrix} \\, .\n",
    "$$\n",
    "Suppose we want to solve $Ax = b$ for the right-hand side $b = \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}$.  The *exact* solution is $x = \\frac{1}{1 - 10^{-100}} \\begin{pmatrix} -1 \\\\ 1 \\end{pmatrix}$.\n",
    "\n",
    "**(a)** Solve this $Ax=b$ for $x$ using Gaussian elimination by hand, but pretend you are doing standard floating-point arithmetic: round the result of every arithmetic operation to 16 significant digits.  Perform *no row swaps*: they are not strictly required since $10^{-100} \\ne 0$.  Show that the resulting approximate solution $\\tilde{x}$ is *completely wrong* compared to the exact solution.\n",
    "\n",
    "**(b)** On a computer, Gaussian elimination uses \"partial pivoting\": it swaps rows to make the pivot element as big in magnitude as possible (looking at all the entries under each pivot before doing elimination under each column).  Carry out this algorithm by hand on the same $Ax=b$ above, again rounding the result of each arithmetic operation.  (Note that any row swaps must be performed on both the matrix and the right-hand-side.)  Show that you get a *much more accurate answer*.\n",
    "\n",
    "This simple example helps explain why partial pivoting is crucial to make Gaussian elimination (i.e. LU factorization) robust in floating-point arithmetic.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.3",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}