{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib widget" ] }, { "cell_type": "code", "execution_count": 118, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "np.set_printoptions(linewidth=100)\n", "from joblib import Parallel, delayed\n", "import sympy as sm\n", "import sympy.abc as sbl\n", "from scipy.sparse import diags\n", "from scipy.integrate import RK45, solve_ivp\n", "import scipy.sparse.linalg as la\n", "import scipy.sparse as sp\n", "import time\n", "from functools import lru_cache\n", "import math\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "from tqdm import tqdm_notebook\n", "from ipywidgets import HBox, IntSlider\n", "\n", "plt.ioff()\n", "plt.clf()" ] }, { "cell_type": "code", "execution_count": 119, "metadata": {}, "outputs": [], "source": [ "class Matrices:\n", " \n", " #A *alpha - beta(D *) +beta(E(1 - 0)alpha) = 0\n", " \n", " def __init__(self, N=10):\n", " \n", " self.N = N\n", " self.shape = (N+1, N+1)\n", " self.h = 1/N\n", " self.domain = np.linspace(0, 1, N+1) \n", " self.beta = -1\n", " self.time_step = self.h**2\n", " self.runs = self.N**2\n", " \n", " self.A = (self.h) * sp.csc_matrix(diags([ [1/6 for i in range(N)], \\\n", " [1/3] + [2/3]*(N-1) + [1/3], \\\n", " [1/6 for i in range(N)] ], \\\n", " [1, 0, -1]), dtype=np.float32)\n", " #-\n", " self.D = (1/self.h) * sp.csr_matrix(diags([ [-1 for i in range(N)], \\\n", " [1] + [2]*(N-1) + [1], \\\n", " [-1 for i in range(N)]], \\\n", " [1, 0, -1] ), dtype=np.float32)\n", " #+ \n", " row = np.array([0, 0, N, N])\n", " col = np.array([0, 1, N-1, N])\n", " data = np.array([1, -1, -1, 1 ])\n", " self.E = (1/self.h) *sp.csr_matrix((data, (row, col)), shape=self.shape, dtype=np.float32) \n", " \n", " self.f = lambda k: ((math.pi**2*self.beta - 1)*(math.pi*h*math.cos(math.pi*h*k) - math.sin(math.pi*h*k) + math.sin(math.pi*h*(k - 1)))*math.exp(t)/(math.pi**2*self.h)) \n", " \n", " #Set initial conditions here\n", " self.alpha = np.sin(np.pi* self.domain).reshape(N+1, 1)\n", " \n", " #Set solution here\n", " self.exact_solution = lambda t: np.sin(np.pi*self.domain).reshape(-1, 1).dot(np.exp(t).reshape(1, -1))\n", " \n", "\n", " \n", "\n", " \n", " def matB(self,i):\n", " if i == 0 :\n", " row = np.array([0, 0, 1, 1])\n", " col = np.array([0, 1, 0, 1])\n", " data = np.array([-1/3, 1/3, -1/6, 1/6])\n", " return sp.csr_matrix((data, (row, col)), shape=self.shape,dtype=np.float32) \n", " elif i == self.N:\n", " row = np.array([self.N-1, self.N - 1, self.N, self.N])\n", " col = np.array([self.N-1, self.N , self.N -1, self.N])\n", " data = np.array([-1/6, 1/6, -1/3, 1/3])\n", " return sp.csr_matrix((data, (row, col)), shape=self.shape, dtype=np.float32) \n", " else:\n", " row = np.array([i-1]*3 + [i]*3 + [i+1]*3)\n", " col = np.array([col for col in range(i-1, i+1 + 1)]*3)\n", " data = np.array([-1/6, 1/6, 0, -1/3, 0, 1/3, 0, -1/6, 1/6 ])\n", " return sp.csr_matrix((data, (row, col)), shape=self.shape, dtype=np.float32) \n", " \n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
half | one | $log(e_i / e_j)$ | |
---|---|---|---|
4 | \n", "0.0926753 | \n", "0.153174 | \n", "10000000 | \n", "
8 | \n", "0.0164509 | \n", "0.0272116 | \n", "10000000 | \n", "
16 | \n", "0.00250125 | \n", "0.00413441 | \n", "10000000 | \n", "
32 | \n", "0.000331428 | \n", "0.000520328 | \n", "10000000 | \n", "
64 | \n", "0.00018354 | \n", "0.000213835 | \n", "10000000 | \n", "
128 | \n", "0.00112933 | \n", "0.000590394 | \n", "10000000 | \n", "