ex2a.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2A.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next function, we implement the Romberg algorithm. We create an array with size (n,m) which we will fill using the following method.\n",
    "\n",
    "We first calculate R(0,0) using the formula, $$R(0,0) = \\frac{1}{2}  (b-a)  [f(a) + f(b)]$$\n",
    "\n",
    "We then calculate R(n,0) using the formula,\n",
    "$$R(n,0) = \\frac{1}{2}  R(n-1,0) + h  \\sum_{k=1}^{2^{n-1}}f[a + (2k - 1)h]$$\n",
    "where $n ≥ 1$ and $h = \\frac{b-a}{2^n}$.\n",
    "\n",
    "Finally we calculate R(n,m) using the formula,\n",
    "$$R(n,m) = R(n,m-1) + \\frac{1}{4^m - 1}  [R(n,m-1) - R(n - 1, m - 1)]$$\n",
    "\n",
    "For $n ≥ 1$ and $m ≥ 1$\n",
    "\n",
    "All these values are stored in the matrix we created in the beginning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmbrg(a: float, b: float, n: int, m: int, f):\n",
    "    array = np.zeros((n, m))\n",
    "    array[0][0] = 0.5 * (b - a) * (f(a) + f(b))\n",
    "\n",
    "    for i in range(n):\n",
    "        if i != 0:\n",
    "            sum = 0\n",
    "            for k in range((2 ** (i - 1))):\n",
    "                sum += f(a + ((2 * (k + 1)) - 1) * ((b - a) / (2 ** i)))\n",
    "            array[i][0] = (0.5 * array[i - 1][0]) + (\n",
    "                (((b - a) / (2 ** i))) * sum\n",
    "            )\n",
    "\n",
    "        for j in range(min(m, i + 1)):\n",
    "            if i != 0 and j != 0:\n",
    "                array[i][j] = array[i][j - 1] + (1 / ((4 ** j) - 1)) * (\n",
    "                    array[i][j - 1] - array[i - 1][j - 1]\n",
    "                )\n",
    "\n",
    "    return array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the following integral $\\int_{0}^{1} \\frac{4}{1 + x^2}dx$ using the Romberg algorithm. As we can see from our computed matrix, there is a minor difference if we compare the elements from the matrix from the book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_2a = lambda x: 4 / (1 + (x ** 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main():\n",
    "    a=0\n",
    "    b=1\n",
    "    n=5\n",
    "    print(rmbrg(a,b,n,n, f_2a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3.         0.         0.         0.         0.        ]\n",
      " [3.1        3.13333333 0.         0.         0.        ]\n",
      " [3.13117647 3.14156863 3.14211765 0.         0.        ]\n",
      " [3.13898849 3.1415925  3.14159409 3.14158578 0.        ]\n",
      " [3.14094161 3.14159265 3.14159266 3.14159264 3.14159267]]\n"
     ]
    }
   ],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    main()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2A.2.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our experiment we tried to compute the integral $f(x)=\\frac{1}{x+1} $ using the Romberg algorithm. We know that the integral of $f(x)$ is $\\ln(x+1)$. We'll show the resemblance by plotting the integral and the integrated function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_a2 = lambda x: 1 / (x + 1)\n",
    "other_F = lambda x: rmbrg(0, x, 20, 20, f_a2)[-1, -1]\n",
    "F = lambda x: np.log(x + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prog():\n",
    "    array = np.linspace(0, 10, 100)\n",
    "    values_F = np.zeros(100)\n",
    "    values_OTHER_F = np.zeros(100)\n",
    "\n",
    "    for i in range(100):\n",
    "        values_F[i] = F(array[i])\n",
    "        values_OTHER_F[i] = other_F(array[i])\n",
    "        #print(i + 1)\n",
    "\n",
    "    plt.plot(array, values_F)\n",
    "    plt.plot(array, values_OTHER_F)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAghUlEQVR4nO3deXxV9Z3/8dcnBMIqi2EPSUjYQiAkIUDiQnGroFisS8VqXVtal1adtk5bOzrtzHRmfu04XVAQrSKVuiOioihqi6iBBAhLQgJhTSDse/bkfn9/kDqURrjAvTm5976fj0ceuTnnm3PfV5K3J+eec77mnENEREJflNcBREQkMFToIiJhQoUuIhImVOgiImFChS4iEiZU6CIiYeK0hW5mA8zsYzNbb2aFZvZAM2MmmNlhMyto+ng0OHFFROTLRPsxpgH4oXNupZl1AVaY2QfOuaKTxn3inJsc+IgiIuKP0xa6c64CqGh6fNTM1gP9gZML/YzExsa6xMTEc9mEiEjEWbFixT7nXM/m1vmzh/4FM0sEMoBlzazOMbPVwE7gR865wlNtKzExkfz8/DN5ehGRiGdm275snd+FbmadgdeBB51zR05avRJIcM4dM7OrgPnA4Ga2MQ2YBhAfH+/vU4uIiB/8OsvFzNpyvMznOufmnbzeOXfEOXes6fFCoK2ZxTYzbpZzLss5l9WzZ7N/MYiIyFny5ywXA/4IrHfOPf4lY/o0jcPMxjZtd38gg4qIyKn5c8jlQuBbwFozK2ha9jMgHsA5NxO4AbjHzBqAamCq020cRURalD9nuSwF7DRjpgPTAxVKRETOnK4UFREJEyp0EZEwcUbnoYuIyNmpr6tl85pPOVj0MZ0HZjHi4ikBfw4VuohIENTX1bJp9SccLPqYThW5DKpey1CrBeDzmttAhS4i0jo11NexafVSDhR++EWBD2sq8K1R8azteTXRSeNJyLycnD4DgpJBhS4ichYaGxrYUpjLvrUf0GHHZyRXrWWoVQP/V+Btk8eTOPqrJPbqT2ILZFKhi4j4wfl8bC9ZRUXBImLKlpJcVcAgKhkEbI/qT2HslUQnjydx9JUk9o5rkQI/mQpdRORLVGwroXzFe9jWJSQeySeBQyQAO60Xxd0nEDVwPAlZE4nvl0hruDuVCl1EpMnh/bvZtHwh9aV/of+BZcS5CvoC++jG1i6j2Zw4nrjMifQbOIx+XodthgpdRCJWTXUlpfmLObp+MbG7PyO5YROZ5jjmOlDaKZ3yuFvpkzGRhKGZxEa1/st2VOgiEjGcz8eWojz2FLxLx/IlDK5eywiro961YWO7FJYlTKP7iCtITh9PersYr+OeMRW6iIS1fbvK2Lr8bdymjxh4eDlJHCIJ2BY1gNW9ptB+2OUMGnMlw8/r7nXUc6ZCF5GwUldbw8b8DzlS+B69di8luXEzscBBurC5Sxabky4lYczVJMQlk+B12ABToYtIyNu5tYSy5Qtot/UjhlSuJNVqjh9GiRlO7oD7OD99EskjL2B0mzZeRw0qFbqIhJya6ko2LHuPqqJF9Nu3lHjfDvpx/HTCdbETaTf0CgaNu4rhXXt4HbVFqdBFJCTs3FpC2bL5tN/6IUOqVpFmddS6tpR0GMXO+Jvpl3UNAwal0S8EzkYJFhW6iLRK9XW1bMhbzNG179B3zxISfGX0A3ZYb9b0vIb2wycyZNwk0jp18Tpqq6FCF5FW48CeHWz6bD5tShcx6FgeqVRR59pQ0j6N3Phv0G/M1xgwKI3+EbwXfioqdBHxjPP52Lo+j915b9Kt/COG1Bczxhx76U5x90uIHjaJITmTGRkGpxS2BBW6iLSoutoaSpa9S9Wat4jft4SB7GUgsLHNIJYlfIeeo6eQNCKHnmF+RkowqNBFJOgO79/Nxk/fwDa8y9Cjyxhp1VS7dhR3yqIs+T6Scq5jcL8EBnsdNMSp0EUkKCq2lbDts9fosmURQ2vXkmU+9tGNoh6X0W7ENQzLmUxGx85exwwrKnQRCYi/3Sdl9/LX6LVjMcmNm+kLbI0aQF7/b9Ej81oGZ3yFWB1KCRoVuoicNV9jIxvyP+TQynnE7/mIJLebRGeUtBtO7sAHicu+kcRBIzyZ7CESqdBF5IzU19VSnLuQqtXzSd7/F4ZxiDoXzfoOGZQP+h5JF91ISpDmzJRTU6GLyGnV1lSx/tMF1K+dz5BDSxhJJVUuhuIu2WwdNpmhF9/AqAi7zL41UqGLSLNqqo6xful8Gte9wbDDn5Ju1RyhIyVdLyI6dQopF11Lpt7UbFVU6CLyhZqqYxR98ga+dW+QcuRTMqyGQ3SmqMelxKR9nZQLrmFMTHuvY8qXUKGLRLia6krWf/IGjeveIOXwUjKthoN0ofD8K+iQfj3Dsq9ibAjO3hOJVOgiEai+rpaiTxdQV/Aqww4tIcOqOURnCntcToeMG0jJuZqxbdt5HVPOkApdJEI0NjRQvGwRx1a8xNADHzGKYxyhI8XdJxAz6gZSLrxGe+IhToUuEsacz0fp6qXsz51L0u73SeUAVS6GovMuInrUjaRcNIUx7Tt6HVMC5LSFbmYDgDlAH8AHzHLO/e6kMQb8DrgKqALucM6tDHxcEfFH2cbVlC/5E3Hl7zDY7STBRVPYaSxlI24gZfwNZHXu6nVECQJ/9tAbgB8651aaWRdghZl94JwrOmHMJGBw08c4YEbTZxFpIft2lVH68Rx6bJrPkIYN9HfG+pg0lg/9LkMvuYWMHj29jihBdtpCd85VABVNj4+a2XqgP3BioU8B5jjnHJBrZt3MrG/T94pIkFRXHqXw4xdpu+4VUqtXkG0+NrVJInfQQyRdcjup/Qd6HVFa0BkdQzezRCADWHbSqv5A2QlflzctU6GLBJivsZH1ue9RmfcCww9+TJZVs4tY8vp/i74X30ZyShbJXocUT/hd6GbWGXgdeNA5d+Tk1c18i2tmG9OAaQDx8fFnEFNEykvXUf6XZ0koX0AqeznmOlDU/RI6jbmVlOyJ9NFdDCOeX4VuZm05XuZznXPzmhlSDpx4N544YOfJg5xzs4BZAFlZWf9Q+CLy944dOcj6xXPovP5lUuoL6eeMdR0y2ZH6Y1Iv/SZjNUGynMCfs1wM+COw3jn3+JcMWwDcb2YvcfzN0MM6fi5ydpzPR1Hue1Qtm03qob8wxmrZHtWfzwfeT9Jld5EWpwMq0jx/9tAvBL4FrDWzgqZlPwPiAZxzM4GFHD9lsZTjpy3eGfCkImFuz44tbPrgKQZse4NUt4ujrgNrz7+SrhfcwdDMS4jXTPdyGv6c5bKU5o+RnzjGAfcFKpRIpKivq2Xtx6/QpmAOI6ry6GWOwnaj2DXiQVIvu5VxOqQiZ0BXiop4oLx0HWUfPsXgigVkcog99GB53B3EX/YdUpNSvY4nIUqFLtJC6mprWLv4BWLW/IkRtQX0cVGs7ZRN+ejbGTH+OnrpZlhyjlToIkG2Y/N6tn/wJEMr5jOaI1TQk88T72HQV79HRr9Er+NJGFGhiwRBQ33d8WPjK58jrSaf3i6KtZ1yKB9zFyPGf52+OmdcgkCFLhJA+3aVsfHd6SRte5UM9rOHHnweP43kK+8lQ5fhS5Cp0EXOkfP5KMlbTOXSGYw88ldyrJG1MRnsyPwFaZfepGPj0mJU6CJnqabqGGvefYbzi55nWONmjtCRlb1voN/l9zJySLrX8SQCqdBFztCuslK2vPs7UnbOYyzH2BKVyLLURxk56dtk6z7j4iEVuogf/nZYpeqT6aQd/YSeOFZ3voiYC+5heM4kBuoqTmkFVOgip1BfV8vqRbPpuvpphjVs5AidyO97MwmTHiAzYajX8UT+jgpdpBmHD+yl6O3fk7z5BbI4wPao/iwb/nNGXjVNh1Wk1VKhi5xg55Ziti/8DWl7FpBjtayLSadi3H8z8ivXE69zx6WVU6GLABtXLeHIh/9D+tG/EksUq7tdTo/L/4kRI7O9jibiNxW6RCzn87H2r/No8/nvSa1bzVHXgeX9biF58o8Yo4uAJASp0CXiNNTXUbBoNt1XPUla4xb20IPcQQ+Res0PyOnaw+t4ImdNhS4Ro6a6ktVvz2BA0Syy3G62RQ1g+ah/J/2q75Ad097reCLnTIUuYa/y6CHWvvm/DCqdzTgOsSF6CKvGPsqoy24mQW90ShhRoUvYOnxwH0Xzf03KthfI5hjrYtLZdfEPSb1gMqYLgSQMqdAl7Bzev5v1b/w3qWV/JseqKeiYQ4dLH2ZE1qVeRxMJKhW6hI2DeysofuM/SdvxMtlWw8ou4+n61Z+SnnaB19FEWoQKXULeoX27WD/vV6TteJlx1LLqvAn0mPQImcPHeB1NpEWp0CVkHT64j6LX/4O0sj9/UeSxV/0Lo1NGex1NxBMqdAk5x44cZO3r/0XqtjnkUMXKLuPpcfVjjE7J8jqaiKdU6BIyaqqOUfDG/zB04zPkcIRVHS/gvEmPkanL80UAFbqEgIb6OlbM/wMDC6eTzQHWxmSy98rHyMic4HU0kVZFhS6tlvP5WPne8/TK+3+Mczspjk5h7yXTGXnh1V5HE2mVVOjSKhV++g5tP/5XRjdsYGtUPAUXzGDUZVN1QZDIKajQpVXZUpTHkbd+xqjq5ezmfJaP+ndGX3MPidH6URU5Hf2WSKuwb+c2Nr/6M0YfeIdK60hu8gOkX/8wYzt29jqaSMhQoYunqiuPUvDyvzFq22zSaSCvz02kfOOXZJ/f2+toIiFHhS6e8DU2suLtp0hY9WtyOMDKLuPpfd1/kZ2U6nU0kZClQpcWV5y3mKhFP2VMwwY2Rg9m/xUzyRx3pdexREKeCl1azJ4dW9j+8o/JOvIBe+lOXvqvGH3N94jSPclFAuK0hW5mzwKTgT3OuRHNrJ8AvAlsaVo0zzn3ywBmlBBXW1PFypf/g1Gbn2YkPj6Pu5O0qf/KmC7dvI4mElb82UOfDUwH5pxizCfOuckBSSRhZfXHr9JjyaPkuJ2s6nQBvW54nJykFK9jiYSl0xa6c26JmSW2QBYJIxXbStj18kNkVH1KmfVjzYRnyZhwvdexRMJaoI6h55jZamAn8CPnXGFzg8xsGjANID4+PkBPLa1JXW0NK176Jembn6Yb8HnSfWTe9HMGtO/odTSRsBeIQl8JJDjnjpnZVcB8YHBzA51zs4BZAFlZWS4Azy2tSFHue3R8/0fk+MpY1elC+k79HTnxzf4oiEgQnHOhO+eOnPB4oZk9aWaxzrl957ptCQ2H9++m5E8PMvbQQiroScHFT5Fx2VSvY4lEnHMudDPrA+x2zjkzGwtEAfvPOZm0es7nY8U7T5O04t/JcJV83u820m/9FX07dfE6mkhE8ue0xReBCUCsmZUDjwFtAZxzM4EbgHvMrAGoBqY653Q4JcxVbCth74v3kVWTR0n0UA5d+wdyRozzOpZIRPPnLJebT7N+OsdPa5QI4GtsZPkr/01a8W85DyN32D8z5saHaaO7IYp4Tr+F4rftGwqoevUesuuLWNNhDD1vfoLshKFexxKRJip0Oa3GhgbyXvwl6aVPUmvtyEv/FVlfu0eTTYi0Mip0OaXtGwqofmUa2Q0lrOp0IQNuncGYfglexxKRZqjQpVm+xkaWv/wr0kt+R43FkJ/1a0Zf9W3tlYu0Yip0+QcV20o4MPfbZNetoaBjNnHfmkWW9spFWj0VunzB+XzkL5jBsFX/Rhcgb9Qvybr2+9orFwkRKnQBjl/tuem57zDm2F9Z3y6Vrt98ljEDh3kdS0TOgApdKPz0HXp+cD8j3WE+T7qfsbf8QueVi4Qg/dZGsPq6WvJn/5hxO+awI6ovW6+dQ076xV7HEpGzpEKPUDu3lnD0hdvIaShmeY/JjLj7STp27up1LBE5Byr0CLTyvdkMyv0pXZxjxbjHGXvV3V5HEpEAUKFHkNqaKgqeuZ9x+15nQ/QQOn1zDqM1HZxI2FChR4idW4qpnHsr4xo2ktv7ZjLv+i3tYtp7HUtEAkiFHgEKFr9I0tIf0hnHqgueIPurt3odSUSCQIUexhobGlj+3I/I2fEcpW2S6XDLXDJ0iEUkbKnQw9TBvRWUP/NNcmpXsrz7ZNKmzaJ9h05exxKRIFKhh6HS1Z/S6Y3bGeIOsTztF4y9/kGvI4lIC1Chh5n8t2eRmvcIR6wL26+dx9iM8V5HEpEWokIPE40NDSz/4wPkVLxAUbsR9Lr7JQb3GeB1LBFpQSr0MHD08AE2z5xKTvUylsVeR+a0mbRtF+N1LBFpYSr0ELdjcyH1L9xEauMOlqU+wrhvPOx1JBHxiAo9hK1ftog+795NFD6Kr3iecRd9zetIIuIhFXqIyn/rKdLyf8buqF5wyyuMGDTS60gi4jEVeohxPh+5s39CzvanKIwZSdx3X6fr+b29jiUirYAKPYTU19Wy6sk7yDm0kLyuVzLq3jm6H4uIfEGFHiIqjx5i0xM3MLYmj9y4uxl3128016eI/B0VegjYv7ucA7OmMLxhM8vT/pXs6x/yOpKItEIq9FZu55ZifHOuJc63n8KvzGDspVO9jiQirZQKvRXbvG4Z5732DdpSz7bJLzJqzOVeRxKRVkyF3koVL/+Afgtvo4b2HLxpAcNSsryOJCKtnAq9FVq75E2SP/wOB6J60OaON0lMGOp1JBEJATpNopUp+ODPDP3wLna36Uv7775PX5W5iPjptIVuZs+a2R4zW/cl683Mfm9mpWa2xswyAx8zMqxY+Bwjlt7H1rZJ9LjvfWL7xHsdSURCiD976LOBiadYPwkY3PQxDZhx7rEiT/5bT5G+7CE2tkuh7/cX6epPETljpy1059wS4MAphkwB5rjjcoFuZtY3UAEjQd78J8jI/2eKY0YS/4OFdOnaw+tIIhKCAnEMvT9QdsLX5U3L/oGZTTOzfDPL37t3bwCeOvTlzX+C0aseoah9OkkPLKRTl25eRxKREBWIQrdmlrnmBjrnZjnnspxzWT179gzAU4e2/AUzGb3qEQrbpzP4gbfp0KmL15FEJIQFotDLgRPnOosDdgZgu2FtxTvPkLHiJ6yPSWPQD96ifcfOXkcSkRAXiEJfANzWdLZLNnDYOVcRgO2GrZWL/sSo5T+mJGYEA3/wlvbMRSQgTnthkZm9CEwAYs2sHHgMaAvgnJsJLASuAkqBKuDOYIUNB2v/Oo8Rnz3IprZDSPj+23Ts3NXrSCISJk5b6M65m0+z3gH3BSxRGFu/bBGDPvouZdED6HPv23oDVEQCSleKtpDS1Z8St/B29rbpSbdpb9O1h94UFpHAUqG3gPLSdXR/42YqrRMxd73F+b3jvI4kImFIhR5k+3Ztx+ZeRxQ+am9+nd5xyV5HEpEwpUIPoiOH9nP46a/R3XeIPZPnkDA03etIIhLGVOhBUl9Xy7YZ1xPfsJ3SCU8yNOtSryOJSJhToQeB8/koePJ2RtauYlX6L0i75AavI4lIBFChB0Hu7J8w5tC7fD7g24z9+ve9jiMiEUKFHmD5C2aQs/0p8rpeSfadv/Y6johEEBV6ABXnLSZtxc8pbDeSUffOwaL0n1dEWo4aJ0AqtpXQ85272BMVS/9pr9Eupr3XkUQkwqjQA6Dy6CGqn/8GbamncerLdIvt43UkEYlAKvRz5Hw+SmbeSkLjNrZd8oTONRcRz6jQz9GyOT8ns/IT8oY8xMivXOd1HBGJYCr0c7D641cZu+VJ8s+7nHE3/4vXcUQkwqnQz1J56ToG/vUBtkQnkvrd2TqjRUQ8pxY6CzVVx6j78y34MDrc+pJmHBKRVkGFfhbWPD2NJN9Wtn3lt/QbOMzrOCIigAr9jOXNn87Yg+/wef87GHXJjV7HERH5ggr9DGwpymPEql9Q2G4kY+7QZf0i0rqo0P1UXXkUe+1OKq0jve+YS3Tbdl5HEhH5Oyp0P6159j4SfWVUXPK/xPZL8DqOiMg/UKH7YdWi5xm3/01y+9yii4dEpNVSoZ/GrrJSkj7/KRujB5N55+NexxER+VIq9FPwNTay/093Ee0aaT/1Od1BUURaNRX6KSx/+T9JrVtNYdpPGTBopNdxREROSYX+JbaVFJBe8lsKOmQz5us/8DqOiMhpqdCbUV9XS92r36ba2hN3+9O6T4uIhAQ1VTPy5z7K4IaNbBn3b8T2ifc6joiIX1ToJ9lSlMforU+T3+UyMifd6XUcERG/qdBP0NjQQP28e6m0TiTf9oTXcUREzogK/QR5L/8HQxo2sGnMo3Tv2dfrOCIiZ0SF3mTH5kJGbZhOQcccRk+62+s4IiJnzK9CN7OJZlZiZqVm9pNm1k8ws8NmVtD08WjgowaP8/k4+PK9NNCGvt98Qme1iEhIij7dADNrAzwBXAGUA3lmtsA5V3TS0E+cc5ODkDHoVrzzNFm1BSxLfYRxcclexxEROSv+7IqOBUqdc5udc3XAS8CU4MZqOYcP7iNxxa/YED2ErOv+yes4IiJnzZ9C7w+UnfB1edOyk+WY2Woze9fMUgOSrgUUz32Y7u4wNvlx2kSf9g8WEZFWy59Ct2aWuZO+XgkkOOdGAX8A5je7IbNpZpZvZvl79+49o6DBsLHgE7L2ziO/53UMTr/Y6zgiIufEn0IvBwac8HUcsPPEAc65I865Y02PFwJtzSz25A0552Y557Kcc1k9e/Y8h9jnztfYiHv7nzhoXUm5VdPJiUjo86fQ84DBZjbQzNoBU4EFJw4wsz5mZk2PxzZtd3+gwwbSirdmMKRhA1sy/pnzup3vdRwRkXN22oPGzrkGM7sfWAS0AZ51zhWa2fea1s8EbgDuMbMGoBqY6pw7+bBMq1F59BCJBb+hJHoooyd/1+s4IiIB4de7gE2HURaetGzmCY+nA9MDGy141rz0GDkc5MDEZ4lq08brOCIiARFxV9Ds3FJMZvlc8s+7gqFZl3odR0QkYCKu0He9/mMaiSL+Jr0RKiLhJaIKvThvMZnHlrA64XZ69R/odRwRkYCKmEJ3Ph++9x9jH91Iu/ERr+OIiARcxBT6mr+8wvD6dWxKuZdOXbp5HUdEJOAiotAbGxo4b+l/Um59yfz6g17HEREJiogo9JXvzGKgbyu7s35E23YxXscREQmKsC/0utoa4lY9TmmbZDImao5QEQlfYV/oqxZMpy97qbr4EV1EJCJhLawLva62hoTCGRRHpzBy/Ne9jiMiElRhXeirFkynD/uov/hhTSsnImEvbFvuxL3zERdf63UcEZGgC9tC1965iESasGw67Z2LSCQKy0IvWPg0fdhH7YU/1N65iESMsGs7X2MjvdY+xeaoRNK+cr3XcUREWkzYFfqav7xKoq+MA+n3aO9cRCJK2DVeu9w/sIuejNJVoSISYcKq0IvzFjO8fh1bh96pe7aISMQJq0Kv/vhxDtGZtGvu9zqKiEiLC5tCL9u4mlGVn1EcdxMdO3f1Oo6ISIsLm0Lf+f7vaaANgyY/5HUUERFPhEWhHztykNQ977Cm6yXE9hngdRwREU+ERaEXvvsUna2a8ybo2LmIRK6QL3Tn89Gn5E9siB7CkMwJXscREfFMyBf6uqVvkuAr5/BInXcuIpEt5Au9Mfcp9tOVtCvv8DqKiIinQrrQd24pJq0ylw1x1xPTvqPXcUREPBXShb5t8UwckDRRb4aKiIRsoTc2NJC8403WdRxD77hkr+OIiHguZAt93ZJ59OIAvvRbvY4iItIqhGyhN66YwwHOI3XCTV5HERFpFfwqdDObaGYlZlZqZj9pZr2Z2e+b1q8xs8zAR/0/+3aVMfLYZ2zofTXtYtoH86lERELGaQvdzNoATwCTgOHAzWY2/KRhk4DBTR/TgBkBzvl3Shc/Q1trpO8l04L5NCIiIcWfPfSxQKlzbrNzrg54CZhy0pgpwBx3XC7Qzcz6BjgrcPzK0L6bX6M4OoWEYUH9Q0BEJKT4U+j9gbITvi5vWnamYwKiJP9DEnzlHB1+czA2LyISsvwpdGtmmTuLMZjZNDPLN7P8vXv3+pOvWWvaZzH8itvP+vtFRMJRtB9jyoET70kbB+w8izE452YBswCysrL+ofD9MWzsFTD2irP5VhGRsObPHnoeMNjMBppZO2AqsOCkMQuA25rOdskGDjvnKgKcVURETuG0e+jOuQYzux9YBLQBnnXOFZrZ95rWzwQWAlcBpUAVoFsfioi0MH8OueCcW8jx0j5x2cwTHjvgvsBGExGRMxGyV4qKiMjfU6GLiIQJFbqISJhQoYuIhAkVuohImLDjJ6h48MRme4FtZ/ntscC+AMYJBXrNkUGvOTKcy2tOcM71bG6FZ4V+Lsws3zmX5XWOlqTXHBn0miNDsF6zDrmIiIQJFbqISJgI1UKf5XUAD+g1Rwa95sgQlNccksfQRUTkH4XqHrqIiJwk5Ar9dBNWhxszG2BmH5vZejMrNLMHvM7UEsysjZmtMrO3vc7SUsysm5m9ZmbFTf/eOV5nCiYze6jpZ3qdmb1oZmE547uZPWtme8xs3QnLepjZB2a2selz90A8V0gVup8TVoebBuCHzrkUIBu4LwJeM8ADwHqvQ7Sw3wHvOeeGAaMI49dvZv2BHwBZzrkRHL8191RvUwXNbGDiSct+AnzonBsMfNj09TkLqULHvwmrw4pzrsI5t7Lp8VGO/5IHZb7W1sLM4oCrgWe8ztJSzOw8YDzwRwDnXJ1z7pCnoYIvGuhgZtFAR5qZ5SwcOOeWAAdOWjwFeL7p8fPAtYF4rlAr9BabjLo1MrNEIANY5nGUYPst8DDg8zhHS0oC9gLPNR1qesbMOnkdKlicczuA3wDbgQqOz3L2vrepWlTvv83q1vS5VyA2GmqF7tdk1OHIzDoDrwMPOueOeJ0nWMxsMrDHObfC6ywtLBrIBGY45zKASgL0Z3hr1HTMeAowEOgHdDKzW71NFfpCrdD9mow63JhZW46X+Vzn3Dyv8wTZhcDXzGwrxw+pXWpmL3gbqUWUA+XOub/99fUaxws+XF0ObHHO7XXO1QPzgAs8ztSSdptZX4Cmz3sCsdFQK3R/JqwOK2ZmHD+uut4597jXeYLNOfdT51yccy6R4/++Hznnwn7PzTm3Cygzs6FNiy4DijyMFGzbgWwz69j0M34ZYfwmcDMWALc3Pb4deDMQG/VrTtHW4ssmrPY4VrBdCHwLWGtmBU3LftY0z6uEl+8Dc5t2VjYTxpOtO+eWmdlrwEqOn8m1ijC9YtTMXgQmALFmVg48BvwX8IqZ3c3x/7ndGJDn0pWiIiLhIdQOuYiIyJdQoYuIhAkVuohImFChi4iECRW6iEiYUKGLiIQJFbqISJhQoYuIhIn/DySoTr7MHz6nAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    prog()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll compare the accuracy of the Romberg algorithm with emperical data given that:\n",
    "$x = 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true_value=1.0986122886681098\n",
      "i=1: error=0.23472104466522348\n",
      "i=2: error=0.012498822443001156\n",
      "i=3: error=0.0006469705911495982\n",
      "i=4: error=1.825969788793458e-05\n",
      "i=5: error=2.290550196182295e-07\n",
      "i=6: error=1.137817173813005e-09\n",
      "i=7: error=2.0761170560490427e-12\n",
      "i=8: error=1.1102230246251565e-15\n",
      "i=9: error=2.220446049250313e-16\n",
      "i=10: error=2.220446049250313e-16\n",
      "i=11: error=1.1102230246251565e-15\n",
      "i=12: error=2.220446049250313e-16\n",
      "i=13: error=8.881784197001252e-16\n",
      "i=14: error=0.0\n",
      "i=15: error=2.4424906541753444e-15\n",
      "i=16: error=6.661338147750939e-16\n",
      "i=17: error=6.661338147750939e-16\n",
      "i=18: error=6.8833827526759706e-15\n",
      "i=19: error=6.8833827526759706e-15\n",
      "i=20: error=2.1316282072803006e-14\n",
      "i=21: error=4.463096558993129e-14\n",
      "i=22: error=2.3092638912203256e-14\n",
      "i=23: error=3.2862601528904634e-14\n",
      "i=24: error=3.397282455352979e-14\n",
      "i=25: error=4.884981308350689e-15\n"
     ]
    }
   ],
   "source": [
    "def comparison():\n",
    "    x=2\n",
    "    true_value = np.log(x+1)\n",
    "    print(f\"{true_value=}\")\n",
    "    for i in range(1, 26):\n",
    "        error = np.abs(true_value - rmbrg(0, x, i, i, f_a2)[-1, -1])\n",
    "        print(f\"{i=}: {error=}\")\n",
    "        \n",
    "comparison()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see the error decreases as the matrix grows in size. However strangly enough the error increases after the matrix reaches a size of 14x14, which we believe is caused by floating point errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2A.2.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to show why the function $f(x)=\\sqrt{x}$ is a \"bad\" function to integrate using the Romberg algorithm on $[0,1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.5        0.         0.         0.         0.        ]\n",
      " [0.60355339 0.63807119 0.         0.         0.        ]\n",
      " [0.64328305 0.65652626 0.6577566  0.         0.        ]\n",
      " [0.65813022 0.66307928 0.66351615 0.66360757 0.        ]\n",
      " [0.6635812  0.66539819 0.66555278 0.66558511 0.66559287]]\n"
     ]
    }
   ],
   "source": [
    "def f_a3(x: float) -> float:\n",
    "    return x ** 0.5\n",
    "\n",
    "def a2_3():\n",
    "    a = 0.0\n",
    "    b = 1.0\n",
    "    n = 5\n",
    "    print(rmbrg(a, b, n, n, f_a3))\n",
    "    \n",
    "a2_3()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see the function converges very slowly to the true value $(\\frac{2}{3})$. This is because the Romberg algorithm is more accurate, when the function has many continuous derivates. However the derivative of $f(x)$ is not continuous on point $x=0$, that is why we consider this function bad."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2A2.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to compute the following $\\int_{0}^{\\pi} x\\cos(3x)dx$ using the Romberg algorithm for a 6x6 matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-4.9348022   0.          0.          0.          0.          0.        ]\n",
      " [-2.4674011  -1.64493407  0.          0.          0.          0.        ]\n",
      " [-0.36134253  0.340677    0.47305107  0.          0.          0.        ]\n",
      " [-0.24981116 -0.21263404 -0.24952144 -0.26099085  0.          0.        ]\n",
      " [-0.22876078 -0.22174398 -0.22235131 -0.22192004 -0.22176682  0.        ]\n",
      " [-0.22383559 -0.22219387 -0.22222386 -0.22222184 -0.22222302 -0.22222347]]\n"
     ]
    }
   ],
   "source": [
    "def f_a6(x: float) -> float:\n",
    "    return x * np.cos(3 * x)\n",
    "\n",
    "def a2_6():\n",
    "    a=0.0\n",
    "    b=np.pi\n",
    "    n=6\n",
    "    print(rmbrg(a, b, n, n, f_a6))\n",
    "    \n",
    "a2_6()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see the true value of the integrated function is indeed close to $-\\frac{2}{9}$."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "3cbebda45f5c9dced8dbb8fde8f28e3f7a9ec6a4a4afc8fb9ba30fdd41085d7f"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}