notes

uni notes
git clone git://popovic.xyz/notes.git
Log | Files | Refs

prb12_p.ipynb (3070B)


      1 {
      2  "cells": [
      3   {
      4    "cell_type": "code",
      5    "execution_count": 556,
      6    "id": "3793f001",
      7    "metadata": {},
      8    "outputs": [],
      9    "source": [
     10     "import numpy as np\n",
     11     "from numpy import testing as testing\n",
     12     "from scipy import integrate"
     13    ]
     14   },
     15   {
     16    "cell_type": "code",
     17    "execution_count": 532,
     18    "id": "38a584af",
     19    "metadata": {},
     20    "outputs": [],
     21    "source": [
     22     "def u_n(n):\n",
     23     "    if n%2 == 0:\n",
     24     "        N = np.array([i for i in range(1, n, 2)])\n",
     25     "        return np.r_[-1/N[::-1], 1/N]\n",
     26     "    elif n%2 == 1:\n",
     27     "        N = np.array([i for i in range(1, n, 2)])\n",
     28     "        return np.r_[1/N, 0, -1/N[::-1]]\n",
     29     "    return 0\n",
     30     "\n",
     31     "def h_n(n):\n",
     32     "    if n%2 == 0:\n",
     33     "        return n//2\n",
     34     "    elif n%2 == 1:\n",
     35     "        return 0\n",
     36     "        \n",
     37     "def t_kl(n, k, l):\n",
     38     "    if (l-k)%n == h(n):\n",
     39     "        return 1\n",
     40     "    elif (l-k)%n == h(n) - 1:\n",
     41     "        return -1\n",
     42     "    else:\n",
     43     "        return 0\n",
     44     "        \n",
     45     "\n",
     46     "def fejer_nodes_weights(l):\n",
     47     "    n = 2**l\n",
     48     "    x = np.array([np.cos(k*np.pi/n) for k in range(n)])\n",
     49     "    u = u_n(n)\n",
     50     "    T = np.array([[t_kl(n, k, l) for l in range(n)] for k in range(n)])\n",
     51     "    w = np.fft.ifft(T@u)\n",
     52     "    return x, w"
     53    ]
     54   },
     55   {
     56    "cell_type": "code",
     57    "execution_count": 569,
     58    "id": "ff3bee84",
     59    "metadata": {},
     60    "outputs": [],
     61    "source": [
     62     "poly = lambda x: 10*x**4 + 2*x**3 - x\n",
     63     "tan2 = lambda x: np.tan(x**2)\n",
     64     "gauss = lambda x: 1/np.sqrt(2*np.pi)*np.exp(x**2/2)\n",
     65     "logistic = lambda x: 1/(1+np.exp(-x))\n",
     66     "gompertz = lambda x: np.exp(-np.exp(0.1-0.02*x))"
     67    ]
     68   },
     69   {
     70    "cell_type": "code",
     71    "execution_count": 570,
     72    "id": "0251393b",
     73    "metadata": {},
     74    "outputs": [
     75     {
     76      "name": "stdout",
     77      "output_type": "stream",
     78      "text": [
     79       "(3.9999999999999996+0j)\t4.000000000000001\n",
     80       "(0.7947051658104725+0j)\t0.796828889194331\n",
     81       "(0.9534356038045235+0j)\t0.953438269251261\n",
     82       "(1+0j)\t1.0\n",
     83       "(0.6623136871346119+0j)\t0.662313687134612\n"
     84      ]
     85     }
     86    ],
     87    "source": [
     88     "x, w = fejer_nodes_weights(3)\n",
     89     "\n",
     90     "print(poly(x)@w, integrate.quad(poly, -1, 1)[0], sep='\\t')\n",
     91     "print(tan2(x)@w, integrate.quad(tan2, -1, 1)[0], sep='\\t')\n",
     92     "print(gauss(x)@w, integrate.quad(gauss, -1, 1)[0], sep='\\t')\n",
     93     "print(logistic(x)@w, integrate.quad(logistic, -1, 1)[0], sep='\\t')\n",
     94     "print(gompertz(x)@w, integrate.quad(gompertz, -1, 1)[0], sep='\\t')"
     95    ]
     96   }
     97  ],
     98  "metadata": {
     99   "kernelspec": {
    100    "display_name": "Python 3",
    101    "language": "python",
    102    "name": "python3"
    103   },
    104   "language_info": {
    105    "codemirror_mode": {
    106     "name": "ipython",
    107     "version": 3
    108    },
    109    "file_extension": ".py",
    110    "mimetype": "text/x-python",
    111    "name": "python",
    112    "nbconvert_exporter": "python",
    113    "pygments_lexer": "ipython3",
    114    "version": "3.10.5"
    115   }
    116  },
    117  "nbformat": 4,
    118  "nbformat_minor": 5
    119 }