notes

commit 33ce29db834efa19b28dde2f3bfc92e957cc6b93
parent a1e896da15085d307d9cf066aeb70fa02f25096b
Author: miksa <milutin@popovic.xyz>
Date:   Mon, 13 Jun 2022 07:05:40 +0200

prog 12

Diffstat:
A
num_ana/prog/prb12_p.ipynb
 | 
119
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

1 file changed, 119 insertions(+), 0 deletions(-)
diff --git a/num_ana/prog/prb12_p.ipynb b/num_ana/prog/prb12_p.ipynb
@@ -0,0 +1,119 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 556,
+   "id": "3793f001",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from numpy import testing as testing\n",
+    "from scipy import integrate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 532,
+   "id": "38a584af",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def u_n(n):\n",
+    "    if n%2 == 0:\n",
+    "        N = np.array([i for i in range(1, n, 2)])\n",
+    "        return np.r_[-1/N[::-1], 1/N]\n",
+    "    elif n%2 == 1:\n",
+    "        N = np.array([i for i in range(1, n, 2)])\n",
+    "        return np.r_[1/N, 0, -1/N[::-1]]\n",
+    "    return 0\n",
+    "\n",
+    "def h_n(n):\n",
+    "    if n%2 == 0:\n",
+    "        return n//2\n",
+    "    elif n%2 == 1:\n",
+    "        return 0\n",
+    "        \n",
+    "def t_kl(n, k, l):\n",
+    "    if (l-k)%n == h(n):\n",
+    "        return 1\n",
+    "    elif (l-k)%n == h(n) - 1:\n",
+    "        return -1\n",
+    "    else:\n",
+    "        return 0\n",
+    "        \n",
+    "\n",
+    "def fejer_nodes_weights(l):\n",
+    "    n = 2**l\n",
+    "    x = np.array([np.cos(k*np.pi/n) for k in range(n)])\n",
+    "    u = u_n(n)\n",
+    "    T = np.array([[t_kl(n, k, l) for l in range(n)] for k in range(n)])\n",
+    "    w = np.fft.ifft(T@u)\n",
+    "    return x, w"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 569,
+   "id": "ff3bee84",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "poly = lambda x: 10*x**4 + 2*x**3 - x\n",
+    "tan2 = lambda x: np.tan(x**2)\n",
+    "gauss = lambda x: 1/np.sqrt(2*np.pi)*np.exp(x**2/2)\n",
+    "logistic = lambda x: 1/(1+np.exp(-x))\n",
+    "gompertz = lambda x: np.exp(-np.exp(0.1-0.02*x))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 570,
+   "id": "0251393b",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(3.9999999999999996+0j)\t4.000000000000001\n",
+      "(0.7947051658104725+0j)\t0.796828889194331\n",
+      "(0.9534356038045235+0j)\t0.953438269251261\n",
+      "(1+0j)\t1.0\n",
+      "(0.6623136871346119+0j)\t0.662313687134612\n"
+     ]
+    }
+   ],
+   "source": [
+    "x, w = fejer_nodes_weights(3)\n",
+    "\n",
+    "print(poly(x)@w, integrate.quad(poly, -1, 1)[0], sep='\\t')\n",
+    "print(tan2(x)@w, integrate.quad(tan2, -1, 1)[0], sep='\\t')\n",
+    "print(gauss(x)@w, integrate.quad(gauss, -1, 1)[0], sep='\\t')\n",
+    "print(logistic(x)@w, integrate.quad(logistic, -1, 1)[0], sep='\\t')\n",
+    "print(gompertz(x)@w, integrate.quad(gompertz, -1, 1)[0], sep='\\t')"
+   ]
+  }
+ ],
+ "metadata": {
+  "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.10.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}