907 lines
27 KiB
Plaintext
907 lines
27 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Populating the interactive namespace from numpy and matplotlib\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%pylab inline"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 112,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import pandas\n",
|
|
"import sympy as Ω\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import math\n",
|
|
"import pylab as p\n",
|
|
"import numpy as np\n",
|
|
"import decimal as d\n",
|
|
"from scipy.misc import derivative\n",
|
|
"import pickle\n",
|
|
"from tqdm import tqdm_notebook as tqdm\n",
|
|
"\n",
|
|
"from dataclasses import dataclass\n",
|
|
"\n",
|
|
"import scipy.sparse\n",
|
|
"import scipy.sparse.linalg\n",
|
|
"\n",
|
|
"import functools"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"Ω.init_printing()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"l_float = 20\n",
|
|
"T_float = 100\n",
|
|
"\n",
|
|
"D_float = 0.6\n",
|
|
"H_float = 4"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 113,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"with open('/Users/abra/projects/study/7semestr/degtyarev/roots128.pickle', 'rb') as file:\n",
|
|
" roots128 = pickle.load(file)\n",
|
|
"\n",
|
|
"roots = roots128\n",
|
|
"\n",
|
|
"def u_elem128(x, t, alpha_n):\n",
|
|
" x = np.float128(x)\n",
|
|
" t = np.float128(t)\n",
|
|
" alpha_n = np.float128(alpha_n)\n",
|
|
" return -4*np.pi**2*(alpha_n*np.cos(alpha_n*x/10) + 40*np.sin(alpha_n*x/10))*(alpha_n*np.sin(2*alpha_n) - 40*np.cos(2*alpha_n) + 40)*np.exp(-3*alpha_n**2*t/500)/(5*(alpha_n**2 - np.pi**2)*(alpha_n**2*(4*alpha_n + np.sin(4*alpha_n)) - 80*alpha_n*(np.cos(4*alpha_n) - 1) + 6400*alpha_n - 1600*np.sin(4*alpha_n)))\n",
|
|
"\n",
|
|
"@functools.lru_cache(maxsize=100000)\n",
|
|
"def u128(x, t, root_n=len(roots)):\n",
|
|
" return sum(u_elem128(x, t, roots[:root_n]))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"@dataclass\n",
|
|
"class Grid:\n",
|
|
" I: int\n",
|
|
" K: int\n",
|
|
" \n",
|
|
" xs: ndarray\n",
|
|
" ts: ndarray\n",
|
|
" \n",
|
|
" h_x: float\n",
|
|
" h_t: float\n",
|
|
" \n",
|
|
" def __init__(self, I, K):\n",
|
|
" self.I = I\n",
|
|
" self.K = K\n",
|
|
"\n",
|
|
" self.xs = linspace(0, l_float, I + 1)\n",
|
|
" self.ts = linspace(0, T_float, K + 1)\n",
|
|
"\n",
|
|
" self.h_x = l_float / I\n",
|
|
" self.h_t = T_float / K\n",
|
|
" \n",
|
|
" self.gamma = D_float * self.h_t / (self.h_x ** 2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def ug_mse(ug1, ug2):\n",
|
|
" return np.mean((ug1 - ug2) ** 2)\n",
|
|
"\n",
|
|
"def ug_mae(ug1, ug2):\n",
|
|
" return np.max(np.abs(ug1 - ug2))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 64,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"n = 5000\n",
|
|
"\n",
|
|
"np.random.seed(0)\n",
|
|
"\n",
|
|
"a = np.random.uniform(0, 0.1, size=n-1)\n",
|
|
"b = np.ones(n)\n",
|
|
"c = np.random.uniform(0, 0.1, size=n-1)\n",
|
|
"\n",
|
|
"z = np.random.uniform(0, 0.1, size=n)\n",
|
|
"\n",
|
|
"A = np.zeros((n, n))\n",
|
|
"\n",
|
|
"for i in range(n):\n",
|
|
" A[i, i] = b[i]\n",
|
|
" \n",
|
|
" if i > 0:\n",
|
|
" A[i, i - 1] = a[i - 1]\n",
|
|
" \n",
|
|
" if i < n - 1:\n",
|
|
" A[i, i + 1] = c[i]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 65,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"1.52 s ± 78.3 ms per loop (mean ± std. dev. of 5 runs, 5 loops each)\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%timeit -n 5 -r 5\n",
|
|
"np.linalg.solve(A, z)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 66,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"A_sp = scipy.sparse.diags([a, b, c], [-1, 0, 1], format='csc')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 67,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"3.1 ms ± 430 µs per loop (mean ± std. dev. of 5 runs, 5 loops each)\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%timeit -n 5 -r 5\n",
|
|
"scipy.sparse.linalg.spsolve(A_sp, z)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## теория"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def calc_theoretical(grid):\n",
|
|
" ug = zeros((grid.K + 1, grid.I + 1), dtype=float)\n",
|
|
" \n",
|
|
" for k, t in enumerate(grid.ts):\n",
|
|
" for i, x in enumerate(grid.xs):\n",
|
|
" ug[k, i] = u128(x, t, 100)\n",
|
|
" \n",
|
|
" return ug"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"58.89331758385291\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"grid = Grid(10, 50)\n",
|
|
"print(sum(calc_theoretical(grid)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## простейшая явная схема"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def calc_explicit(grid):\n",
|
|
" if grid.gamma > 1/2:\n",
|
|
" raise Exception('Bad gamma!')\n",
|
|
" \n",
|
|
" ug = zeros((grid.K + 1, grid.I + 1), dtype=float)\n",
|
|
"\n",
|
|
" ug[0, :] = sin(np.pi * grid.xs / l_float) ** 2 * 4 / 10\n",
|
|
"\n",
|
|
" for k in range(1, grid.K + 1):\n",
|
|
" # middle\n",
|
|
" ug[k, 1 : grid.I] = (\n",
|
|
" grid.gamma * ug[k - 1, 0 : grid.I - 1] +\n",
|
|
" (-2 * grid.gamma + 1) * ug[k - 1, 1 : grid.I] + \n",
|
|
" grid.gamma * ug[k - 1, 2 : grid.I + 1]\n",
|
|
" )\n",
|
|
"\n",
|
|
" # left\n",
|
|
" ug[k, 0] = ug[k, 1] / (H_float * grid.h_x + 1)\n",
|
|
"\n",
|
|
" # right\n",
|
|
" ug[k, grid.I] = ug[k, grid.I - 1] / (H_float * grid.h_x + 1)\n",
|
|
" \n",
|
|
" return ug"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"58.72948448630509\n",
|
|
"4.3733803061715146e-07\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"grid = Grid(10, 50)\n",
|
|
"ug = calc_explicit(grid)\n",
|
|
"print(sum(ug))\n",
|
|
"print(ug_mse(ug, calc_theoretical(grid)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## простейшая неявная схема"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 76,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def calc_implicit(grid):\n",
|
|
" ug = zeros((grid.K + 1, grid.I + 1), dtype=float)\n",
|
|
"\n",
|
|
" ug[0, :] = sin(np.pi * grid.xs / l_float) ** 2 * 4 / 10\n",
|
|
"\n",
|
|
" for k in range(1, grid.K + 1):\n",
|
|
" diags = {\n",
|
|
" -1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" 0: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" +1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" }\n",
|
|
" b = np.zeros(grid.I + 1, dtype=float)\n",
|
|
"\n",
|
|
" # left\n",
|
|
" diags[ 0][0] = H_float + 1 / grid.h_x\n",
|
|
" diags[+1][0] = -1 / grid.h_x\n",
|
|
"\n",
|
|
" # middle\n",
|
|
" for i in range(1, grid.I):\n",
|
|
" diags[-1][i] = grid.gamma\n",
|
|
" diags[ 0][i] = - 2 * grid.gamma - 1\n",
|
|
" diags[+1][i] = grid.gamma\n",
|
|
" b[i] = -ug[k - 1, i]\n",
|
|
"\n",
|
|
" # right\n",
|
|
" diags[ 0][grid.I] = H_float + 1 / grid.h_x\n",
|
|
" diags[-1][grid.I] = -1 / grid.h_x\n",
|
|
"\n",
|
|
" A = scipy.sparse.diags([\n",
|
|
" diags[-1][1:], \n",
|
|
" diags[ 0], \n",
|
|
" diags[+1][:-1],\n",
|
|
" ], [-1, 0, 1], format='csc')\n",
|
|
" \n",
|
|
" ug[k, :] = scipy.sparse.linalg.spsolve(A, b)\n",
|
|
" \n",
|
|
" return ug"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 69,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"59.49556361037571\n",
|
|
"4.006492537092588e-06\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"grid = Grid(10, 50)\n",
|
|
"ug = calc_implicit(grid)\n",
|
|
"print(sum(ug))\n",
|
|
"print(ug_mse(ug, calc_theoretical(grid)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## схема Кранка-Николсон"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 79,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def calc_kn(grid):\n",
|
|
" ug = zeros((grid.K + 1, grid.I + 1), dtype=float)\n",
|
|
"\n",
|
|
" ug[0, :] = sin(np.pi * grid.xs / l_float) ** 2 * 4 / 10\n",
|
|
"\n",
|
|
" for k in range(1, grid.K + 1):\n",
|
|
" diags = {\n",
|
|
" -1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" 0: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" +1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" }\n",
|
|
" b = np.zeros(grid.I + 1, dtype=float)\n",
|
|
"\n",
|
|
" # left\n",
|
|
" diags[ 0][0] = H_float + 1 / grid.h_x\n",
|
|
" diags[+1][0] = -1 / grid.h_x\n",
|
|
"\n",
|
|
" # middle\n",
|
|
" for i in range(1, grid.I):\n",
|
|
" diags[-1][i] = grid.gamma\n",
|
|
" diags[ 0][i] = - 2 * grid.gamma - 2\n",
|
|
" diags[+1][i] = grid.gamma\n",
|
|
" b[i] = -(\n",
|
|
" ug[k - 1][i - 1] * grid.gamma +\n",
|
|
" ug[k - 1][i] * (-2 * grid.gamma + 2) +\n",
|
|
" ug[k - 1][i + 1] * grid.gamma\n",
|
|
" )\n",
|
|
"\n",
|
|
" # right\n",
|
|
" diags[ 0][grid.I] = H_float + 1 / grid.h_x\n",
|
|
" diags[-1][grid.I] = -1 / grid.h_x\n",
|
|
"\n",
|
|
" A = scipy.sparse.diags([\n",
|
|
" diags[-1][1:], \n",
|
|
" diags[ 0], \n",
|
|
" diags[+1][:-1],\n",
|
|
" ], [-1, 0, 1], format='csc')\n",
|
|
" \n",
|
|
" ug[k, :] = scipy.sparse.linalg.spsolve(A, b)\n",
|
|
" \n",
|
|
" return ug"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 15,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"0.048762671501371764\n",
|
|
"5.115545125775158e-07\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"grid = Grid(10, 50)\n",
|
|
"ug = calc_kn(grid)\n",
|
|
"print(mean(ug[-1]))\n",
|
|
"print(ug_mse(ug, calc_theoretical(grid)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## простейшая неявная схема повышенного порядка аппроксимации"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 174,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def calc_implicit_improved(grid):\n",
|
|
" ug = zeros((grid.K + 1, grid.I + 1), dtype=float)\n",
|
|
"\n",
|
|
" ug[0, :] = sin(np.pi * grid.xs / l_float) ** 2 * 4 / 10\n",
|
|
"\n",
|
|
" for k in range(1, grid.K + 1):\n",
|
|
" diags = {\n",
|
|
" -1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" 0: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" +1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" }\n",
|
|
" b = np.zeros(grid.I + 1, dtype=float)\n",
|
|
"\n",
|
|
" # left\n",
|
|
" diags[ 0][0] = (\n",
|
|
" 2 * H_float * grid.gamma * grid.h_x + \n",
|
|
" 2 * grid.gamma + \n",
|
|
" 1\n",
|
|
" )\n",
|
|
" diags[+1][0] = - 2 * grid.gamma\n",
|
|
" b[0] = ug[k - 1, 0]\n",
|
|
"\n",
|
|
" # middle\n",
|
|
" for i in range(1, grid.I):\n",
|
|
" diags[-1][i] = grid.gamma\n",
|
|
" diags[ 0][i] = - 2 * grid.gamma - 1\n",
|
|
" diags[+1][i] = grid.gamma\n",
|
|
" b[i] = -ug[k - 1, i]\n",
|
|
"\n",
|
|
" # right\n",
|
|
" diags[ 0][grid.I] = (\n",
|
|
" 2 * H_float * grid.gamma * grid.h_x + \n",
|
|
" 2 * grid.gamma + \n",
|
|
" 1\n",
|
|
" )\n",
|
|
" diags[-1][grid.I] = - 2 * grid.gamma\n",
|
|
" b[grid.I] = ug[k - 1, grid.I]\n",
|
|
"\n",
|
|
" A = scipy.sparse.diags([\n",
|
|
" diags[-1][1:], \n",
|
|
" diags[ 0], \n",
|
|
" diags[+1][:-1],\n",
|
|
" ], [-1, 0, 1], format='csc')\n",
|
|
" \n",
|
|
" ug[k, :] = scipy.sparse.linalg.spsolve(A, b)\n",
|
|
" \n",
|
|
" return ug"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 177,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"59.486724613775095\n",
|
|
"4.014426575945732e-06\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"grid = Grid(10, 50)\n",
|
|
"ug = calc_implicit_improved(grid)\n",
|
|
"print(sum(ug))\n",
|
|
"print(ug_mse(ug, calc_theoretical(grid)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Кранк-Николсон повышенного порядка аппроксимации"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 174,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def calc_kn_improved(grid):\n",
|
|
" ug = zeros((grid.K + 1, grid.I + 1), dtype=float)\n",
|
|
"\n",
|
|
" ug[0, :] = sin(np.pi * grid.xs / l_float) ** 2 * 4 / 10\n",
|
|
"\n",
|
|
" for k in range(1, grid.K + 1):\n",
|
|
" diags = {\n",
|
|
" -1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" 0: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" +1: np.zeros(grid.I + 1, dtype=float),\n",
|
|
" }\n",
|
|
" b = np.zeros(grid.I + 1, dtype=float)\n",
|
|
"\n",
|
|
" # left\n",
|
|
" #diags[ 0][0] = H_float + 1 / grid.h_x\n",
|
|
" #diags[+1][0] = -1 / grid.h_x\n",
|
|
" \n",
|
|
" diags[ 0][0] = (\n",
|
|
" 2 * H_float * grid.gamma * grid.h_x + \n",
|
|
" 2 * grid.gamma + \n",
|
|
" 1\n",
|
|
" )\n",
|
|
" diags[+1][0] = - 2 * grid.gamma\n",
|
|
" b[0] = ug[k - 1, 0]\n",
|
|
"\n",
|
|
" # middle\n",
|
|
" for i in range(1, grid.I):\n",
|
|
" diags[-1][i] = grid.gamma\n",
|
|
" diags[ 0][i] = - 2 * grid.gamma - 2\n",
|
|
" diags[+1][i] = grid.gamma\n",
|
|
" b[i] = -(\n",
|
|
" ug[k - 1][i - 1] * grid.gamma +\n",
|
|
" ug[k - 1][i] * (-2 * grid.gamma + 2) +\n",
|
|
" ug[k - 1][i + 1] * grid.gamma\n",
|
|
" )\n",
|
|
"\n",
|
|
" # right\n",
|
|
" diags[ 0][grid.I] = H_float + 1 / grid.h_x\n",
|
|
" diags[-1][grid.I] = -1 / grid.h_x\n",
|
|
"\n",
|
|
" A = scipy.sparse.diags([\n",
|
|
" diags[-1][1:], \n",
|
|
" diags[ 0], \n",
|
|
" diags[+1][:-1],\n",
|
|
" ], [-1, 0, 1], format='csc')\n",
|
|
" \n",
|
|
" ug[k, :] = scipy.sparse.linalg.spsolve(A, b)\n",
|
|
" \n",
|
|
" return ug"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 177,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"59.486724613775095\n",
|
|
"4.014426575945732e-06\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"grid = Grid(10, 50)\n",
|
|
"ug = calc_kn_improved(grid)\n",
|
|
"print(sum(ug))\n",
|
|
"print(ug_mse(ug, calc_theoretical(grid)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## сравнение"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 128,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAO4AAAAPBAMAAAAR2O42AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMpndu3bvImbNiRBUq0Qb3U6NAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADlklEQVQ4Eb2TT4hbVRTGf8nkz5u8JBN0J8jEwT906rSDFUrFYlDRjdqhIGqlTIRS0Rlq6KLKbPo2yuCmEV0IVYnakTEgBksLFiRv5a4kDnQYhdKg6LbS/zNJTc8596UdU9deeCffO/d835dz3r0A98gDC/UXiO2erNcNuXBq3zbwJz8skK6fCOGH5S+lsr6voAS4d6VerxkafT32qAENqpdZ0srhZVK3WY/rdrpIK/D6/f5NQxYokgsZDf2dHCZ7Az7jZIUfO/Eq/NyAY1LeJjsHqX6/JiKaBNHzv2CLwqFlUo4FyQO6m+wxNh2vQNuQhWxAusoH8C2zIZfxymSqPE2i5K9cEItHIMdPU10Y/e4VMbOk6SWKPDvkqa8mZSx5ee8TTeUP0axmZEAFQ+61hD/Ni9Bq/NXx10lIc938dS1nVnxrINNMim/OcpY0vYvSwH8sk6oZS9jmK1WtQMKf8jjUCvJXGpmAm3AxROc81ibfS5S1wlmQn7nbV/Wes5q7QiSlLLxK5Os/oXUleQxpuHBjGf+S+ErhyRrNIvmrYx8tyWGLfBOCrN/FvYUoaXrr9+3twJmXZaab1kBKWZzC+aa/3ipv0h+GLGT6x0lvwHiNhccaHK0R22h+Sioc+C7jfBON9D8Cdfiq568HzJPfyYIk76xICmXJRxrM+ZsONK1KERLO/nGtk5Z+xZfMIY4W1bdLXE6wWuD/LUH7lfWMPO6ji2+/wcONI6H/gG0NQiRlrGzhtm9K5F6zGkWk5rw2s3OD4XCw4eY8jSe9mW+yKnWR7/aKJZ3eZblkndX3X7U/9/EOWbsK8jfdJzPWGZxvukKiZ3fPkIVUhVhPz1UrfBf2hHKusr2RKp5cZfMdKUW+v8L3HUs6vV3iGwhxaJkUxjp7/vy1c7I91lVfX3YMWdChr+lFPNaQwe0JEyW8bq58p99mWyq037dgu/Qmf8bprWq/G0OuOCmMJXsHdD9VYuQ6MZmEIRcq8Lxe9i1a81TFKxOvxuT7lqN+x4vCVF/5lRbdELR2XL4vb8ueJDctk8JYkr0i5w8voDVDVnwNWcgfx5thJPTfYJHsVXiS0wV+43QYWbSKQlffw3hvCrThix65mj8vZ8X/9z1yUhgLTvTX2A2/T7wkV2he2IYs3D+5TU7D0v4C+YkHO3ohf5GDPSFh6p3PAzgSQnx1fY305IqM2SVVj6+mCiQfWhS1zcukjLU5+3/iW7ckZCmz5Ei9AAAAAElFTkSuQmCC\n",
|
|
"text/latex": [
|
|
"$$4.3733803061715146e-07$$"
|
|
],
|
|
"text/plain": [
|
|
"4.3733803061715146e-07"
|
|
]
|
|
},
|
|
"execution_count": 128,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"ug_mse(calc_theoretical(grid), calc_explicit(grid))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 169,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def print_comparison(Is, Ks, E, i_p, k_p, i_d, k_d):\n",
|
|
" if (\n",
|
|
" not 0 <= i_p < len(Is) or\n",
|
|
" not 0 <= k_p < len(Ks) or\n",
|
|
" not 0 <= i_p + i_d < len(Is) or\n",
|
|
" not 0 <= k_p + k_d < len(Ks) or\n",
|
|
" E[k_p, i_p] < 0 or\n",
|
|
" E[k_p + k_d, i_p + i_d] < 0\n",
|
|
" ):\n",
|
|
" return\n",
|
|
" \n",
|
|
" error_ratio = E[k_p, i_p] / E[k_p + k_d, i_p + i_d]\n",
|
|
" \n",
|
|
" I1 = Is[i_p]\n",
|
|
" I2 = Is[i_p + i_d]\n",
|
|
" K1 = Ks[k_p]\n",
|
|
" K2 = Ks[k_p + k_d]\n",
|
|
" \n",
|
|
" print(f'h_x x{I2//I1} h_t x{K2//K1} @ ({I1:5},{K1:5}): {error_ratio}')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 170,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def compare_errors(Is, Ks, calc_function):\n",
|
|
" E = np.zeros((len(Ks), len(Is)))\n",
|
|
"\n",
|
|
" for k_p in range(len(Ks)):\n",
|
|
" for i_p in range(len(Is)):\n",
|
|
" grid = Grid(Is[i_p], Ks[k_p])\n",
|
|
"\n",
|
|
" try:\n",
|
|
" error = ug_mae(calc_function(grid), calc_theoretical(grid))\n",
|
|
" except:\n",
|
|
" error = -1\n",
|
|
"\n",
|
|
" E[k_p, i_p] = error\n",
|
|
" \n",
|
|
" for k_p in range(len(Ks)):\n",
|
|
" for i_p in range(len(Is)):\n",
|
|
" print_comparison(Is, Ks, E, i_p, k_p, +1, +1)\n",
|
|
" \n",
|
|
" for k_p in range(len(Ks)):\n",
|
|
" for i_p in range(len(Is)):\n",
|
|
" print_comparison(Is, Ks, E, i_p, k_p, +2, +1)\n",
|
|
" \n",
|
|
" for k_p in range(len(Ks)):\n",
|
|
" for i_p in range(len(Is)):\n",
|
|
" print_comparison(Is, Ks, E, i_p, k_p, +1, +2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 171,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"h_x x2 h_t x4 @ ( 16, 256): 1.919077767270527\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"compare_errors([16, 32, 64, 128, 256], [256, 1024], calc_explicit)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 172,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"h_x x2 h_t x2 @ ( 16, 16): 1.8945839194805507\n",
|
|
"h_x x2 h_t x2 @ ( 32, 16): 1.8480652483857702\n",
|
|
"h_x x2 h_t x2 @ ( 64, 16): 1.825718428756318\n",
|
|
"h_x x2 h_t x2 @ ( 128, 16): 1.8200879583854386\n",
|
|
"h_x x2 h_t x2 @ ( 16, 32): 2.011003900865954\n",
|
|
"h_x x2 h_t x2 @ ( 32, 32): 1.8852164156444766\n",
|
|
"h_x x2 h_t x2 @ ( 64, 32): 1.8445734893796497\n",
|
|
"h_x x2 h_t x2 @ ( 128, 32): 1.8438755692851074\n",
|
|
"h_x x2 h_t x2 @ ( 16, 64): 2.17619657916294\n",
|
|
"h_x x2 h_t x2 @ ( 32, 64): 1.9562383707853712\n",
|
|
"h_x x2 h_t x2 @ ( 64, 64): 1.912181287104185\n",
|
|
"h_x x2 h_t x2 @ ( 128, 64): 1.9058730801460957\n",
|
|
"h_x x2 h_t x2 @ ( 16, 128): 2.3986401861806734\n",
|
|
"h_x x2 h_t x2 @ ( 32, 128): 2.066009495516749\n",
|
|
"h_x x2 h_t x2 @ ( 64, 128): 1.9741429057474853\n",
|
|
"h_x x2 h_t x2 @ ( 128, 128): 1.955899186659026\n",
|
|
"h_x x4 h_t x2 @ ( 16, 16): 1.9431575629889533\n",
|
|
"h_x x4 h_t x2 @ ( 32, 16): 1.8496006164556726\n",
|
|
"h_x x4 h_t x2 @ ( 64, 16): 1.8261475842007147\n",
|
|
"h_x x4 h_t x2 @ ( 16, 32): 2.0776417503976705\n",
|
|
"h_x x4 h_t x2 @ ( 32, 32): 1.8918649575362794\n",
|
|
"h_x x4 h_t x2 @ ( 64, 32): 1.845407456579902\n",
|
|
"h_x x4 h_t x2 @ ( 16, 64): 2.326489502477712\n",
|
|
"h_x x4 h_t x2 @ ( 32, 64): 1.9755444903443895\n",
|
|
"h_x x4 h_t x2 @ ( 64, 64): 1.9125944713395193\n",
|
|
"h_x x4 h_t x2 @ ( 16, 128): 2.8326076010600714\n",
|
|
"h_x x4 h_t x2 @ ( 32, 128): 2.1104815578650395\n",
|
|
"h_x x4 h_t x2 @ ( 64, 128): 1.9752019588093679\n",
|
|
"h_x x2 h_t x4 @ ( 16, 16): 3.45714271998841\n",
|
|
"h_x x2 h_t x4 @ ( 32, 16): 3.396912364589268\n",
|
|
"h_x x2 h_t x4 @ ( 64, 16): 3.3648762818042726\n",
|
|
"h_x x2 h_t x4 @ ( 128, 16): 3.3552270369444783\n",
|
|
"h_x x2 h_t x4 @ ( 16, 32): 3.6798635241236157\n",
|
|
"h_x x2 h_t x4 @ ( 32, 32): 3.569646703206942\n",
|
|
"h_x x2 h_t x4 @ ( 64, 32): 3.514763487476464\n",
|
|
"h_x x2 h_t x4 @ ( 128, 32): 3.512604694405925\n",
|
|
"h_x x2 h_t x4 @ ( 16, 64): 3.8072301037478553\n",
|
|
"h_x x2 h_t x4 @ ( 32, 64): 3.780516278348026\n",
|
|
"h_x x2 h_t x4 @ ( 64, 64): 3.7380285131624063\n",
|
|
"h_x x2 h_t x4 @ ( 128, 64): 3.726890300676519\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"compare_errors([16, 32, 64, 128, 256], [16, 32, 64, 128, 256], calc_implicit)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 173,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"h_x x2 h_t x2 @ ( 16, 16): 1.7867407245342732\n",
|
|
"h_x x2 h_t x2 @ ( 32, 16): 2.0659339289789247\n",
|
|
"h_x x2 h_t x2 @ ( 64, 16): 2.206577164218172\n",
|
|
"h_x x2 h_t x2 @ ( 128, 16): 2.205398898130605\n",
|
|
"h_x x2 h_t x2 @ ( 16, 32): 1.0260466758467883\n",
|
|
"h_x x2 h_t x2 @ ( 32, 32): 1.985994321847054\n",
|
|
"h_x x2 h_t x2 @ ( 64, 32): 2.0974055723617107\n",
|
|
"h_x x2 h_t x2 @ ( 128, 32): 2.1318052621473274\n",
|
|
"h_x x2 h_t x2 @ ( 16, 64): 1.8854175495471748\n",
|
|
"h_x x2 h_t x2 @ ( 32, 64): 1.5501374544674071\n",
|
|
"h_x x2 h_t x2 @ ( 64, 64): 2.0605293935478826\n",
|
|
"h_x x2 h_t x2 @ ( 128, 64): 2.1241455078974423\n",
|
|
"h_x x2 h_t x2 @ ( 16, 128): 2.2073516657031678\n",
|
|
"h_x x2 h_t x2 @ ( 32, 128): 1.8202146022973351\n",
|
|
"h_x x2 h_t x2 @ ( 64, 128): 1.8615590435162546\n",
|
|
"h_x x2 h_t x2 @ ( 128, 128): 2.093178235555532\n",
|
|
"h_x x4 h_t x2 @ ( 16, 16): 1.5416112730228488\n",
|
|
"h_x x4 h_t x2 @ ( 32, 16): 1.9411758005230835\n",
|
|
"h_x x4 h_t x2 @ ( 64, 16): 2.132375626592411\n",
|
|
"h_x x4 h_t x2 @ ( 16, 32): 0.672923482269376\n",
|
|
"h_x x4 h_t x2 @ ( 32, 32): 1.8096548816821525\n",
|
|
"h_x x4 h_t x2 @ ( 64, 32): 2.0030692793517577\n",
|
|
"h_x x4 h_t x2 @ ( 16, 64): 2.09685316770308\n",
|
|
"h_x x4 h_t x2 @ ( 32, 64): 1.3513796667001667\n",
|
|
"h_x x4 h_t x2 @ ( 64, 64): 1.9355394149338634\n",
|
|
"h_x x4 h_t x2 @ ( 16, 128): 5.089679627450866\n",
|
|
"h_x x4 h_t x2 @ ( 32, 128): 2.07031910687416\n",
|
|
"h_x x4 h_t x2 @ ( 64, 128): 1.8247920519286784\n",
|
|
"h_x x2 h_t x4 @ ( 16, 16): 5.410544492760973\n",
|
|
"h_x x2 h_t x4 @ ( 32, 16): 4.755334696107873\n",
|
|
"h_x x2 h_t x4 @ ( 64, 16): 4.925531450048806\n",
|
|
"h_x x2 h_t x4 @ ( 128, 16): 4.865081194253896\n",
|
|
"h_x x2 h_t x4 @ ( 16, 32): 1.4301343984858494\n",
|
|
"h_x x2 h_t x4 @ ( 32, 32): 4.6940709147810455\n",
|
|
"h_x x2 h_t x4 @ ( 64, 32): 4.742894620222057\n",
|
|
"h_x x2 h_t x4 @ ( 128, 32): 4.7415271368204355\n",
|
|
"h_x x2 h_t x4 @ ( 16, 64): 1.4883710757235502\n",
|
|
"h_x x2 h_t x4 @ ( 32, 64): 2.537069293873301\n",
|
|
"h_x x2 h_t x4 @ ( 64, 64): 4.399957273891053\n",
|
|
"h_x x2 h_t x4 @ ( 128, 64): 4.733335280216436\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"compare_errors([16, 32, 64, 128, 256], [16, 32, 64, 128, 256], calc_kn)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 181,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"h_x x2 h_t x2 @ ( 32, 128): 2.038390823677269\n",
|
|
"h_x x2 h_t x2 @ ( 64, 128): 1.9736276882133699\n",
|
|
"h_x x2 h_t x2 @ ( 32, 256): 2.1614872398503246\n",
|
|
"h_x x2 h_t x2 @ ( 64, 256): 2.025328231523804\n",
|
|
"h_x x2 h_t x2 @ ( 32, 512): 2.373745419372439\n",
|
|
"h_x x2 h_t x2 @ ( 64, 512): 2.0878950429120877\n",
|
|
"h_x x2 h_t x2 @ ( 32, 1024): 2.713172301556474\n",
|
|
"h_x x2 h_t x2 @ ( 64, 1024): 2.185936854335897\n",
|
|
"h_x x4 h_t x2 @ ( 32, 128): 2.1195416132190292\n",
|
|
"h_x x4 h_t x2 @ ( 32, 256): 2.330981592342115\n",
|
|
"h_x x4 h_t x2 @ ( 32, 512): 2.7365187229137757\n",
|
|
"h_x x4 h_t x2 @ ( 32, 1024): 3.543382713251337\n",
|
|
"h_x x2 h_t x4 @ ( 32, 128): 3.8282183811233788\n",
|
|
"h_x x2 h_t x4 @ ( 64, 128): 3.8442015881699336\n",
|
|
"h_x x2 h_t x4 @ ( 32, 256): 3.9146871029344097\n",
|
|
"h_x x2 h_t x4 @ ( 64, 256): 3.921190229218736\n",
|
|
"h_x x2 h_t x4 @ ( 32, 512): 3.973114349424527\n",
|
|
"h_x x2 h_t x4 @ ( 64, 512): 3.958967998390918\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"compare_errors([32, 64, 128], [128, 256, 512, 1024, 2048], calc_implicit_improved)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"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.7.2"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|