273 lines
114 KiB
Plaintext
273 lines
114 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 68,
|
||
|
|
"id": "1dc137ec",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import numpy as np\n",
|
||
|
|
"\n",
|
||
|
|
"def f(x):\n",
|
||
|
|
" return np.sin(np.pi * x)\n",
|
||
|
|
" #return 1 - 4*(x-0.5)**2\n",
|
||
|
|
"\n",
|
||
|
|
"def nonuniform_grid(num_points,left_pt,right_pt,is_non_unif=True):\n",
|
||
|
|
" grid = np.array([(float(is_non_unif)/(num_points+1-0.5*(num_points+2))**2)*((i+1)-0.5*(num_points+2))**2 + 1 for i in range(num_points+1)])\n",
|
||
|
|
" grid *= (right_pt-left_pt)/sum(grid)\n",
|
||
|
|
" final_grid = np.zeros(num_points+2)\n",
|
||
|
|
" final_grid[0] = left_pt\n",
|
||
|
|
" for i in range(1,num_points+1) :\n",
|
||
|
|
" final_grid[i] = final_grid[i-1]+grid[i-1]\n",
|
||
|
|
" final_grid[num_points+1] = right_pt\n",
|
||
|
|
" return grid, final_grid"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 69,
|
||
|
|
"id": "6391d3e7",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAIUCAYAAAAHco0LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACNq0lEQVR4nOzdd3hUZfrG8e+TTksoKfTeEnoRsCtNsKHuuqvr2l3XVbc3t7m9/Nzurqtrr6urrroWFGl2ei8JEEJJIKRCIEDqvL8/ZnBjDBAgkzPl/lxXrmRmzjlzTwJz5jnnPc9rzjlERERERETk1MV4HUBERERERCRSqMASERERERFpISqwREREREREWogKLBERERERkRaiAktERERERKSFqMASERERERFpISqwJOSZ2TtmdksrPddXzKzIzCrNrEtrPOepMrNrzOztBrfPNLMtgddwmYfRmsXMfmhmDx/j8e1mNrU1M4mIhKPG+wOvmVkbM3vNzCrM7IUW2uYGMzvvKI+dZ2YFLfE8LSGwH+7vdQ5pfSqwJCQEPkQfDrwZFZnZY2bW/gS30dfMnJnFnWSGeOBPwHTnXHvnXNnJbKe1Oeeecc5Nb3DXL4C/B17DKx7Fajbn3G+cc61SQIuINIeZfcHMlgf2SYVm9qaZneV1ruNpYn/gtc8CGUAX59yVLbFB59ww59w7LbGtYAvsh/O8ziGtTwWWhJJLnHPtgbHAacCPW/n5M4AkYEMrP29L68NJvoaTLU5PVms/n4jI8ZjZt4C/AL/Bv1/oDfwDmOVhrOMK0ffTPsBm51zdqW4oRF+fSJNUYEnIcc7tAt4Ehjd+zMxizOzHZrbDzIrN7EkzSwk8/F7g+77AUcfTm1g/0cz+Yma7A19/Cdw3GNjUYP0FTax75AzZ9Wa208xKzexHx9t24LHzzKzAzL4dyF1oZjc2WPcTwyDN7AYz+6DBbWdmtwWG/u01s/vMzBova2Zbgf7Aa4HfQaKZdTezV82s3MxyzexLDbb7MzN70cyeNrP9wA2BLL8ys48C23jNzLqY2TNmtt/MlplZ36P9/czsusDfp8zMftJwiN9Rnu9nZvZ0g/WvbbD+j472PCIiLS2wP/kFcIdz7iXn3EHnXK1z7jXn3HcDyzTnvf57Dd7rLzOzC81sc+B9+IcNnu/Ie+K/zeyAma00s1ENHr/LzLYGHttoZpc3eOwGM/vQzP5sZuXAzxrtDyzwWLH5h+itNbPhR15nYP9ZEni//bGZxTTY7gdm9ofA/mabmc08xu8sM7Df2Gf+4XuXBu7/OXA38PnAvuTmJtZtY2ZPBJ4nO/B7K2jw+HYz+76ZrQUOmllco31KGzN7PLD+RvwHZ4+W81i/j8fN7AEzmxv4Xb9rZn0arPtXM8sP7ANXmNnZDR6LNf9Q9yN/pxVm1ivwmDOzgQ2e4z4zeyOw3BIzG9BgO9PNbFMg2z8CGTS6I0ypwJKQE3hjuhBY1cTDNwS+zsdfSLQH/h547JzA946B0/KLmlj/R8AkYDQwCpgA/Ng5txkY1mD9yceIeBYwBJgC3G1mmcfadoP1ugIpQA/gZuA+M+t0jOdp7GL8O49RwOeACxov4JwbAOwkcDbQOVcNPAsUAN3xD9f4jZlNabDaLOBFoCPwTOC+q4BrA1kHAIuAx4DOQDbw06YCmlkW/iO91wDdGrzehpp6vobr3x947u5AF6Dn0X4hIiIt7HT8IxlePsYyzXmvT8L/3nc38BDwRWAccDb+/UbD63JmAS/gf3/9F/CK+YesA2wNrJMC/Bx42sy6NVh3IpAHpAO/bpRzOv794mD877efB44Mff9bYJv9gXOB64AbG6w7Ef9Bx1TgHuARM/9BvYYCOV8D3g5k+CrwjJkNcc79FP9ZwH8H9kePNF4f/76kbyDHtMDvqbGrgYvw75sbnwn7Kf591AD8+8Trm1j/iGP9PsC/3/pl4DWv5pP7p2X4/95H/kYvmFlS4LFvBTJeCCQDNwGHjpLhavx/x05ALoG/mZml4t8v/gD/fm8TcMYxXouEOBVYEkpeMbN9wAfAu/jfmBu7BviTcy7POVeJ/83oKmv+0IFrgF8454qdcyX43+iuPcGcP3fOHXbOrQHW4N/BNmfbtYHHa51zs4FK/IVac/3OObfPObcTWIj/zf6YAsXqWcD3nXNVzrnVwMONci1yzr3inPM55w4H7nvMObfVOVeB/2ziVufcvMDO7QVgzFGe8rPAa865D5xzNfg/XLhGyzT1fA3Xf905916gOPwJ4Dve6xQRaSFdgNLjDGlrznv9r51ztcBz+D+w/9U5d8A5twH/EO6RDZZf4Zx7MbD8n/AXZ5MAnHMvOOd2B94v/w1swV/QHbHbOfc351xdE++ntUAHYChgzrls51yhmcXiLy5+EMi0Hfhjo9ewwzn3kHOuHngC/wGzjCZ+F5PwH+j8nXOuxjm3AHgdfyHRHJ8DfuOc2+ucKwDubWKZe51z+U28viPr/9o5V+6cyz/K+kc0+fto8PgbDfY9PwJOP3Imyjn3tHOuLPB7/iOQyP/237fgP1C7yfmtOcY13C8555YG/n09w//24xcCGwJnTesCr2PPMV6LhDgVWBJKLnPOdXTO9XHO3X6UN9PuwI4Gt3cAcTT9xt+UptbvfoI5G77pHcK/c2nOtssa7bQbrnsqz3ss3YFy59yBRrkanlXKb2K9ogY/H27i9tGeu3vD7TnnDvHJI4RHe76jrX+wifVFRIKlDEg9zkG75rzX1wd+PrIfO9Z7aMP3PB//G3FwZMj16sDwu334h86nNrVuY4Fi5+/AfUCRmT1oZsmB9ROaeA0N9wt7GmznyNmYpt73uwP5gdxH29axfOI9n6ZfT7P3GXzyNX3CMX4fn3qewAHccv73d/h2YAhjReDvkML//g698J9pbI5jfX5o+PwO/78DCVMqsCTc7MZ/0ewRvYE6/DuvxmdKmrv+7iBma+62DwJtG9zu2oKZOptZh0a5djW43ZzfW3MV0mBIn5m1wX9EuKFjPV8h/p3VkfXbNrG+iEiwLAKqgMuOsUxL70cavufF4H8P3R24Bugh4E78Xfg6AuuBhkP1jvn+7Zy71zk3Dv8Q+MHAd4FS/GdzGr+GXZ/ewnHtBnoduX7rJLb1iX0GDX4XDTR7nxF47qM6yu/jU89t/i7GnfH/Hc4Gvo//bFmnwN+hgv/9HfLxD1E8FY33nYaGx4c1FVgSbp4Fvmlm/QJvgEfGd9cBJfiHkx1rzolngR+bWVpgzPPdwNPHWP5Es53stlcDV5hZ28AFsZ+6GPhkBIZMfAT81sySzGxkYNvPHHvNk/YicImZnWFmCfiHznxq3P5x1r/YzM4KrP8L9D4lIq0kMCz6bvzXyF4WeE+ON7OZZnZPYLGW3o+MM7MrAmfNvgFUA4uBdviLixIA8zdG+lTzp6Mxs9PMbGLgOqmD+AvH+sDZteeBX5tZh0Ah962TfA1LAtv+XuD3dB5wCf6hkc3xPPADM+tkZj3wF5MnouH6PfFfA9ako/0+GixyYYN9zy+BJYF9aAf8B3JLgDgzuxv/tVZHPAz80swGmd9IO/F5NN8ARgT+zcUBd9ByB1rFA/rgIuHmUeAp/B0Dt+F/g/wqfDyM4dfAh4HhFJOaWP9XwHJgLbAOWBm4ryWcyrb/DNTgPxP3BC1bAF2N/yLi3fgv3P6pc25uC27/Y4HrC76Kf+daCBwAivF/YGju+nfgv4i4ENiLhkmISCtyzv0Jf8HxY/wfqvPxf/B/JbBIS+9H/ov/mqi9+K+DuiJwre5G/NdGLcK/bxgBfHgC203GfwZsL/6hc2XAHwKPfRV/kZGH/7rnf+Hfv56QwLW2lwIz8Z8Z+wdwnXMup5mb+AX+9/htwDz8B9matb8I+Dn+17YNf6ONp46x7LF+H+D/HfwU/9DAcfivtQOYg/9a5M2B9ar45LDEP+Ev9N4G9gOPAG1O4DXgnCsFrsTfUKQMyML/b+xEfhcSQsw/zFNEpOUFzjLuAwY557Z5HEdEJKSY2c+Agc65prrnRR0z+wpwlXPu3FZ+3seBAudca8+/2aTAkMsC4Brn3EKv88iJ0xk
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 864x540 with 2 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"\n",
|
||
|
|
"num_points = 100\n",
|
||
|
|
"h_vals, grid = nonuniform_grid(num_points,0,1)\n",
|
||
|
|
"#grid = np.linspace(0,1,num_points+2)\n",
|
||
|
|
"#h_vals = np.diff(grid)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.figure(figsize=(12,7.5))\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot 1\n",
|
||
|
|
"plt.subplot(1,2,1)\n",
|
||
|
|
"plt.plot(grid)\n",
|
||
|
|
"plt.title(\"Plot of nonuniform grid\")\n",
|
||
|
|
"plt.xlabel(\"Grid point index\")\n",
|
||
|
|
"plt.ylabel(\"Domain value\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot 2\n",
|
||
|
|
"plt.subplot(1,2,2)\n",
|
||
|
|
"plt.plot(grid[:-1],h_vals)\n",
|
||
|
|
"plt.title(\"Comparison of grid spacing\")\n",
|
||
|
|
"plt.xlabel(\"Point in domain\")\n",
|
||
|
|
"plt.ylabel(\"Grid spacing\")\n",
|
||
|
|
"\n",
|
||
|
|
"plt.tight_layout()\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 70,
|
||
|
|
"id": "ac600ea8",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"# Construct the coefficient matrix\n",
|
||
|
|
"A = np.zeros((num_points,num_points))\n",
|
||
|
|
"\n",
|
||
|
|
"for i in range(num_points) :\n",
|
||
|
|
" A[i,i] = -2./(h_vals[i]*h_vals[i+1])\n",
|
||
|
|
" \n",
|
||
|
|
"for i in range(num_points-1) :\n",
|
||
|
|
" A[i+1,i] = 2./(h_vals[i+1]*(h_vals[i+1]+h_vals[i+2]))\n",
|
||
|
|
" A[i,i+1] = 2./(h_vals[i+1]*(h_vals[i]+h_vals[i+1]))"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 71,
|
||
|
|
"id": "7343d63c",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([[ -9343.45827382, 4718.44170984, 0. , ...,\n",
|
||
|
|
" 0. , 0. , 0. ],\n",
|
||
|
|
" [ 4813.73429872, -9724.66514761, 4910.93084888, ...,\n",
|
||
|
|
" 0. , 0. , 0. ],\n",
|
||
|
|
" [ 0. , 5010.04771244, -10121.16880196, ...,\n",
|
||
|
|
" 0. , 0. , 0. ],\n",
|
||
|
|
" ...,\n",
|
||
|
|
" [ 0. , 0. , 0. , ...,\n",
|
||
|
|
" -10121.16880196, 5010.04771244, 0. ],\n",
|
||
|
|
" [ 0. , 0. , 0. , ...,\n",
|
||
|
|
" 4910.93084888, -9724.66514761, 4813.73429872],\n",
|
||
|
|
" [ 0. , 0. , 0. , ...,\n",
|
||
|
|
" 0. , 4718.44170984, -9343.45827382]])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 71,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Check the resulting matrix's form\n",
|
||
|
|
"A"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 72,
|
||
|
|
"id": "be8ce9d4",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"# Solve the discrete system\n",
|
||
|
|
"\n",
|
||
|
|
"bc_left = 0\n",
|
||
|
|
"bc_right = 0\n",
|
||
|
|
"\n",
|
||
|
|
"rhs = f(grid[1:-1])\n",
|
||
|
|
"rhs[0] -= (2/(h_vals[0]*(h_vals[0]+h_vals[1])))*bc_left\n",
|
||
|
|
"rhs[-1] -= (2/(h_vals[-1]*(h_vals[-2]+h_vals[-1])))*bc_right\n",
|
||
|
|
"\n",
|
||
|
|
"u = np.linalg.solve(A,rhs)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 73,
|
||
|
|
"id": "ec5cd1a5",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAHVCAYAAADRmWzIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABdZ0lEQVR4nO3dd3xV9f3H8dcne0IgCXvKEhQn7joQFXCh1jpbrVpHqz87tHVvW0cdrVVrrVVbt1XrqCIu3FoFByjIVPYMK2SP7++P7yUEZCQkud873s/HI4/k3nuSvJN77znv+73fc4455xARERERkW2TEjqAiIiIiEg8U6EWEREREWkBFWoRERERkRZQoRYRERERaQEVahERERGRFlChFhERERFpARVqEQnGzPqYmTOztG38/lPN7LXWzhUtZtbLzNaaWWoTlm3R/yryM35qZu9v6/dH43eZ2VgzO30zt23z/8C8h8xspZl90tzvb4kt/U0ikhhUqEWkxczsB2b2oZmtNrMVZvaBme3Ryr/je2XKOfeYc+6w1vw90eScm+ucy3PO1bX0Z5nZtWb2aGvkCsk5N9o59882+NE/AA4Fejjn9myDnw9s+n5ow79JRGLENo90iIgAmFk74L/Az4GngQxgf6AqZK5YZ2Zpzrna0DmSSG/gO+dcWeggIpJ4NEItIi01EMA594Rzrs45V+Gce805NwnAzFLM7Eozm2NmS83sX2bWflM/yMy+M7NDGl1uPNr3buTzqsg0iX02nlZgZvua2aeRkfJPzWzfRre9bWY3REbPS83sNTMr2kyOqWZ2ZKPLaWa23Mx2M7MsM3vUzErMbFXk93Te2j+p0Qj7WWY2F3hr41F3M+trZu9G8r1hZvdsYtT5VDObG8lzReT7RgGXAydG/jdfbiXLbZGpD9+a2ehG158R+dtLzWy2mZ3b6LaDzGy+mV0UuR8XmdkZjW4vNLMXzWxNZEpFvy38/s3+DyP3088iX6dGsi43s9nAERv9nPZm9o9IlgVmdqNtYvqMmZ0FPADsE/n/XLfxYyeynDOz/pGvH478/1+O/D/+Z2b9Gi27g5m9HnlHZomZXb65+2Gjv2mzz4dGj4fTN76PRSS2qVCLSEtNB+rM7J9mNtrMOmx0+08jH8OB7YA84O5t+D0HRD4XRKZJfNT4RjPrCLwM3AUUAncAL5tZYaPFTgHOADrhR9Iv3szvegI4udHlkcBy59xnwOlAe6Bn5PecB1Q04+84EBgc+Zkbexz4JPJzrwV+sollfgAMAkYAV5vZYOfcq8AfgKci/5udt/D79wKmAUXArcA/zMwity0FjgTa4f9Pd5rZbo2+twv+b+8OnAXc0+j+vgeoBLoCZ0Y+Nqep/8OzI3l2BYYBx290+z+BWqB/ZJnDgJ9t/EOcc/+I/I6PIv+fa7aQrbGTgeuADsBM4PcAZpYPvAG8CnSL/P43m3g//JStPx++dx83Ma+IBKJCLSIt4pxbgy8ADvg7sCwyUrlu1PZU4A7n3Gzn3FrgMuAka8HOdZtxBDDDOfeIc67WOfcE8A1wVKNlHnLOTXfOVeCnp+yymZ/1OHC0meVELp8SuQ6gBl8C+0dG5CdG/gdNda1zriySoYGZ9QL2AK52zlU7594HXtzE918XeRfgS+BLYEvleVPmOOf+Hpm3/U98Ae4M4Jx72Tk3y3nvAK/hp++sUwNc75yrcc69AqwFBkVGhX8YyV7mnPsq8rM3p6n/wxOAPznn5jnnVgA3rbsh8vgaDfwq8juXAncCJzXz/7ElzznnPolMzXmM9Y+XI4HFzrnbnXOVzrlS59z/mvgzm/J8aOl9LCJRpkItIi3mnJvqnPupc64HsCN+1O5PkZu7AXMaLT4Hv//GVqdJNNPGv2fd7+re6PLiRl+X40cHv8c5NxOYChwVKdVHs75QPwKMA540s4VmdquZpTcj57wt5F/hnCvfyrJN+hvMbP/ItIO1Zvb1pr6/0e/Ki3zPaDP7ODKNYRVwOH4ke52SjeZ9r/v9xfj7tHHeje+Lxpr6P+y2hZ/ZG0gHFkWmjawC/oZ/96G1bO5/3ROYtY0/synPhybdxyISO1SoRaRVOee+AR7GF2uAhfjys04v/Nv0Szbx7WVATqPLXRr/6K386o1/z7rftWAr37c566Z9jAGmREo2kdHZ65xzQ4B98aOVpzXj527u71gEdGw0Kg6+uG3Tz3XOvReZdpDnnNtha99sZpnAs8BtQGfnXAHwCmBb+r6IZfj7tHHeXpsN2vT/4aIt/Mx5+B1fi5xzBZGPdk35WyM2eKyZWZctLLuxeWx+jnhzH6dbej6ISJxQoRaRFjGz7SM7qvWIXO6JL6IfRxZ5Avi1+R3u8lg/x3RTR7j4Av/2d7qZbTxndhlQj593uimvAAPN7BTzOxGeCAzBH4FkWzyJn5P7c9aPTmNmw81saGSawxr89IUWH/bOOTcHmABca2YZZrYPG05X2ZolQB8z29b1egaQSaQcm99ZsUmHJIxMH3kOnz3HzIbg50lvUjP+h08DF5pZj8hc7Usb/c5F+Ckpt5tZu8jOfv3M7MAm/bV+KsUOZraLmWXh56w31X+BLmb2KzPLNLN8M9srctvW7ofmPB9EJE6oUItIS5Xid3T7n5mV4Yv0V8BFkdsfxL/F/y7wLX7Htf/bzM+6Cj/ytxK/M1hDkY1MT/g98EHkLf69G3+jc64EP9J5EVAC/A440jm3fFv+qEhh+wg/gvpUo5u6AM/gi+BU4B3gUQAzu8/M7tuW3xdxKrAPPv+Nkd/b1MMP/jvyucTMPmvuL3bOlQIX4kvsSvy88U3N4d6cC/BTExbj36F4aAvLbvZ/uJG/46eGfAl8hi/tjZ2GfyEwJZL5Gfyc8K1yzk0HrsfvXDgDaPJJaCL/q0PxL3gWR75/eOTmrd0PzXk+iEicMOe29u6UiIiEYGZPAd8046gUIiISgEaoRURihJntEZm2kBI5pvEY4PnAsUREZCt0pkQRkdjRBT+toRCYD/zcOfd52EgiIrI1mvIhIiIiItICmvIhIiIiItICKtQiIiIiIi0Q13Ooi4qKXJ8+fULHEBEREZEEN3HixOXOueJN3RbXhbpPnz5MmDAhdAwRERERSXBmNmdzt2nKh4iIiIhIC6hQi4iIiIi0gAq1iIiIiEgLxPUcahERERGJrpqaGubPn09lZWXoKG0iKyuLHj16kJ6e3uTvUaEWERERkSabP38++fn59OnTBzMLHadVOecoKSlh/vz59O3bt8nfpykfIiIiItJklZWVFBYWJlyZBjAzCgsLmz36rkItIiIiIs2SiGV6nW3521SoRURERERaQIVaRERERKQFVKhFREREJGFUVFRw4IEHUldXt9llqqurOeCAA6itrW2V36lCLSIiIiIJ48EHH+S4444jNTV1s8tkZGQwYsQInnrqqVb5nSrUIiIiIhJ39t57b7777jsAFixYwLBhwwB47LHHGDNmTMNyw4cP5/XXXwfgyiuv5MILLwTgmGOO4bHHHmuVLDoOtYiIiIhsk+te+popC9e06s8c0q0d1xy1wxaXcc4xd+5cevfuDcCkSZMYOnQo1dXVzJ49mz59+qzPeN11XH311SxdupTPP/+cF198EYAdd9yRTz/9tFUyR2WE2sweNLOlZvbVZm43M7vLzGaa2SQz2y0auUREREQk/sycOZO+ffs2HOJuXaFevnw5BQUFGyx7wAEH4Jzjjjvu4Mknn2yYCpKamkpGRgalpaUtzhOtEeqHgbuBf23m9tHAgMjHXsBfI59FREREJEZtbSS5rUyePJmhQ4c2XJ4wYQLnnnsu2dnZ3zspy+TJk1m0aBFFRUXk5+dvcFtVVRVZWVktzhOVEWrn3LvAii0sMgb4l/M+BgrMrGs0somIJCznoLYKylfA6gWwat6GH6sX+Ntqq/yyIiJxYsWKFWRnZwMwdepUXn75ZYYOHUqHDh2oq6trKNWLFi3i1FNP5YUXXiA3N5dx48Y1/IySkhKKi4tJT09vcZ5YmUPdHZjX6PL8yHWLwsQREYlhFSth1dxIKZ4Hq+fD2qVQXgLly31JrlgJNeXg6pv2My0
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 864x540 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Add boundary conditions\n",
|
||
|
|
"\n",
|
||
|
|
"U = np.append(np.array([bc_left]),np.append(u,bc_right))\n",
|
||
|
|
"RHS = np.append(np.array([bc_left]),np.append(rhs,bc_right))\n",
|
||
|
|
"\n",
|
||
|
|
"plt.figure(figsize=(12,7.5))\n",
|
||
|
|
"plt.plot(grid,U,label=\"$u(x)$\")\n",
|
||
|
|
"plt.plot(grid,RHS,label=\"$f(x)$\")\n",
|
||
|
|
"plt.title(\"Solution vs. right-hand side function\")\n",
|
||
|
|
"plt.xlabel(\"$x$\")\n",
|
||
|
|
"plt.ylabel(\"$u(x)$\")\n",
|
||
|
|
"plt.legend()\n",
|
||
|
|
"\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 74,
|
||
|
|
"id": "6a80c620",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAHVCAYAAAA3sKeMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABUH0lEQVR4nO3dd3iUVf7+8fukExJSSKEFAoTeIVRRsGDBgg3brm1trGWta11Xd93VXdt37d21ixVFxV4QpIbeIfTQElpITyZzfn8k+mMRJIFMzpT367pymZl5MrnjE8Kdh8+cY6y1AgAAAFB3Ya4DAAAAAIGGEg0AAADUEyUaAAAAqCdKNAAAAFBPlGgAAACgnijRAAAAQD1RogGgERljlhhjRvr4cxxpjFnhy88BAKHOsE40ANSfMWadpHRJ1ZJKJE2SdJ21triBP8fl1tpvGuo5AQANgyvRAHDoTrXWxknqL2mgpL84zuMXjDERdbmvvs8BAP6EEg0Ah8lau0nS55J6SpIx5rTasY3dxpgfjDHdfj7WGLPOGHNc7fv3GmPeNca8Zowpqv2Y7NrHXpfUVtInxphiY8ytxphXjTE31z7e2hhjjTFX197OMsbsNDVGGmPy9vqctxljNtV+jhXGmGNr7w8zxtxujFltjNlRmyX5QF+nMeYUY8z82q9rmjGm9z5f123GmIWSSmrzWGPMZcaYDZK+q/18fzHGrDfG5Nd+3Qm1H5+57/ENcnIAwEco0QBwmIwxGZJGS5pnjOks6W1JN0hKVc2YxyfGmKgDfPhpksZLSpQ0UdKTkmStvVDSBtVe7bbWPihpsqSRtR83QtKa2v9K0lGSpth9ZvSMMV0kXStpoLU2XtIJktbVPvwnSafXPkcrSbskPXWAr7G/pJclXSWpuaTnJE00xkTvddj5kk6u/Vo8e+XsVvt5L6l9O1pSB0lxP3+9e9n7eADwWwFZoo0xL9dexVjcQM9XXXt1Zb4xZmJDPCeAkPCRMWa3pKmqKbj3SzpX0mfW2q+ttVWSHpbURNKwAzzHVGvtJGtttaTXJfX5jc83WdKRxpgw1ZTmByUdUfvYiNrH91UtKVpSd2NMpLV2nbV2de1jV0m6y1qbZ62tkHSvpLMPMEpxhaTnrLUzrbXV1tpXJVVIGrLXMY9bazdaa8v2uu9ea21J7X2/k/SotXZN7ez4HZLO2+fz7X08APitgCzRkl6RdGIDPl+ZtbZv7dtpDfi8AILb6dbaRGttO2vt1bXFr5Wk9T8fYK31StooqfUBnmPrXu+XSoo50DxwbfktltRX0pGSPpW0ufZq835LtLU2VzVXxe+VlG+MGW+MaVX7cDtJE2rHM3ZLWqaa0p2+n0/fTtLNPx9be3xG7df7s437+bi97/uf/ze170fs8/n29xwA4HcCskRba3+UtHPv+4wxHY0xXxhj5hhjphhjujqKByC0bVZN4ZQkGWOMasrmpkN4rv0tnzRZ0tmSompnsSdLukhSkqT5+30Sa9+y1g6vzWUl/bv2oY2STqr9ReDnt5ja593XRkn/3OfYWGvt2wfJu/d9//P/RjUz3x5J2w7yHADgdwKyRB/A86pZXmqApFskPV2Pj40xxuQYY2YYY073SToAoeJdSScbY441xkRKulk1Yw/TDuG5tqlmdnhvk1Uz4/xj7e0fJF2nmrGQ6n2fwBjTxRhzTO3scrmkMtVcbZakZyX90xjTrvbYVGPMmANkeUHSOGPM4NoXLzY1xpxsjImvx9fztqQbjTHtjTFxqhl/ecda6znIxwGA3wmKJYRqfxgPk/RezUUfSTUzgDLGnCnp7/v5sE3W2p9fuNLWWrvZGNNBNa8gX7TXzCAA1Jm1doUx5veSnlDNCMd81bw4sPIQnu4BSU8YYx6U9A9r7cOqKdHx+v8leqqk2L1u7yta0r9U82K9KtWU+StrH3tMkpH0Ve2IR76kdyR9vJ+vK8cYc4VqXgjYSTVlfOpvfN79eVk1Ix0/SoqR9KVqfgEAgIATsJutGGMyJX1qre1pjGkmaYW1tmUDPO8rtc/7/uE+FwAAAIJTUIxzWGv3SFprjBkr1cwgGmN+6xXuvzDGJP28RJMxJkU1r3Rf6rOwAAAACHgBWaKNMW9Lmi6pizEmzxhzmWqWTrrMGLNA0hJJB5rr21c3STm1H/e9pH9ZaynRAAAAOKCAHecAAAAAXAnIK9EAAACAS5RoAAAAoJ4Cbom7lJQUm5mZ6ToGAAAAgtycOXO2W2tT9/dYwJXozMxM5eTkuI4BAACAIGeMWX+gxxjnAAAAAOqJEg0AAADUEyUaAAAAqCdKNAAAAFBPlGgAAACgnijRAAAAQD1RogEAAIB6okQDAAAA9USJBgAAAOqJEg0AAADUEyUaAAAAqCdKNAAAAFBPlGgAAACgnijRAAAAQD1RogEAAIB6inAdAACAYGatVVlVtUorq+W1VtZK1qrmfUnREWFqGhWhmMgwGWNcxwVQR5RoAADqyeu12lZUrvU7SrV5d5kKiiqU//PbnnLtKq1UcblHRRUelVR45LUHf05jpCaR4YqNilBCkwilxccovVm00prFKC0+Wi0SYtQ+pak6pMSpSVS4779IAL+JEg0AwAGUV1Vr1bZirdhWpBVb92jt9hKt21GqjTtLVeHx/s+xsVHhSouPVlp8TdmNj4lUXHREzVtMhGKjwhVmjIyRwoxRmPn5c3hVWlmtskqPSiurVVJZrcKySm3bU6E5G3Zp254KVe71uYyRWic2UcfUOGWlxalri3j1b5ekDilNuZINNCJKNAAAkkoqPFq8qVAL8nZrwcZCLduyR+t2lPxyFTk6IkztU5qqY2pTHdM1TW2TY9WueaxaJzZRWrMYxUX75q9Ua632lHm0ubBMawpKtLqgWLn5xVpdUKyZa3eovKqmYCc0iVS/tonq3zZJ/dsmKTszSTGRXLEGfIUSDQAISZt2l2nW2h2atXan5m3YrZXbin4pzG2SmqhHq2Y6pU8rdW0Rry4t4pXZvKnCwxr/Sq8xRgmxkUqIjVS3ls3+5zGv12p1QbHmbtiluet3a97GXZq8skDW1pT+YR2b65iuaRrZJU0ZybGNnh0IZsbaOgxq+ZHs7Gybk5PjOgYAIMBsKSzTlFXbNWP1Ds1cu1ObdpdJkuJjItSvbZL6ZiSqb0aCerdJVEpctOO0h66wrEpzN+zS5BUF+n5FvtbvKJUkdUqL0zHd0jSmT2t1b9XsIM8CQJKMMXOstdn7fYwSDQAIRmWV1ZqxdoemrNyuKasKtCq/WJLUvGmUBrVP/uWta4tmTq4wNwZrrdZuL9F3y/P1w4oCzVizQx6vVdcW8Tqzf2uN6dta6c1iXMcE/BYlGgAQErYWluvb5dv07bJ8/ZS7XRUer6IjwjSofbKO6pSqIzunqEt6fMi+AG9XSaU+XbhZH8zdpPkbdyvMSEdkpWhsdoZO6tlCkeFsHwHsjRINAAhaufnF+mLxFn25ZJsWbSqUJGUkN9GxXdN1TNc0DWqfzAvs9mNNQbEmzNukCfM2KW9XmVo0i9HFwzJ1waC2SoiNdB0P8AuUaABAUFmxtUiTFm3R54u3aOW2mjGN/m0TdVz3dB3XLV2d0uJC9mpzfXm9VpNXFuilqWs1NXe7mkSGa2x2G116RHu1T2nqOh7gFCUaABDwNu4s1cQFmzVx/mat2FYkY6SBmcka3bOFTuzZUi0SmO09XMu27NFLU9dq4vzNqvJ6dUL3Frr5+M7qlB7vOhrgBCUaABCQdpdW6pMFm/XR/M2as36XJCm7XZJO69tKJ/ZsobR4irMv5BeV6/Xp6/Xfn9aptNKjM/u30Q3HdVKbJJbJQ2ihRAMAAoan2qspq7br/Tl5+nrpNlVWe9UlPV5j+rXSqb1bsd5xI9pZUqlnfsjVq9PXS1b63ZC2uuborIBeAhCoD0o0AMDvrdteovGzN+qDuXkqKKpQctMojenbSmcPaKMerRJcxwtpm3eX6fFvV+ndnI1qEhmuq0Z01JVHdeAFmwh6lGgAgF+q8FTryyXbNH7WBk1bvUPhYUZHd0nT2Ow2OrpLmqIiWHLNn6wuKNbDX67Q54u3qkN
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 864x540 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Check error\n",
|
||
|
|
"\n",
|
||
|
|
"def error_check(approx,true,grid) :\n",
|
||
|
|
" error_grid = np.append(0.5*(grid[:-1]+grid[1:]),grid[-1])\n",
|
||
|
|
" return np.sqrt(np.multiply(error_grid,(approx-true)**2))\n",
|
||
|
|
"\n",
|
||
|
|
"true_sol = -(1/np.pi**2)*np.sin(np.pi*grid)\n",
|
||
|
|
"\n",
|
||
|
|
"error = error_check(U,true_sol,grid)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.figure(figsize=(12,7.5))\n",
|
||
|
|
"plt.plot(grid,error)\n",
|
||
|
|
"plt.title(\"Pointwise error\")\n",
|
||
|
|
"plt.ylabel(\"Error\")\n",
|
||
|
|
"plt.xlabel(\"$x$\")\n",
|
||
|
|
"\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 75,
|
||
|
|
"id": "6ae3c077",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"0.000910607929436465"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 75,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# L2 error\n",
|
||
|
|
"sum(error)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"id": "aee347c7",
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": []
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"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.10.4"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|