{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fun Math" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gauss Quadrature\n", "\n", "\n", "\n", "The following is a simplied case of the fundamental theorem from wikipedia, and its proof. \n", "\n", "If $p_n$ is a polynomial of degree $n$ such that\n", "\\begin{eqnarray}\n", "\\int_{-1}^1 x^k p_n(x)\\,dx=0\\qquad \\forall k=0, 1, 2, \\ldots, n-1, \n", "\\end{eqnarray}\n", "and $\\{x_i\\}$ are roots of $p_n$, then there is a set of weights $\\{w_i\\}$ such that\n", "$$\n", "\\int_{-1}^1 h(x)\\,dx = \\sum_{i=1}^n h(x_i)w_i\n", "$$\n", "for all polynomial $h$ of degree $2n-1$. \n", "\n", "\n", "To prove, first write \n", "$$\n", "\\underbrace{h(x)}_{\\text{of degree } 2n-1} = \\underbrace{p_n(x)}_{n}\\underbrace{q(x)}_{n-1}+\\underbrace{r(x)}_{n-1}. \n", "$$ Thus\n", "$$\n", "\\int_{-1}^1 h(x)\\,dx = \\underbrace{\\int_{-1}^1 p_n(x)q(x)\\,dx}_{= 0} + \\int_{-1}^1 r(x)\\,dx. \n", "$$\n", "The first term of the right hand side vanishes.\n", "On the other hand, we also have\n", "$$\n", "\\sum_{i=1}^n h(x_i)w_i = \\underbrace{\\sum_{i=1}^n p_n(x_i)q(x_i)w_i}_{=0} + \\sum_{i=1}^n r(x_i)w_i\n", "$$\n", "since $\\{x_i\\}$ are roots of $p_n(x)$. \n", "It remains to prove there is a set of weights $\\{w_i\\}$ such that\n", "$$\n", "\\int_{-1}^1 r(x)\\,dx = \\sum_{i=1}^n r(x_i)w_i\n", "$$\n", "for all polynomial $r(x)$ of degree $n-1$, which is true if and only if\n", "$$\n", "\\int_{-1}^1 x^k\\,dx = \\sum_{i=1}^n x_i^k w_i\\qquad\\forall k=0, 1, 2, \\ldots, n-1. \n", "$$\n", "Set $I_k = \\int_{-1}^1 x^k\\,dx$. The above can be written as a linear system\n", "\n", "\n", "$$\n", "\\begin{pmatrix}\n", "1 & 1 & 1 & \\cdots & 1 \\\\\n", "x_1 & x_2 & x_3 & \\cdots & x_n \\\\\n", "x_1^2 & x_2^2 & x_3^2 & \\cdots & x_n^2 \\\\\n", "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "x_1^{n-1} & x_2^{n-1} & x_3^{n-1} & \\cdots & x_n^{n-1}\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "w_1\\\\\n", "w_2\\\\\n", "w_3\\\\\n", "\\vdots\\\\\n", "w_n\n", "\\end{pmatrix}\n", "=\n", "\\begin{pmatrix}\n", "I_0\\\\\n", "I_1\\\\\n", "I_2\\\\\n", "\\vdots\\\\\n", "I_{n-1}\n", "\\end{pmatrix}. \n", "$$\n", "\n", "\n", "The existence of $\\{w_i\\}$ is obvious since the matrix is invertible. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## An Interesting Property of 13\n", "\n", "If you reverse its digits and then find square, or if you reverse the digits of its square, the results are the same: \n", "\\begin{align*}\n", "13^2 &= 169, \\\\\n", "31^2 &= 961. \n", "\\end{align*}\n", "Other examples of numbers with this property are 12 and below: " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n = 20101001\n", "\n", "int(str(n**2)[::-1]) == int(str(n)[::-1])**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following two examples are a bit surprising: " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n = 20211001\n", "\n", "int(str(n**2)[::-1]) == int(str(n)[::-1])**2" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n = 20211002\n", "\n", "int(str(n**2)[::-1]) == int(str(n)[::-1])**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A Somewhat Counterintuitive Result\n", "\n", "Let $X_t$ be the solution of $dX_t = X_t dW_t$ with $X_0 = 1$. Then \n", "\\begin{align*}\n", "X_t = e^{-\\frac{t}{2} + W_t} \\rightarrow 0 \\qquad\\text{ as }t\\rightarrow \\infty, \n", "\\end{align*}\n", "but $X_t$ is a martingale so $E[X_t] = X_0 = 1$ for all $t>0$. \n", "\n", "\n", "Although this might seen counterintuitive at first glance, if one think about a large but finite value of $t$, this is less so. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Closed Form Expression of the Fibonacci Numbers\n", "\n", "\n", "令 $\\lambda_1, \\lambda_2$ 為二次方程式 $a x^2 + b x + c = 0$ 的兩根，則\n", "\\begin{align*}\n", "a \\lambda_1^2 + b \\lambda_1 + c = 0.\n", "\\end{align*}\n", "等號兩邊同乘 $\\lambda_1^n$ 可以得到 \n", "\\begin{align*}\n", "a \\lambda_1^{n+2} + b \\lambda_1^{n+1} + c \\lambda_1^n = 0.\n", "\\end{align*}\n", "同樣的，對 $\\lambda_2$ 也可以得到\n", "\\begin{align*}\n", "a \\lambda_2^{n+2} + b \\lambda_2^{n+1} + c \\lambda_2^n = 0.\n", "\\end{align*}\n", "現在把上面兩個式子作線性組合：\n", "\\begin{align*}\n", "0 &= c_1(a \\lambda_1^{n+2} + b \\lambda_1^{n+1} + c \\lambda_1^n) + c_2(a \\lambda_2^{n+2} + b \\lambda_2^{n+1} + c \\lambda_2^n)\\\\\n", "&= a \\left(c_1\\lambda_1^{n+2}+c_2\\lambda_2^{n+2}\\right) + b \\left(c_1\\lambda_1^{n+1}+c_2\\lambda_2^{n+1}\\right) + c \\left(c_1\\lambda_1^{n}+c_2\\lambda_2^{n}\\right).\n", "\\end{align*}\n", "於是可以得到以下結論：數列 $F_n = c_1\\lambda_1^{n}+c_2\\lambda_2^{n}$ 滿足遞迴式\n", "\\begin{align*}\n", "\\begin{cases}\n", "a F_{n+2} + b F_{n+1} + c F_{n} = 0\\\\\n", "F_0 = c_1 + c_2\\\\\n", "F_1 = c_1\\lambda_1 + c_2\\lambda_2\n", "\\end{cases}\n", "\\end{align*}\n", "其中 $\\lambda_1, \\lambda_2$ 為 $a x^2 + b x + c = 0$ 的兩根。$a x^2 + b x + c = 0$ 就是這個遞迴式的特徵方程。\n", "現在代入費式數列的遞迴定義 $a=1, b=c=-1$ 可以求得\n", "\\begin{align*}\n", "\\lambda_j = \\frac{1\\pm \\sqrt{5}}{2}, \\qquad j=1, 2, \n", "\\end{align*}\n", "再把 $\\lambda_1, \\lambda_2$ 代入\n", "\\begin{align*}\n", "\\begin{cases}\n", "F_0 = c_1 + c_2 = 0\\\\\n", "F_1 = c_1\\lambda_1 + c_2\\lambda_2 = 1\n", "\\end{cases}\n", "\\end{align*}\n", "可以求得 $c_1 = 1/\\sqrt{5}, c_2 = -1/\\sqrt{5}$，所以\n", "\\begin{align*}\n", "F_n = \\frac{1}{\\sqrt{5}}\\left(\\left(\\frac{1+\\sqrt{5}}{2}\\right)^n - \\left(\\frac{1-\\sqrt{5}}{2}\\right)^n\\right). \n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Solution Formula to the Riccati Equation \n", "\n", "Assuming $pq\\neq 0$ and $p\\neq q$, the Riccati equation \n", "\\begin{align*}\n", "\\begin{cases}\n", "y'(t) = a(y(t) - p)(y(t) - q)\\\\\n", "y(0) = 0\n", "\\end{cases}\n", "\\end{align*}\n", "has a closed form \"solution\"\n", "\\begin{align*}\n", "y(t) = pq\\left(\\frac{1-e^{-a(p-q)t}}{p-qe^{-a(p-q)t}}\\right)\n", "\\end{align*}\n", "that people use all the time, although it does not always hold. \n", "Define \n", "\\begin{align*}\n", "t^* = -\\frac{1}{a}\\left(\\frac{\\log p-\\log q}{p-q}\\right). \n", "\\end{align*}\n", "If $t^*$ is real positive, the closed form formula is only the solution to the Riccati equation for $t\\in[0, t^*)$. This is because when $t$ approaches $t^*$, the denominator in the closed form formula goes to 0, making the solution infinity. \n", "Here is a numerical example: When $a=-1, p=2, q=1$, below is a plot of the closed form formula: " ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEMCAYAAAA1VZrrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlm0lEQVR4nO3deZhcZZn+8e/T+753pztJZw9ZCCRAJ2GHAA4RZVNQXBAHNTKg4+g4g8tPnc0Z1FHHGdeMMrgAiuwquwIRMCSdEEISsq+dPensvXc/vz+qEjrQSVdIVZ061ffnuuqq7fSpJ4em7zrPec97zN0RERE5noygCxARkdSnsBARkX4pLEREpF8KCxER6ZfCQkRE+qWwEBGRfqVMWJhZvZk9a2avm9lSM/tM9PUKM3vazFZF78uDrlVEZKCxVDnPwszqgDp3X2hmxcAC4Brgo0Czu99hZl8Ayt399uAqFREZeFJmz8Ldt7r7wujjA8DrwBDgauDn0cV+TiRAREQkiVJmz6I3MxsBzAEmARvdvazXe3vc/S2tKDObBcwCKCwsPGv8+PHJKVZE3rZNe1o41N7N+Nrit72OzXtb2dfSycTBJXGsbGBasGDBLnev7uu9lAsLMysCnge+7u4PmtneWMKit4aGBm9sbExwpSJysj533yLmrWvmhdsvedvr+L8X1/HPv1vG/C9fRnVxbhyrG3jMbIG7N/T1Xsq0oQDMLBt4ALjb3R+Mvrw9ejzj8HGNHUHVJyJxFofvqqcMiuyVrNp+4ORXJseUMmFhZgb8DHjd3b/T661HgZuij28CHkl2bSKSOGYn9/Nja4oAWLXjYByqkWPJCrqAXs4DbgReM7NF0de+BNwB3GdmHwM2AtcHU56IpKLq4lxK87NZqT2LhEqZsHD3F4Bjfce4NJm1iEhyxOOIqZkxrraYJVv2x2Ftciwp04YSkYHJjvkdMXYNw8tZunkfrR3dcahI+qKwEJHQmzqigq4e55VNe4IuJW0pLEQkMPEaun/m8HLMYP46hUWiKCxEJFAnOxoKoDQ/m3GDimnc0HzyK5M+KSxEJC1MG1nBwg176OruCbqUtKSwEJHAxHP+iKkjKjjU0c3rWzWENhEUFiISqDh0oYBIWADMW69WVCIoLEQkLdSW5lFfkU+jwiIhFBYiEph4z2M6dXgF89c3x22UlbxBYSEigbJ4DIeKmjqygl0HO1i361Dc1ikRCgsRSRtTR0SuXtC4XudbxJvCQkQCE+9m0ejqIsoLsnWQOwEUFiISqPg1oSItrYYRFTrInQAKCxFJK9NGVLB+dwtb9rYGXUpaUViISGASMWrpkgk1ADy9bHvc1z2QKSxEJFjx7EMROW4xtqaIJ5Zsi++KBziFhYiknZmTanl53W6aD3UEXUraUFiISGASderc5afW0uPwjFpRcZNSYWFmd5rZDjNb0uu1fzKzzWa2KHq7IsgaRSS+4tyFAuDUwSUMKcvniaVqRcVLSoUFcBcws4/Xv+vuU6K3x5Jck4iEjJkxc1ItL6zaxcH2rqDLSQspFRbuPgfQAGmRgSKBUzjNnFRLR3cPzy7fkbgPGUBSKiyO41NmtjjapioPuhgRiZ94zg3V25nDyqkqylUrKk7CEBY/AkYDU4CtwLf7WsjMZplZo5k17ty5M4nliUgqysww3jFxEM8u30FbZ3fQ5YReyoeFu29392537wH+F5h2jOVmu3uDuzdUV1cnt0gReVs8kX0oIq2olo5u/rxqV0I/ZyBI+bAws7peT68FlhxrWREJn8Q0oSLOGVVJeUE29y/YlMBPGRiygi6gNzO7F7gYqDKzJuBrwMVmNoXIobD1wCeDqk9E4ivR1yjKycrgfVPr+d85a9m6r5W60vzEfmAaS6k9C3f/gLvXuXu2uw9195+5+43ufpq7n+7uV7n71qDrFJH4SdDx7SM+PH04Dtz78sbEflCaS6mwEBGJt/qKAmaMq+He+Zvo6OoJupzQUliISGCSdansG88Zzs4D7TypYbRvm8JCRAJlCT3EHXHR2GrqK/L55dwNCf+sdKWwEJG0l5FhfHj6cOata2bFtgNBlxNKCgsRCUyiz7Po7X0N9eRkZfDLueuT9pnpRGEhIoFK9Giow8oLc7h68mDuX9DEjv1tyfnQNKKwEJEB41OXjKGr2/n+s6uDLiV0FBYiEphkjYY6bHhlIe+fWs+98zayqbkluR8ecgoLERlQPn3JWDLM+O4zK4MuJVQUFiIyoNSW5nHTuSN46JXNrNyukVGxUliISGCS3IU64paLRlOYk8V3ntLeRawUFiISqERd/Oh4Kgpz+PgFI3li6TYWbNDFOWOhsBCRAekTF4xiSFk+tz/wGu1dujhSfxQWIhKYZI+G6q0wN4uvXzuJ1TsO8oNn1wRXSEgoLEQkUMlvQr3h4nE1XHvGEH703GpNA9IPhYWIDGhfefdEivOyuf2BxXT3BLirk+IUFiISoOD/OFcU5vC1KyeyaNNefvbC2qDLSVkKCxEJVACDod7iqsmD+auJg/jmEytoXK/RUX1RWIjIgGdmfOv6yQwpz+e2exay80B70CWlnJQKCzO708x2mNmSXq9VmNnTZrYqel8eZI0iEj9BjoZ6s9L8bH70obPY19rJp+9dSFe3LsHaW0qFBXAXMPNNr30B+KO7jwX+GH0uImkiFdpQh00cXMLXrzmNuWub+eaTK4IuJ6WkVFi4+xzgzQ3Dq4GfRx//HLgmmTWJyMDy3rOG8uGzhzF7zlpdhrWXrKALiMEgd98K4O5bzaymr4XMbBYwC2DYsGFJLE9E3q4U6kId5WtXnsq2fW189ZEllOZnc9XkwUGXFLiU2rM4Ge4+290b3L2huro66HJEJEYW6Gl5fcvOzOD7HzyTqcMr+NxvFvHcih1BlxS4MITFdjOrA4je67+aiCRcXnYmP/1oA2MHFfM3v1rI3LW7gy4pUGEIi0eBm6KPbwIeCbAWEYkjT6XhUH0oycvmFzdPY3BZHjfdOY9nlm0PuqTApFRYmNm9wF+AcWbWZGYfA+4A3mFmq4B3RJ+LSJpIpdFQfakuzuW3t5zLuNpiPvmrBTy4sCnokgKRUge43f0Dx3jr0qQWIiLSS0VhDvd84mxm/aKRz933KrsPdvDxC0YGci2OoKTUnoWIDCyp3YQ6WlFuFnd+dCpXnFbL1x97nb+/71XaOgfOdTAUFiISqDB9N8/LzuT7HziTz152Cg++spnrfvwSTXtagi4rKRQWIhKYFD++3aeMDOMzl43lZzc1sGFXC1f+zws8PQAOfCssRCRYIe37XzphEI986jzqSvP5xC8auf3+xRxs7wq6rIRRWIiIvE2jqot4+Lbz+JuLR3Pfgk2883tz0vZ8DIWFiAQmhF2ot8jJyuD2meP5zaxzALhh9lw+d98idh1Mr2nOFRYiEqhwNqHeatrICp76u4u49eLR/O7VLVz67ef55dwNaTPVucJCRCRO8nMy+ceZ43n8Mxcwoa6Yrzy8hMv/aw5PLt2W8mer90dhISKBCfsf0GMZU1PMvZ84m5/ceBYOfPKXC7jux3/hxdW7QvtvVliISKBCOhiqX2bG5afW8tTfXci/X3saTXta+NBPX+Y9P3qJPy3fHrrQSKnpPkRE0k1WZgYfnD6M95w5hPsXNPGj59Zw812NjK8t5ubzRnLVlMHkZWcGXWa/tGchIpIEedmZfPjs4Tz3DxfzretOB+AfH1jMuXf8if98cgWbmlP7THDtWYhIoNK0C3VM2ZkZXN9Qz3VnDWXu2mbufHEdP3huNT94bjXnj6niA9OGcemEGnKzUmtvQ2EhIhIAM+Oc0ZWcM7qSzXtb+W3jJu6bv4lb715IaX42V5xWxzVTBjN1RAUZGcFHqsJCRAITsmO8CTOkLJ+/u+wUPn3JWF5YvYuHX9nMI4s2c++8jdSV5nH5qbW8c1ItDSMqyAwoOBQWIhKogXRNiP5kZhgXnVLNRadU09LRxdPLtvP7xVu5Z95G7nppPVVFuVwyvppLJwzi/DFVFOYm70+4wkJEJAUV5GRx9ZQhXD1lCAfbu3h2+Q6eWLqNx5ds477GJnKyMpg+soLzx1RxwdhqJtQVJzR4FRYiEhhPi9mhEq8oN4srJw/mysmD6ezuYf76Zv70+g7mrNrJfzy+nP94fDlVRTlMH1XJOaMqOXtUJaOrC+MaHqEJCzNbDxwAuoEud28ItiIRiQc1oU5MdmYG546u4tzRVQBs29fGn1ft5MXVu/jL2t38YfFWACoLczhreDlTR1Rw5vByTh1cclLnc4QmLKJmuPuuoIsQEUkVtaV5XN9Qz/UN9bg7G3a3MHftbuav30Pjhmaeil6YKTvTmFBXwpT6Mk4bUsrpQ8sYXV1IVmZsp9uFLSxEJI1oNFR8mRkjqgoZUVXIDdOGAbBjfxsLN+5l0aa9LNq0hwcWNPGLv2wAID87k/F1xUysK2Hi4JLjrjtMYeHAU2bmwE/cfXbvN81sFjALYNiwYQGUJyJvhwZDJVZNSR4zJ9Uyc1ItAN09zrpdh3ht814WN+1j2Zb9PPrqFu5+eeNx1xOmsDjP3beYWQ3wtJktd/c5h9+MhsdsgIaGBn1fERHpQ2aGMaamiDE1RVx7xlAgMvtv055Whn3j2D8Xmrmh3H1L9H4H8BAwLdiKRORkqQ2VGsyM+oqC4y4TirAws0IzKz78GPgrYEmwVYlIPJjGQ4VCWNpQg4CHomOGs4B73P2JYEsSERk4QhEW7r4WmBx0HSISXzopLzxC0YYSkTSmLlQoKCxERKRfCgsRCYxGQ4WHwkJEAqUuVDgoLEREpF8KCxEJjLpQ4aGwEJFAaW6ocFBYiEhwtGsRGgoLEQmUpvsIB4WFiIj0S2EhIoHRdB/hobAQkUDpAHc4KCxERKRfCgsRCYym+wiPEw6L6IWIMhNRjIgMPGpDhUO/YWFmGWb2QTP7g5ntAJYDW81sqZl9y8zGJr5MEREJUix7Fs8Co4EvALXuXu/uNcAFwFzgDjP7cAJrFJE0pS5UeMRypbzL3L3TzBa6+5mHX3T3ZjPb7O7vNbPsBNYoImlMJ+WFQyx7Ftea2R1AsZlNeNPxitkA7t6ZkOp6MbOZZrbCzFab2RcS/XkiIvKGWPYsXgTygI8D3wHGmdleYAvQmrjS3hANqB8A7wCagPlm9qi7L0vG54tIYriGQ4VGv2Hh7puBX5jZGnd/EcDMKoCRRA52J8M0YLW7r41+/q+Bq4E+w2Jjcwuf+80icrMzyM3KPHKfl51BXvR5XlYmednR17LfeJwffXz4Pjcrg4wM7SaLJIpGQ4VDv2FhZuYRLx5+zd2bgeY3L5OgGgGGAJt6PW8Cpr+pzlnArMPPv3vDGQksR0Ti6VcfD7oC6U8sbahnzewB4BF333j4RTPLAc4HbiIyYuquhFQY/bg+XjsqnNx9NtFjKA0NDd7Y2HjUwj09Tkd3D22d3bR1Ru+7ej3uPPpxa/R5a2c3rR1dvR5H3m/p6D7yvKWzK3Iffe1EYjMvO4PCnCwKcyO3otxMio48jt7yIvfFeVkU5WZH7vOyKMnLojgv8jw/OxPTVzQJmWt/+CJFuVn88mPT+19YEu54f0NiCYuZwM3AvWY2EthL5BhGJvAU8F13X3TSVR5fE1Df6/lQIsdMYpaRYeRlRFpLieTutHf10NLRzaH2LlqjwXKovYuWjm5aOro42N5FS3s3hzoirx1s7+JQexeH2iPL7T7UwYbdLRxo7+JgW2Qd/cnMsCPhUZKfRUleduQWfVxWkE1pfjYl+dmUFeRQmp9NWX7k9ZK8bLXaROS4Yjlm0Qb8EPhhdIhsFdDq7nsTXFtv84Gx0bDaDNwAfDCJnx8zMztyDKSiMCcu6+zq7uFQNFQOtnVxsL2T/W1dHGjr4kBbJwfautjf2nnk+b7o4zU7D7KvtZP9bZ20dfYcp2Yozc+mvCCHsoLIfXlBDhWF2ZQX5lBRkEN5YQ5VRTlUFOZSUZhDSV6W9mREBpBY9iwAMLNngL9391cTWE+f3L3LzD4FPElkj+ZOd1+a7DqCkpWZQWl+BqX5b/90lvau7khwtEbCZF9rJ3tbDt862NPSyZ6WDva2dLJ9fxvLt+5nT0vnMfdqsjONysJcKotyqCrKjdyKc6guyqWmJI/qolyqi3MZVJJLUa6CRfqmwVDhEXNYAP8IfNfMNgBfcvetCaqpT+7+GPBYMj8zneRmZVJTnElNcd4J/VxrRzfNLR00H+xg96F2mg910Hyog10HO9h9sJ1dB9vZdbCDldsPsOtgO53db/2/vyAnk5riXAaV5FFbmkdt9L6uNI/a0nwGl+ZRVZSrVtgApS8S4RBzWLj7QuASM3sv8ISZPQh8092Tcq6FBCM/J5MhOfkMKcvvd1l3Z39rFzsOtLHjQHvkfn87Ow60s31/5PErG/eybX8bHV1Ht8WyM4260nwGl+UxpKyAIeX5DI3e6ssLqCvNIytTkySLBOVE9iywyFeAFcCPgH8DPmFmX3T3XyaiOAkXM6O0IJvSgmzGDio+5nLuzp6WTrbua2Xr3ja27mtly742tuxtZfOeVl5as4tt+9uOalFkZRhDyvMZVlHA8MoCRlQWRm5VhQyrKCAnS0ESRupChceJHLN4ARgFLCUygeBHiZyU9xkzu8DdZx3nx0WOMDMqCnOoKMzh1MGlfS7T0dXDtn1tbNrTwqbmFjbtaWFjcysbdx/i0UVb2N/WdWTZDIP6igJGVxcxurqQMTVFjKkpZuygIkryNG1ZqlMTKhxOZM/iFmBpHyfffdrMXo9jTSLkZGUwrLKAYZUFfb6/51AH63cfYt2uyG3tzkOs2XmQF1bvOqrFVVeaxymDihlfW8z4umLG15YwpqaIbLW0RE7IiRyzWHKct98Vh1pEYlZeGBnOe8aw8qNe7+5xmva0sGr7QVbuOMCq7QdZse0Af1mzm47uSIjkZGZwSm0Rp9aVMmloKZOHljK+tkStrCBoOFRonNAxi2M5PGeTSNAyM4zhlYUMryzksomDjrze2d3D+l2HWLZ1P8u27Gfplv08uWwbv2mMzCKTk5nBhMElnFFfxpnDyzlzWBlDyvI1UicJtInDIS5hIZLqsjMzGDuomLGDirl6yhAgcqC9aU8ri5v2sbhpL4s27eU38zdx10vrAagtyWPqyAqmjaxg+sgKxtYUKTxkwFJYyIBlZtRXFFBfUcC7Tq8DImfLL992gAUb9jB/fTMvr93N716NzCxTVZTLuaMrOW9MJReeUk1daf/DieX41IQKD4WFSC9ZmRlMGlLKpCGl3HTuCNydjc0tzF27m5fW7ObF1bt5NBoe4wYVc9G4amaMq2HqiHKdB/I2aV8tHBQWIsdh9sYxkPdPHYa7s3L7QZ5fuYPnV+7krhfXM3vOWkrzs7lkfA3vmDiIGeNqyM9J7ISV6ULHt8NDYSFyAsyMcbXFjKstZtaFoznU3sWfV+3kqWXb+dPyHTz0ymbyszO5dEIN7z59MBePq074TMdhp+NA4aCwEDkJhblZzJxUx8xJdXR19zBvXTO/f20rTyzZxu8Xb6UkL4srJw/mPWcO5cxhZfrDKKGlsBCJk6zMDM4dU8W5Y6r4l6tO5aU1u3lwYRMPLGzi7pc3Mrq6kA9OH851Zw6ltEBnlgO4DnGHhsJCJAGyMjO48JRqLjylmoPtXTy2eCv3zt/Iv/5+Gd98YjlXTR7MzeePZEJdSdClBk77WuGgsBBJsKLcLN43tZ73Ta1n6ZZ93P3yRh5+ZTO/XdDEBWOrmHXhKM4fU6UWlaQ0jfUTSaJTB5fy79eexktfuIR/uHwcy7cd4MafzeOq77/In5Zv561Tr6W3AfbPDTWFhUgAygpyuG3GGF64fQbfeO9p7Gnp4Oa7Grnmhy/x/MqdQZeXVNqhCgeFhUiAcrMyef/UYTz7+Yu54z2nsetAOzfdOY+P3DmPldsPBF2eyBEpHxZm9k9mttnMFkVvVwRdk0i8ZWdmcMO0Yfzp8xfx/941gVc27uGd3/szX3l4CftaOoMuL2HUhgqPlA+LqO+6+5ToTdfhlrSVm5XJxy8YxfP/MIMPTR/GPfM2cul3nuex17am8fEM9aHCICxhITKgVBTm8C9XT+KR286jtjSXW+9eyKxfLmDbvragS5MBKixh8SkzW2xmd5pZeV8LmNksM2s0s8adOwfWAUJJX5OGlPLwrefxpSvG8+dVO5n5vTk8vWx70GXFTbruK6WjlAgLM3vGzJb0cbsa+BEwGpgCbAW+3dc63H22uze4e0N1dXXyihdJsKzMDGZdOJrH/vYChpTl84lfNPJPjy6lrbM76NLiQqOhwiElTspz98tiWc7M/hf4fYLLEUlJo6qLePDWc/nG4yu488V1zFvXzOyPnMXQ8r6vUy4STymxZ3E8ZlbX6+m1wPGuBS6S1nKzMvnqlRP52U0NbNrTwjU/eJGFG/cEXdbblr4H7dNPyocF8E0ze83MFgMzgM8GXZBI0C6dMIiHbj2Pwtwsbpg998gFmcJIXahwSIk21PG4+41B1yCSisbUFPHwrefxyV8t4G/vfYWte1v55EWjgy5L0lQY9ixE5BjKC3P41cemc+XkwfzH48v5nz+uCrokSVMpv2chIseXk5XBf71/CtkZxrefXklnj/PZy8aGZhbbkJQ54CksRNJAZobxresnk5lh/PcfV9HT43z+8nFBlyVpRGEhkiYyM4xvvPd0MjOM7z+7mkGledx49vCgyzouDYYKD4WFSBrJyDD+7ZpJ7DjQztceWcLQsnxmjK8JuqzjMo2HCgUd4BZJM1mZGfzPB85g4uASbrtnIUs27wu6JEkDCguRNFSYm8WdN02lLD+bm++az/b9qTkBoWt2qNBQWIikqZqSPP7vr6dxoK2Lz/5mEd09qfmHWaOhwkFhIZLGxtUW889XncpLa3bzkzlrgi5HQkxhIZLmrm8YyrtPr+PbT61MuXmkNBoqPBQWImnOzPj6tadRV5rHZ379CvvbUusyrWpDhYPCQmQAKM3P5ns3nMGWvW386++WBV2OhJDCQmSAOGt4OZ+4YBS/XdDEgg3NQZcD6Ep5YaKwEBlAPn3JGOpK8/jKw0tTZnSUTsoLB4WFyABSmJvFl981gWVb93P3yxuCLkcXPwoRhYXIAPOu0+o4b0wl//nkCnYfbA+6HF39KCQUFiIDjJnxz1edSktHN994YnnQ5UhIKCxEBqAxNcXcfP5IfrugidU7DgRWh5pQ4ZESYWFm15vZUjPrMbOGN733RTNbbWYrzOzyoGoUSTe3XDSavKxMfvhssGd2qwsVDikRFsAS4D3AnN4vmtlE4AbgVGAm8EMzy0x+eSLpp6Iwhw9OH8Yjr25h4+6WoMuRFJcSYeHur7v7ij7euhr4tbu3u/s6YDUwLbnViaSvWReOItOMHwc1b5T6UKGREmFxHEOATb2eN0Vfewszm2VmjWbWuHPnzqQUJxJ2g0ryuL5hKPc3NrFtXzDTmIflWuEDXdLCwsyeMbMlfdyuPt6P9fFan99F3H22uze4e0N1dXV8ihYZAG65aDTd7syeszboUiSFJe2yqu5+2dv4sSagvtfzocCW+FQkIgD1FQVcM2UI98zbwG0zRlNZlJu0z1YXKjxSvQ31KHCDmeWa2UhgLDAv4JpE0s7fXDyats4efj1/U/8Lx5maUOGQEmFhZteaWRNwDvAHM3sSwN2XAvcBy4AngNvcvTu4SkXS05iaIqaPrOC3jZs0BYf0KSXCwt0fcveh7p7r7oPc/fJe733d3Ue7+zh3fzzIOkXS2fUN9azf3ULjhuRdIEnBFB4pERYiErwrTqulMCeT3zYmtxWlwVDhoLAQEQAKcrK44rQ6/rB4Ky0dXUGXIylGYSEiR1zfUM+hjm4ee21bUj5PTajwUFiIyBFTR5QzorIgqa0odaHCQWEhIkeYGdedNZSX1zVrvig5isJCRI7ynjOHYgb3L0j83oUGQ4WHwkJEjjK4LJ9zR1fy+9e2JuXzNDdUOCgsROQtLh0/iLU7D6kVJUcoLETkLWaMrwHguZU7Evo5rvFQoaGwEJG3GFlVyPDKAp5bkfjp/tWECgeFhYj0aca4Gl5as4u2Tk3HJgoLETmGi8dV09bZw9y1uxP2GRoNFR4KCxHp09mjKsnNykh8K0p9qFBQWIhIn/KyMzl3dCXPrUjsQW4JB4WFiBzTjPE1rN/dwrpdhxKyfrWhwkNhISLHdPEp0SG0Cdy7MPWhQkFhISLHNKyygFHVhTybhCG0ktoUFiJyXBefUsPctbtp7dAQ2oEsJcLCzK43s6Vm1mNmDb1eH2FmrWa2KHr7cZB1igxE54+tpKOrh8VNexOyfk0NFQ5ZQRcQtQR4D/CTPt5b4+5TkluOiBx2+tAyABY37WP6qMq4rlvX4A6PlAgLd38dNPukSCqqKsplSFk+ixK1Z5GQtUq8pUQbqh8jzewVM3vezC4IuhiRgWhKfVnC2lASDkkLCzN7xsyW9HG7+jg/thUY5u5nAJ8D7jGzkmOsf5aZNZpZ486dGrkhEk+nDy1lU3Mruw+2x3W9akKFR9LaUO5+2dv4mXagPfp4gZmtAU4BGvtYdjYwG6ChoUG/gyJxNLm+DIDFm/cxY1xNXNet7nM4pHQbysyqzSwz+ngUMBZYG2xVIgPPpCGlmMGrm/YGXYoEJCXCwsyuNbMm4BzgD2b2ZPStC4HFZvYqcD9wi7s3B1WnyEBVlJvFmOoiFjfti+t6NRgqPFJlNNRDwEN9vP4A8EDyKxKRN5tcX8azy3fg7nEduajpPsIhJfYsRCT1TR5ayu5DHWze2xp0KRIAhYWIxOTIQe44tqJ0De7wUFiISEzG15aQk5kR94PcGg0VDgoLEYlJTlYGEwaX8KpOzhuQFBYiErPJQ0t5rWkf3T3xaR9pNFR4KCxEJGaTh5ZxqKObtTsPxm2dakOFg8JCRGI2ub4UgEU6OW/AUViISMxGVRWRn53J61sPxGV96kKFh8JCRGKWkWEMryxgw+5DcVyr+lBhoLAQkRMyorKQ9XENCwkDhYWInJARVYVsam6Ny4gojYYKD4WFiJyQEZUFdHT3sCVO035oNFQ4KCxE5ISMqCoEUCtqgFFYiMgJGVF5OCxa4rA29aHCQmEhIidkUEkuedkZrN8Vnz0LdaHCQWEhIifEzBhRWRjn4bOS6hQWInLCRlQWsi4OexYaDRUeCgsROWHDqwriNnxWo6HCQWEhIidsZGVhXIfPSupLibAws2+Z2XIzW2xmD5lZWa/3vmhmq81shZldHmCZIhI1PDoiasNJjohSFyo8UiIsgKeBSe5+OrAS+CKAmU0EbgBOBWYCPzSzzMCqFBEARkbPtVgXh4PcpvFQoZASYeHuT7l7V/TpXGBo9PHVwK/dvd3d1wGrgWlB1Cgib6gpju/wWUl9WUEX0Iebgd9EHw8hEh6HNUVfewszmwXMij5tN7MlCaswfKqAXUEXkUK0Pd5wUtviq9HbyXgF+LeTXEccDfTfjeHHeiNpYWFmzwC1fbz1ZXd/JLrMl4Eu4O7DP9bH8n22Od19NjA7up5Gd2846aLThLbH0bQ93qBtcTRtj2NLWli4+2XHe9/MbgLeDVzqfmT0dRNQ32uxocCWxFQoIiLHkhLHLMxsJnA7cJW79x5e8Shwg5nlmtlIYCwwL4gaRUQGslQ5ZvF9IBd42iJn6Mx191vcfamZ3QcsI9Keus3du2NY3+zElRpK2h5H0/Z4g7bF0bQ9jsFc59uLiEg/UqINJSIiqU1hISIi/Qp1WJjZzOg0IKvN7At9vG9m9t/R9xeb2ZlB1JksMWyPD0W3w2Ize8nMJgdRZzL0ty16LTfVzLrN7Lpk1pdssWwPM7vYzBaZ2VIzez7ZNSZLDP+flJrZ78zs1ei2+Osg6kw57h7KG5AJrAFGATnAq8DENy1zBfA4kfM1zgZeDrrugLfHuUB59PE703V7xLItei33J+Ax4Lqg6w74d6OMyECSYdHnNUHXHeC2+BLwjejjaqAZyAm69qBvYd6zmAasdve17t4B/JrI9CC9XQ38wiPmAmVmVpfsQpOk3+3h7i+5+57o097TqqSbWH43AD4NPADsSGZxAYhle3wQeNDdNwK4e7puk1i2hQPFFhmaWUQkLLoY4MIcFkOATb2e9zUVSCzLpIsT/bd+jMheVzrqd1uY2RDgWuDHSawrKLH8bpwClJvZc2a2wMw+krTqkiuWbfF9YAKRE4BfAz7j7j3JKS91pcp5Fm9HLFOBxDxdSBqI+d9qZjOIhMX5Ca0oOLFsi/8Cbnf3bkv/q+/Esj2ygLOAS4F84C9mNtfdVya6uCSLZVtcDiwCLgFGEzn/68/uvj/BtaW0MIdFLFOBDKTpQmL6t5rZ6cBPgXe6++4k1ZZssWyLBuDX0aCoAq4wsy53fzgpFSZXrP+v7HL3Q8AhM5sDTCZyyYB0Esu2+GvgDo8ctFhtZuuA8Qzw2SPC3IaaD4w1s5FmlkPkuhePvmmZR4GPREdFnQ3sc/etyS40SfrdHmY2DHgQuDENvzH21u+2cPeR7j7C3UcA9wO3pmlQQGz/rzwCXGBmWWZWAEwHXk9ynckQy7bYSGQPCzMbBIwD1ia1yhQU2j0Ld+8ys08BTxIZ4XCnR6YHuSX6/o+JjHK5gsh1MFqIfGNISzFuj68ClUQuIgXQ5Wk4w2aM22LAiGV7uPvrZvYEsBjoAX7q7mk3zX+Mvxv/CtxlZq8RaVvd7u4DedpyQNN9iIhIDMLchhIRkSRRWIiISL8UFiIi0i+FhYiI9EthISIi/VJYiIhIvxQWIiLSL4WFSBKZ2VAze3/QdYicKIWFSHJdCqT1RbgkPekMbpEkMbPziczBtBc4AFzr7usCLUokRgoLkSSKzr/0+XScd0nSm9pQIsk1DlgRdBEiJ0phIZIkZlZJZJr8zqBrETlRCguR5BlJ+l58S9KcwkIkeZYDVWa2xMzODboYkROhA9wiItIv7VmIiEi/FBYiItIvhYWIiPRLYSEiIv1SWIiISL8UFiIi0i+FhYiI9Ov/A+Z96z/Zls6yAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "from pandas import DataFrame\n", "\n", "a = -1\n", "p = 2\n", "q = 1\n", "\n", "t = np.arange(0, 1, 0.01)\n", "y = p*q*((1-np.exp(-a*(p-q)*t))/(p-q*np.exp(-a*(p-q)*t)))\n", "\n", "ax = DataFrame({'t': t, 'y': y}).set_index('t').plot(legend=None)\n", "ax.set(ylim=(-20, 20), xlim=(0, 0.99), xlabel='$t$', ylabel='$y(t)$')\n", "ax.axhline(0, color='k', lw=1)\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example we have $t^* = \\log 2 \\approx 0.693147$. One can see from the plot that only before this critical time the closed form formula is the solution to the initial value problem. \n", "\n", "If $t^*$ is complex or if it is real negative then the closed form formula gives the solution for all $t\\ge 0$. When $p$ and $q$ are both real, the critical value $t^*$ is real positive if and only if $a$ is real negative, which is not the case in the CIR interest rate model. That is why the CIR model never runs into this issue. The Heston model with complex $p$ and $q$, however, is not so lucky. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Newton's Method for Square Root Has Quadratic Convergence\n", "\n", "For square root a simple completing the square will do the trick:\n", "\\begin{align*}\n", "\\frac{x + \\frac{a}{x}}{2} - \\sqrt{a} &= \\frac12\\left(x -2\\sqrt{a} + \\frac{a}{x}\\right)\\\\\n", "&= \\frac12\\left(\\sqrt{x} - \\sqrt{\\frac{a}{x}}\\right)^2\\\\\n", "&= \\frac{1}{2x}\\left(x - \\sqrt{a}\\right)^2, \n", "\\end{align*}\n", "which says the error in the next iteration is the error in this iteration squared divided by $2x$.\n", "\n", "\n", "A more general proof for general functions can be found on [wikipedia.](https://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton's_iterative_method)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checking if an Integer Is a Perfect Square\n", "\n", "The goal here is to prove that Newton's method with integer division (integer division) can be used to determine if an integer is a perfect square. \n", "\n", "This method is discussed (without proof) in a [stackoverflow post](https://stackoverflow.com/questions/2489435/check-if-a-number-is-a-perfect-square), which gives the following procedure. It works for any large integer, for example 123456789**100 which is too large to convert to float for implementations relying on ```math.sqrt```. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n", "CPU times: user 12 ms, sys: 8 ms, total: 20 ms\n", "Wall time: 16.5 ms\n" ] } ], "source": [ "%%time\n", "\n", "def is_square(a):\n", " x = a//2\n", " seen = set([x])\n", " while x*x != a:\n", " x = (x + a//x)//2\n", " if x in seen:\n", " return False\n", " seen.add(x)\n", " return True\n", "\n", "print(is_square(123456789**100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assume that $x\\in\\mathbb N, \\sqrt{a}\\in\\mathbb N$ and $x>\\sqrt{a}$. With integer division, the error after one iteration is really \n", "\\begin{align*}\n", "\\frac{x+\\frac{a - r_1}{x} - r_2}{2} - \\sqrt{a} &= \\frac{x+\\frac{a}{x}}{2} - \\sqrt{a} -\\left( \\frac{r_1}{2x} + \\frac{r_2}{2}\\right)\\\\\n", "&= \\frac{1}{2x}\\left(x - \\sqrt{a}\\right)^2 - \\left( \\frac{r_1}{2x} + \\frac{r_2}{2}\\right), \n", "\\end{align*}\n", "where $r_1$ and $r_2$ are the remainders from the division. Since they must satisfy $r_1 < x$, $r_2\\leq 1$, we have $r_1/(2x) + r_2/2 < 1/2 + 1/2 = 1$.\n", "Now note that this error is an integer, so we can write\n", "\\begin{align*}\n", "\\frac{x+\\frac{a - r_1}{x} - r_2}{2} - \\sqrt{a} &= \\left \\lfloor\\frac{1}{2x}\\left(x - \\sqrt{a}\\right)^2 \\right \\rfloor. \n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation\n", "\n", "This is to confirm the following is correct: \n", "\\begin{align*}\n", "\\text{cov}\\left(\\int_0^T t\\,dW_t, \\int_0^T W_t \\,dt\\right) = \\frac{T^3}{6}, \n", "\\end{align*}\n", "which is 1.3333 if $T=2$ like in the below simulation. \n", "\n", "\n", "Recall that $E(W(t)W(s)) = s\\land t$. Roughly speaking, \n", "\\begin{align*}\n", "E\\left(\\int_0^T t\\,dW_t\\int_0^T W_t dt\\right) &= E\\left(\\int_0^T\\int_0^T s W_t \\,dW_sdt\\right)\\\\\n", "&\\approx E\\left(\\sum_i\\sum_j s_j W(t_i) (W(s_{j+1}) - W(s_j)) \\Delta t\\right)\\\\\n", "&\\approx \\sum_i\\sum_j s_j \\left[E(W(t_i)W(s_{j+1})) - E(W(t_i)W(s_j))\\right] \\Delta t \\\\\n", "&\\approx \\sum_i\\sum_j s_j \\left[ \\Delta s1_{\\{t_i>s_j\\}} \\right] \\Delta t \\\\\n", "&\\approx \\iint\\limits_{t>s} s \\,ds dt = \\frac{T^3}{6}.\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[2.65784801 1.33743398]\n", " [1.33743398 2.68333453]]\n", "CPU times: user 23.2 s, sys: 41.1 ms, total: 23.2 s\n", "Wall time: 24 s\n" ] } ], "source": [ "%%time\n", "import numpy as np\n", "from scipy.stats import norm\n", "from pandas import Series\n", "\n", "n = 10000\n", "T = 2\n", "\n", "sample_size = 10000\n", "sample = []\n", "for i in range(sample_size):\n", " w = Series(norm.rvs(loc=0, scale=np.sqrt(T/n), size=n).cumsum()).shift(1, fill_value=0)\n", " w.index = w.index*T/n\n", " t = w.index\n", " sample.append((w.diff().dropna().values @ t[1:].values, w.sum()*T/n))\n", "\n", "print(np.cov(np.array(sample).T))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inverse Mills Ratio \n", "\n", "We have the conditional expectation\n", "\\begin{align*}\n", "E[X|X>K] &= \\frac{f_X(K)}{1-F_X(K)}, \\\\\n", "E[X|X \\pi^e$。和微分作法一樣可以推出更一般的結論：當指數和底數都大於 $e$ 時，指數比較大者數值就比較大。但像 $2.5^3$ 和 $3^{2.5}$ 這種兩個數分別在 $e$ 的左右時微分作法比不出來，這裡的泰展作法比較好比。不過泰展本來就是在計算數值，有點作弊。\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Square Root of a Complex Number\n", "\n", "If $x\\ge 0$, $y\\ge 0$ then \n", "\\begin{align*}\n", "\\sqrt{x+yi} = \\sqrt{\\frac{\\sqrt{x^2+y^2}+x}{2}} + i \\sqrt{\\frac{\\sqrt{x^2+y^2}-x}{2}}. \n", "\\end{align*}\n", "This can be used to find the numerical values like $\\sqrt{2+17i} = 3.0917 + 2.74929i$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $\\cos(1^\\circ)$ in Closed Form\n", "\n", "用三倍角公式和 $\\cos(60^\\circ)=1/2$ 加上卡當諾公式可以算出 $\\cos(20^\\circ)$ 的 closed form。再利用 $\\sin(18^\\circ)=(-1+\\sqrt 5)/4$ 跟和角公式可以算出 $\\cos(2^\\circ) = \\cos(20^\\circ - 18^\\circ)$。最後用半角公式，可以算出\n", "\n", "\\begin{align*}\n", "\\cos(1^\\circ) = \\sqrt{\\frac{1}{2} \\left(1+\\frac{1}{4} \\sqrt{\\frac{1}{2} \\left(5+\\sqrt{5}\\right)}\n", " \\left(\\sqrt[3]{\\frac{1}{2} \\left(1-i \\sqrt{3}\\right)}+\\sqrt[3]{\\frac{1}{2} \\left(1+i\n", " \\sqrt{3}\\right)}\\right)+\\frac{1}{4} \\left(\\sqrt{5}-1\\right) \\sqrt{1-\\frac{1}{4}\n", " \\left(\\sqrt[3]{\\frac{1}{2} \\left(1-i \\sqrt{3}\\right)}+\\sqrt[3]{\\frac{1}{2} \\left(1+i\n", " \\sqrt{3}\\right)}\\right)^2}\\right)} \n", "\\end{align*}\n", "\n", "但是有 $i$ 在三次根號裡開不出來，所以有這個公式還是沒辦法用古典的方法把數值算出來。$\\cos(1^\\circ)$ 的數值大約是 0.999848。\n", "\n", "為了計算 $\\cos(1^\\circ)$ 的數值，與其設計複數直式開立方來計算上面這個公式的值，還不如設計直式解三次方式程求 $\\cos(20^\\circ)$。有沒有辦法直接直式計算 $\\cos(x^\\circ)$？\n", "\n", "\n", "把 $\\cos(30^\\circ)=\\sqrt{3}/2$ 拿來套半角公式算出 $\\cos(15^\\circ)$，再用 $\\cos(18^\\circ)$ 的值和合角公式就可以算出\n", "\\begin{align*}\n", "\\cos(3^\\circ) = \\frac{1}{8} \\left(\\sqrt{2-\\sqrt{3}} \\left(\\sqrt{5}-1\\right)+\\sqrt{2 \\left(2+\\sqrt{3}\\right) \\left(5+\\sqrt{5}\\right)}\\right)\n", "\\end{align*}\n", "這個是真正的 closed form，沒有任何複數，可以直接用這個公式算出數值。有了這個值，再用卡當諾公式解 $4x^3-3x = \\cos(3^\\circ)$ 就得到\n", "\\begin{align*}\n", "\\cos(1^\\circ) &= \\frac{1}{4} \\left(\\sqrt[3]{\\sqrt{2-\\sqrt{3}} \\left(\\sqrt{5}-1\\right)+\\sqrt{2 \\left(2+\\sqrt{3}\\right) \\left(5+\\sqrt{5}\\right)}-2 \\sqrt{-8+\\sqrt{3}+\\sqrt{15}-\\sqrt{\\frac{1}{2}\n", " \\left(5+\\sqrt{5}\\right)}+\\sqrt{\\frac{1}{2} \\left(25+5 \\sqrt{5}\\right)}}}\\right.\\\\\n", " &+\\left.\\sqrt[3]{\\sqrt{2-\\sqrt{3}} \\left(\\sqrt{5}-1\\right)+\\sqrt{2 \\left(2+\\sqrt{3}\\right)\n", " \\left(5+\\sqrt{5}\\right)}+2 \\sqrt{-8+\\sqrt{3}+\\sqrt{15}-\\sqrt{\\frac{1}{2} \\left(5+\\sqrt{5}\\right)}+\\sqrt{\\frac{1}{2} \\left(25+5 \\sqrt{5}\\right)}}}\\right)\n", "\\end{align*}\n", "這個公式乍看之下沒有複數但實際上是作弊的因為根號裡面有負數。結果就是這個公式的數值還是沒辦法用直式開方法算出來。\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 4 }