Initial Commit

.ipynb_checkpoints/validation_new-checkpoint.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "def A_matrix(n):\n",
+    "    A = np.zeros((n-1, n-1))\n",
+    "    for i in range(0, len(A)):\n",
+    "        for j in range(len(A[i])):\n",
+    "            if j == i - 1:\n",
+    "                A[i][j] = 1\n",
+    "            elif j == i:\n",
+    "                A[i][j] = -2\n",
+    "            elif j == i + 1:\n",
+    "                A[i][j] = 1\n",
+    "            else:\n",
+    "                A[i][j] = 0\n",
+    "\n",
+    "    return A\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "def u_exact(x, t, a, b, T):\n",
+    "    L = b-a\n",
+    "    u_e = np.sin(np.pi*x/L)*np.cos(np.pi*t/T)\n",
+    "\n",
+    "    return u_e\n",
+    "\n",
+    "\n",
+    "\n",
+    "def g_function(x, t, a, b, T, c_squared):\n",
+    "    \n",
+    "    L = b-a\n",
+    "\n",
+    "    g = np.sin(np.pi*x/L)*np.cos(np.pi*t/T)*(c_squared*a**2*np.pi**2/(L**2) - b**2*np.pi**2/(T**2))\n",
+    "\n",
+    "    return g\n",
+    "\n",
+    "\n",
+    "def ibvp_step(x, state, n, deltat, deltax, tk, a, b, A, I, c_squared, T):\n",
+    "    \n",
+    "    gt = g_function(x, tk, a, b, T, c_squared)\n",
+    "    \n",
+    "    gt1 = g_function(x, tk+deltat, a, b, T, c_squared)\n",
+    "\n",
+    "    gt_concat = np.concatenate((np.zeros((n-1)), gt), axis=0)\n",
+    "    gt1_concat = np.concatenate((np.zeros((n-1)), gt1), axis=0)\n",
+    "\n",
+    "\n",
+    "    UL = np.zeros((n-1, n-1))\n",
+    "    UR = I\n",
+    "    \n",
+    "    #####PROBLEM LINE\n",
+    "    LL = c_squared/(deltax**2)*A\n",
+    "    \n",
+    "    LR = np.zeros((n-1, n-1))\n",
+    "\n",
+    "    temp1 = np.concatenate((UL,UR), axis=1)\n",
+    "    temp2 = np.concatenate((LL,LR), axis=1)\n",
+    "    f_matrix = np.concatenate((temp1,temp2), axis=0)\n",
+    "\n",
+    "    I_concat = np.identity(2*n-2)\n",
+    "\n",
+    "    LHS = I_concat - 1/2*deltat*f_matrix\n",
+    "    RHS = state + (1/2*deltat*((np.matmul(f_matrix, state)) + gt_concat + gt1_concat))\n",
+    "\n",
+    "    state_updated = np.linalg.solve(LHS, RHS)\n",
+    "    \n",
+    "    return state_updated\n",
+    "    \n",
+    "\n",
+    "def solve_ibvp(state_initial, deltat, deltax, T, a, b, c_squared):\n",
+    "    n = int(abs((b-a)/deltax))\n",
+    "    nt = int(abs(T/deltat))\n",
+    "    x = np.linspace(a+deltax, b-deltax, num=n-1)\n",
+    "    t = np.linspace(0, T, num=nt+1)\n",
+    "    A = A_matrix(n)\n",
+    "    I = np.identity(n-1)\n",
+    "    \n",
+    "    state_hat = [state_initial]\n",
+    "    for h in range(len(t)-1):\n",
+    "        state_updated = ibvp_step(x, state_hat[h], n, deltat, deltax, t[h], a, b, A, I, c_squared, T)\n",
+    "        state_hat.append(state_updated)\n",
+    "    \n",
+    "    state_hat = np.array(state_hat)\n",
+    "\n",
+    "    return state_hat, x, t"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "done\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "a = 0\n",
+    "b = 1\n",
+    "L = b-a\n",
+    "#c_squared = 68.9e9 / 2710\n",
+    "c_squared = 1\n",
+    "deltat = 0.005\n",
+    "deltax = 0.001\n",
+    "T = 1\n",
+    "n = int(abs((b-a)/deltax))\n",
+    "nt = int(abs(T/deltat))\n",
+    "x = np.linspace(a+deltax, b-deltax, num=n-1)\n",
+    "\n",
+    "ui = np.sin(np.pi*x/L)\n",
+    "vi = np.zeros((n-1))\n",
+    "initial_state = np.concatenate((ui, vi), axis=0)\n",
+    "\n",
+    "\n",
+    "state_hat, x, t = solve_ibvp(initial_state, deltat, deltax, T, a, b, c_squared)\n",
+    "\n",
+    "print('done')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "u_true = []\n",
+    "\n",
+    "for d in range(len(t)):\n",
+    "    utk = u_exact(x, t[d], a, b, T)\n",
+    "    u_true.append(utk)\n",
+    "\n",
+    "u_true = np.array(u_true)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.2825066996978956e+16\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyu0lEQVR4nO3deVxV1f7/8deHSQQRQQFRVFAR5xGHykxF09SyvM2TjdZNy6ZbdvveW90pb3lvaTbXLZsHszQ1UzFNs0HMeWJQVBSZVBRExvX7g1M/BxDwnMM+h/N5Ph7ncfZhL896b8SPi332XkuMMSillGr4vKwOoJRSqn5owVdKKQ+hBV8ppTyEFnyllPIQWvCVUspD+Fgd4FxatGhhoqOjrY6hlFJuY/369bnGmLCq9rl0wY+OjiYpKcnqGEop5TZEZG91+/SUjlJKeQgt+Eop5SG04CullIfQgq+UUh5CC75SSnkIhxR8ERktIrtEJFVEplWxX0Rklm3/ZhHp64h+lVJK1Z7dBV9EvIGXgcuArsANItL1jGaXAbG2xyTgVXv7VUopVTeOuA5/AJBqjNkNICKfAOOB7ae0GQ+8ZyrnYv5JRJqJSKQxJtMB/Z+losLAnLF4Iafv6HYlDLgbSk7Ah9ec/Qd73wh9boLCPPjs1rP3978Duv8B8jNg3j1n779wCsRdBrkp8PWDZ+8f8ih0GAaZm2HJE2fvT/grtB0I+36GxL+dvX/0sxDZE9K+g+9nnL3/8hehRSzs+gbWzj57/4TXITgKtn4B6/539v5r34PA5rDhQ9j40dn7b/oc/ALglzdh21dn7799UeXzD7Mg+dvT9/n6w81fVG6veg52rzp9f0AIXPdB5fbyp2H/utP3N20Ff3izcvubaXBoy+n7m3eAK2ZVbi94APLSTt/fsgdcNr1y+4u74djB0/e36Q8jnq7c/vRmOHHk9P3tL4FLHqvc/uAPUHry9P2dRsFFD1RuvzOWs+jPnv7sQa1+9ioKDyO3L0LkjPrlAI44pdMa2H/K6wzb1+raBgARmSQiSSKSlJOTU+cwx06WctnM1RzKP1lzY6WUcjEHjhZxxewfOFla7vD3FnsXQBGRa4BRxpi7bK9vAQYYY+4/pc0i4FljzBrb60TgMWPM+nO9d3x8vDmfO21vfutndh46zurHhtHYz7vOf14ppaxw7GQpF01fwYUdmvP6LfHn9R4ist4YU+UfdsQIPwNoc8rrKODgebRxmAcSYsktKOajX/Y5qwullHK4OT+kc/xkGfcPj3XK+zui4K8DYkUkRkT8gOuBBWe0WQDcartaZxCQ76zz9wADYkK5oH1zXluV5pRfi5RSytEKist4a80eRnQJp3vrYKf0YXfBN8aUAVOAb4EdwGfGmG0icq+I3GtrthjYDaQCbwL32dtvTR5IiCXneDGf6ChfKeUG3vsxnfyiUqeN7sFBs2UaYxZTWdRP/dprp2wbYLIj+qqtCzo0Z0BMKK+uSuP6AW3x99Vz+Uop11RYXMZbq/cwNC6MXm2aOa2fBn2n7dSEWLKOFfN50v6aGyullEU+/HkvhwtLnDq6hwZe8C/s0Jz4diG8sjKN4jI9l6+Ucj1FJeW88f1uLo5tQb92IU7tq0EXfBHhgYRYMvNPMnd9htVxlFLqLB//so/cghIeSHDu6B4aeMEHuDi2BX3aNuOV79IoKauwOo5SSv3uZGk5r61K44L2zekfHer0/hp8wRcRpibEcuBoEfN+1VG+Usp1fJa0n+zjxfUyugcPKPgAl3QKo1dUMLO/S6W0XEf5SinrFZeV8+rKNPpHhzCovfNH9+AhBV9EmDoilowjRXy54YDVcZRSirnrM8jMP8kDCbFOmSitKh5R8AGGxYXTo3UwL3+XSpmO8pVSFiopq+CV79Lo07YZgzu2qLd+Pabg/3bFzt68E8zf6LRpfJRSqkZfbsjgwNGieh3dgwcVfIARXcLpGtmU2TrKV0pZpLS8gtnfpdIzKpihncLqtW+PKvi/jfL35BaycLPT5m5TSqlqzd94kP2Hi3hgeP2O7sHDCj7ApV0j6NwyiFkrUiivsG8tAKWUqouy8gpe/i6VrpFNSegSXu/9e1zB9/KqHOXvzilk0RYd5Sul6s/Xmw+yJ7ew3s/d/8bjCj7A6G4t6RTRhJcSUyrXv1VKKScrrzC8lJhK55ZBXNo1wpIMHlnwvbyE+4fHkpJdwDdbD1kdRynlAb7edJDduYU8OCIWL6/6H92DhxZ8gDE9IukQFsgsHeUrpZysvMIwa0WKbXTf0rIcdhV8EQkVkWUikmJ7rnJuTxH5n4hki8hWe/pzJG/bufxdWcdZul1H+Uop51m4+SC7cyrP3Vs1ugf7R/jTgERjTCyQaHtdlXeB0Xb25XDjeraifYtAZiam6ihfKeUU5RWGl1akEhcRxOhu1o3uwf6CPx6YY9ueA1xZVSNjzPfAYTv7cjhvL2HK8I7syDzG8h1ZVsdRSjVAi7ZkkppdwP0JHS0d3YP9BT/CGJMJYHu2+8JSEZkkIkkikpSTk2Pv29Xoil6tiG4ewMzEFCqX3lVKKceoqDC8lJhCbHgTxnSPtDpOzQVfRJaLyNYqHuOdEcgY84YxJt4YEx8W5vzbjn28vZg8rCPbDh5jxc5sp/enlPIci7dmkpJdYPm5+9/41NTAGDOiun0ikiUikcaYTBGJBNyyYl7ZpzWzVqQwMzGF4Z3DLbkhQinVsFRUGGYlptAxvAljelg/ugf7T+ksACbaticC8+18P0v4ensxZVhHNmfks3KX808jKaUaviXbDpGcVcD9wzvi7QKje7C/4E8HRopICjDS9hoRaSUii39rJCIfAz8CcSKSISJ32tmvw03oG0XrZo31XL5Sym4VFYaZy1PoEBbIuJ6trI7zO7sKvjEmzxiTYIyJtT0ftn39oDFmzCntbjDGRBpjfI0xUcaYt+0N7mi+tnP5G/cf5fuUXKvjKKXc2LfbDrEr6zj3D491mdE9ePCdtlW5ul8UrYL9mbk8WUf5SqnzUlFhmJmYQvsWgVzey3VG96AF/zR+Pl78cVhHft13lB9S86yOo5RyQ0u3Z7Hz0HHuT3Cdc/e/0YJ/hmvjo2jZ1J+ZiTrKV0rVzW+j+5gWgVzuQufuf6MF/wyNfLz549AOrEs/wo+7dZSvlKq9ZTuy2JF5jCnDOuLj7Xrl1fUSuYDr+rchPKgRsxJTrI6ilHITxlRedx/dPIDxvV1vdA9a8Kvk7+vNvZd04Kfdh/lZR/lKqVpYviObbQePMdlFR/egBb9aNw5sS4smjZi1Qkf5SqlzM8YwMzGZds0DuKpPa6vjVEsLfjUqR/nt+SE1j3XpLjfRp1LKhSTuyGbrAdce3YMW/HO6aWA7WjRpxAvLkq2OopRyUZWj+xTahDZ26dE9aME/p8Z+3tw3tANr0/L4MU3P5Sulzvbdrmy2HMhnyrCO+Lrw6B604NfoxoFtCQ9qxAt6961S6gzGGF5cnkJUSGMm9I2yOk6NtODXwN/Xm8nDOvLLnsOs1VG+UuoUiTuy2ZyRzwPDY11+dA9a8Gvluv5tiAz254VlOspXSlUyxvDfZbYrc/q69rn732jBrwV/X2/uG9aRpL1HWK0zaSqlgG+3ZbE985jbjO5BC36tXRtfOV/+f3WUr5THq6gwvLg8mfYtAl32rtqqaMGvpUY+3kwZXjlfvq6KpZRn+2brIXYeOs7UEbEufd39mexKKiKhIrJMRFJszyFVtGkjIt+JyA4R2SYiU+3p00pX94siKqSxXrGjlAcrt43uO4Y3canVrGrD3v+apgGJxphYINH2+kxlwCPGmC7AIGCyiHS1s19L+Hp78cDwWDZn5JO4wy3Xa1dK2Wnh5oOkZBfw4AjXWs2qNuwt+OOBObbtOcCVZzYwxmQaY361bR8HdgDu8ZF2Fa7q25p2zQP0XL5SHqisvIKZy1OIiwhiTPdIq+PUmb0FP8IYkwmVhR0IP1djEYkG+gA/29mvZX4b5W/PPMa327KsjqOUqkcLNh1kd24hD42MxcvNRvdQi4IvIstFZGsVj/F16UhEmgBfAA8aY46do90kEUkSkaScHNf8cHR871bEtAjkxeXJVFToKF8pT1BWXsHMxBS6Rjbl0q4trY5zXmos+MaYEcaY7lU85gNZIhIJYHuu8sS2iPhSWew/NMbMq6G/N4wx8caY+LCwsLofUT3w8fZiakIsOw8dZ8m2Q1bHUUrVg3m/HmBv3gkeGtnJLUf3YP8pnQXARNv2RGD+mQ1ERIC3gR3GmP/a2Z/LuLxXKzqEBfLCsmTKdZSvVINWUlbBrBUp9GgdzIgu5zxz7dLsLfjTgZEikgKMtL1GRFqJyGJbm4uAW4DhIrLR9hhjZ7+W8/YSpo7oREp2AYu2ZFodRynlRHPXZ5BxpIiHR3aicgzrnnzs+cPGmDwgoYqvHwTG2LbXAO77HTqHsT0imb0ihZnLkxnbI9LtLtFSStWsuKyc2StS6N2mGUPjXPM0c225zy1iLsjbS3hwRCfScgr5etNBq+MopZzgs3X7OZh/0u1H96AF326ju7Wkc8sgZiamUFZeYXUcpZQDnSwtZ/Z3qfSPDuHi2BZWx7GbFnw7edlG+XtyC5m34YDVcZRSDvTxL/vIOlbMQw1gdA9a8B1iVLcIerQOZubyFIrLyq2Oo5RygBMlZbz8XRoDY0K5sIP7j+5BC75DiAiPjorjwNEiPl233+o4SikHmLN2L7kFxTw6Ks7qKA6jBd9BhsS2YEBMKC+tSKWoREf5Srmz/KJSXluVxtC4MPpHh1odx2G04DuIiPCnUXHkHC9mzo/pVsdRStnh7dW7yS8q5dFLG87oHrTgO1T/6FCGxoXx6so0jp0stTqOUuo85BUU8/aaPYztEUn31sFWx3EoLfgO9uilceQXlfLW6j1WR1FKnYdXV6ZRVFrOQyM7WR3F4bTgO1j31sGM6dGSt1fvJq+g2Oo4Sqk6yMwv4r2f9jKhbxQdw5tYHcfhtOA7wcMjO1FUWs5rq9KsjqKUqoNZiakYY5iaEGt1FKfQgu8EHcODuKpPFHN+3Muh/JNWx1FK1UJ6biGfJ+3nhgFtaRMaYHUcp9CC7yQPjojFGMNLK1KsjqKUqoUXlyfj4y1MGdbR6ihOowXfSdqEBnB9/7Z8um4/+/JOWB1HKXUOuw4dZ/6mg9x2YQzhTf2tjuM0WvCdaMrwjnh7CS8uT7Y6ilLqHP6zdBdN/Hy495L2VkdxKi34ThTR1J/bLozmy40HSM46bnUcpVQVNu4/ytLtWdw9pD3NAvysjuNUWvCd7N5LOhDo58N/l+ooXylX9J+luwgN9OOOwTFWR3E6uwq+iISKyDIRSbE9h1TRxl9EfhGRTSKyTUSesadPdxMS6Medg2NYsu0QmzOOWh1HKXWKH9PyWJ2Sy31DO9CkkV0LALoFe0f404BEY0wskGh7faZiYLgxphfQGxgtIoPs7Net3HVxDM0CfHn+211WR1FK2RhjeP7bnbRs6s/Ng9pZHade2FvwxwNzbNtzgCvPbGAqFdhe+toexs5+3UqQvy9ThnVkdUoua1JyrY6jlAKWbs/i131HmToiFn9fb6vj1At7C36EMSYTwPYcXlUjEfEWkY1ANrDMGPNzdW8oIpNEJElEknJycuyM5zpuHtSO1s0a8+8lO6mo8Kj/75RyOWXlFTy3ZCcdwgK5pl+U1XHqTY0FX0SWi8jWKh7ja9uJMabcGNMbiAIGiEj3c7R9wxgTb4yJDwtz7xXiT+Xv683DIzux5UA+i7ZkWh1HKY82d30GaTmF/GlUZ3y8PefalRqP1BgzwhjTvYrHfCBLRCIBbM/ZNbzXUWAlMNr+6O7nyj6t6dwyiBlLd1FSpgueK2WFopJyXlieTJ+2zRjVLcLqOPXK3v/aFgATbdsTgflnNhCRMBFpZttuDIwAdtrZr1vy9hIeH92ZvXkn+GTdPqvjKOWR3lm7h6xjxUwb3blBLExeF/YW/OnASBFJAUbaXiMirURksa1NJPCdiGwG1lF5Dn+hnf26raFxYQyMCWVWYgoFxWVWx1HKoxwpLOHVlWkkdA5nYPvmVsepd3YVfGNMnjEmwRgTa3s+bPv6QWPMGNv2ZmNMH2NMT9upoL85Iri7EhGmXdaZ3IIS3lq92+o4SnmUV1amUlBcxmOjO1sdxRKe82mFC+nTNoQxPVry5ve7yTmui6QoVR8yjpxgztq9/KFvFHEtg6yOYwkt+BZ59NI4TpZVMFunT1aqXrywLAWEBrl0YW1pwbdI+7AmXN+/DR/+vI+9eYVWx1GqQdt56BjzNmRw24XRtG7W2Oo4ltGCb6GpCbH4ensxQydWU8qpnluyiyaNfLhvaAero1hKC76Fwpv6c9fFMXy96SBbMvKtjqNUg/TT7jxW7MzmvqEdG/z0xzXRgm+xSUPaExLgy7+XeOStCUo5lTGG6d9UTpB2+0XRVsexnBZ8iwX5+3L/8FjWpOayKrnhzB2klCtYvOUQG/cf5aGRnjNB2rlowXcBNw1qS9vQAJ5dvINynVhNKYcoLitn+pIddG4ZxNX92lgdxyVowXcBjXy8mXZZZ3YeOs7nSfutjqNUg/De2r3sP1zEn8d0wdvLs6ZQqI4WfBdxWfeWxLcL4T/LknXKBaXsdKSwhJdWpHBJpzCGdGo4s+7aSwu+ixARnhzbhZzjxby+Ks3qOEq5tZm2uaqeHNvF6iguRQu+C+nTNoQrerXizdW7ycwvsjqOUm5pd04BH/y0l+v6t6VThGdOoVAdLfgu5rHRcVQYdP1bpc7T9G920sjHi4c9eAqF6mjBdzFRIQHccVEM8349oDdjKVVHP+3OY+n2LP44tANhQY2sjuNytOC7oPuGdSA00I9/LNqOMXqZplK1UVFh+OeiHUQG+3Pn4PZWx3FJWvBdUFN/Xx4a2Ymf9xxm2fYsq+Mo5RbmbzrAlgP5/GlUHI399CarqthV8EUkVESWiUiK7TnkHG29RWSDiHjsald1cUP/NnQMb8L0b3ZSWq7r3yp1LidLy3l+yS56tA7myt6trY7jsuwd4U8DEo0xsUCi7XV1pgI77OzPY/h4e/HnMZ3ZnVvIhz/ttTqOUi7t7TV7OJh/kifHdsFLb7Kqlr0Ffzwwx7Y9B7iyqkYiEgWMBd6ysz+PMiwunIs6NufFxBSOniixOo5SLin72Ele+S6VkV0jGOSB69TWhb0FP8IYkwlgew6vpt2LwGNAjecmRGSSiCSJSFJOjmdPJiYi/N/YrhwrKuXF5boyllJVmb5kJ6XlhifH6E1WNamx4IvIchHZWsVjfG06EJFxQLYxZn1t2htj3jDGxBtj4sPC9JboLpFNuXFgW97/aS+7Dh23Oo5SLmXDviPM+/UAdwyOIbpFoNVxXF6NBd8YM8IY072Kx3wgS0QiAWzP2VW8xUXAFSKSDnwCDBeRDxx4DA3eIyPjaNLIh78t3KaXaSplU1FhePrr7YQHNWLK8I5Wx3EL9p7SWQBMtG1PBOaf2cAY84QxJsoYEw1cD6wwxtxsZ78eJSTQj4dHduKH1Dy+3aaXaSoF8OWGA2zaf5THR3emSSMfq+O4BXsL/nRgpIikACNtrxGRViKy2N5w6v+7aWBb4iKC+Ofi7ZwsLbc6jlKWKiguY/qSnfRq04yr+uhlmLVlV8E3xuQZYxKMMbG258O2rx80xoypov1KY8w4e/r0VD7eXjx1eVf2Hy7irdW7rY6jlKVmr0gl53gxT1/eVS/DrAO909aNXNixBaO7teTl79J0Nk3lsdJzC/nfmj1M6NuaPm2rvddTVUELvpt5cmwXym0LMyvlif6xaAe+3sK00Z2tjuJ2tOC7mTahAdwzpD3zNx5kXfphq+MoVa++T85h+Y4spgyPJbypv9Vx3I4WfDf0x6EdiAz25+kF23TRc+UxSsoq+NvC7bRrHsAdg6OtjuOWtOC7oQA/H6Zd1pltB4/x6Tpd9Fx5hnd+2ENqdgF/HdeVRj46G+b50ILvpq7o1YoBMaE89+1ODhfqPDuqYcvML2JmYgojukSQ0CXC6jhuSwu+mxIR/j6+OwUny/i3foCrGri/L9xOeYXhqcu7Wh3FrWnBd2NxLYO4c3AMnybtZ/1e/QBXNUzfJ+eweMshpgzrSJvQAKvjuDUt+G7ugYRYIoP9efLLrZTpQimqgSkuK+epBduIaRHIpEt02UJ7acF3c4GNfPjruK7sPHSc937UhVJUw/LGqt3syS3kmSu66Qe1DqAFvwEY3b0ll3QK47/Lksk6dtLqOEo5xP7DJ5j9XSpjerRkSCedKt0RtOA3ACLCM1d0o6S8gn8s0lUkVcPwzNfb8PYS/jJOP6h1FC34DUR0i0DuG9qBrzcd5IfUXKvjKGWX5duzWL4jm6kJsUQGN7Y6ToOhBb8BufeSDrRrHsBf5m+luEynUFbu6URJGU8t2EZseBPuGBxjdZwGRQt+A+Lv683TV3Rjd04hb36vUygr9/TCsmQOHC3iXxN64OutJcqR9LvZwAyLC2dsj0hmrUhld06B1XGUqpOtB/J5e80ebhzYlv7RoVbHaXDsKvgiEioiy0QkxfZc5eTUIpIuIltEZKOIJNnTp6rZU5d3pZGPF0/M20KFTq6m3ERZeQXT5m2meZNGPK5THzuFvSP8aUCiMSYWSLS9rs4wY0xvY0y8nX2qGoQ39efJMV34ec9hPkvSydWUe3h3bTpbDxzjmSu6EdzY1+o4DZK9BX88MMe2PQe40s73Uw5yXf82DIwJ5V+Ld5B9XK/NV65t/+ET/GdpMiO6hHNZ95ZWx2mw7C34EcaYTADbc3g17QywVETWi8ikc72hiEwSkSQRScrJybEznucSEZ6d0IOTZRU8s2C71XGUqpYxhr/M34qXwN/Gd0dE16h1lhoLvogsF5GtVTzG16Gfi4wxfYHLgMkiMqS6hsaYN4wx8caY+LAwvbvOHu3DmvDA8I4s2pLJ8u1ZVsdRqkoLN2eyclcOj46Ko1UzvebemXxqamCMGVHdPhHJEpFIY0ymiEQC2dW8x0Hbc7aIfAkMAL4/z8yqDiYN6cDCzZn8Zf5WBrYPJchfz40q13H0RAnPfL2NXlHB3HpBtNVxGjx7T+ksACbaticC889sICKBIhL02zZwKbDVzn5VLfn5ePHshB4cOnaS57/dZXUcpU7zt6+3c/REKc9O6Im3l57KcTZ7C/50YKSIpAAjba8RkVYistjWJgJYIyKbgF+ARcaYJXb2q+qgT9sQJl4Qzfs/7eWXPTpvvnINiTuymLfhAPcN60jXVk2tjuMRxBjXvU47Pj7eJCXpZfuOUFhcxuiZ3+MlwjdTLybAr8azeUo5Tf6JUi59cRUhAX4smDIYPx+9B9RRRGR9dZe/63fZQwQ28uH5q3uxN+8Ezy3RUzvKWn9ftJ3cghJmXNNLi3090u+0BxnUvjm3XRjNu2vT+Wl3ntVxlIf6bmc2c9dncN/QDnRvHWx1HI+iBd/DPDY6jnbNA3hs7mYKi8usjqM8TH5RKU/M20KniCZMGd7R6jgeRwu+hwnwqzy1s//ICf69ZKfVcZSH+eei7eQUFDPjml66ZKEFtOB7oAExodx+YQzv/biXtWm6WIqqH9/tyuazpAzuGdKenlHNrI7jkbTge6g/jYojpkWgntpR9eJwYQmPzd1Mp4gmPJAQa3Ucj6UF30M19vPm+at7cuBoka6Dq5zKGMOf523h6IkSXryuD/6+eirHKlrwPVh8dCiThrTn41/2sXTbIavjqAbqi18PsGTbIR65NE5vsLKYFnwP98jIOLpGNmXavC06jbJyuP2HT/D0gm0MiAnl7ovbWx3H42nB93B+Pl7MuqE3hcVl/OnzzbjyndfKvZRXGB7+bCMA/722l86V4wK04Cs6hgfx5NgurErOYc7adKvjqAbije93sy79CM9c0Y2okACr4yi04CubWwa1Y1hcGP/6ZifJWcetjqPc3LaD+fx32S7G9GjJhL6trY6jbLTgK6Byhaznru5FUCMfHvh4A8Vl5VZHUm6qsLiM+z/eQEiAH/+8soeuYOVCtOCr34UFNeK5q3uy89BxntcJ1tR5emrBNvbkFvLi9b0JCfSzOo46hRZ8dZqELhHcMqgdb63ZQ+IOXRZR1c2XGzKYuz6D+4d15MIOLayOo86gBV+d5cmxXegS2ZRHPt/EwaNFVsdRbmJ3TgFPfrmVAdGhejeti7Kr4ItIqIgsE5EU23NINe2aichcEdkpIjtE5AJ7+lXO5e/rzSs39aW0rIL7P95AaXmF1ZGUiztZWs6UjzbQyMeLmTf0xsdbx5KuyN6/lWlAojEmFki0va7KTGCJMaYz0AvQe/ldXEyLQP41oQfr9x5hxlI9n6/Obfo3O9meeYwZ1/QiMrix1XFUNewt+OOBObbtOcCVZzYQkabAEOBtAGNMiTHmqJ39qnowvndrbhjQltdX7ea7ndlWx1EuasnWQ7y7Np07LoohoUuE1XHUOdhb8COMMZkAtufwKtq0B3KAd0Rkg4i8JSKB1b2hiEwSkSQRScrJybEznrLXU5d3pXPLIB7+bKOez1dn2Z1TwKOfb6JXVDCPXxZndRxVgxoLvogsF5GtVTzG17IPH6Av8Koxpg9QSPWnfjDGvGGMiTfGxIeFhdWyC+Us/r7evHxTX4rLKpj80a96fb76XWFxGfd+sB5fb+GVm/vpgiZuoMaCb4wZYYzpXsVjPpAlIpEAtueqfu/PADKMMT/bXs+l8j8A5SY6hDXh+at7sWHfUf729Xar4ygXYIzh8S82k5pdwEs39KV1Mz1v7w7sPaWzAJho254IzD+zgTHmELBfRH77fS8B0KrhZsb2jOSeS9rz4c/7+GzdfqvjKIu980M6Czdn8silcQyO1evt3YW9BX86MFJEUoCRtteISCsRWXxKu/uBD0VkM9Ab+Jed/SoL/OnSOAZ3bMH/fbWVTfuPWh1HWeSXPYf51+IdjOwawR8v6WB1HFUH4srT4cbHx5ukpCSrY6hTHCksYdxLa6gwhq/vH0yLJo2sjqTq0aH8k1w+ew2Bft4suH8wTf19rY6kziAi640x8VXt07sjVJ2EBPrx+i39OFxYwuQPf9WbsjxIUUk5d7+XxIniMl6/JV6LvRvSgq/qrHvrYKb/oQc/7znM0wu26aIpHsAYw6NzN7H1YD6zbuhDXMsgqyOp8+BjdQDlnq7qE8WuQwW8tiqNDmFNuGNwjNWRlBPNSkxl0eZMnriss95c5ca04Kvz9tioOPbkFvCPRduJbhHA8M5aCBqiRZszeWF5MhP6tmbSEF2X1p3pKR113ry8hBeu603XVk25/6MN7Dx0zOpIysG2Hsjnkc830rdtM56doIuZuDst+MouAX4+vHVrf5r4+3Dnu0nkHC+2OpJykP2HT3D7u+sIDfDj9Vvi9U7aBkALvrJby2B/3p7Yn8OFJdw1Zx0nSsqsjqTsdPRECbe98wvFpeW8e8cAwoL08tuGQAu+cojurYOZdUMfthzI5z69XNOtnSytvPxy/+Ei3rg1nk4RekVOQ6EFXznMyK4R/OuqHqzclcPjczdTUaGXa7qb8grDQ59uZF36Ef57XS8GtW9udSTlQHqVjnKo6we0Jed4Mf9ZlkxYUCOeGNPF6kiqlowx/H3hdr7Zeoj/G9uFcT1bWR1JOZgWfOVwU4Z3JKegmNe/301YUCPuulgv5XMHLyxL5t216dw5OEb/zhootyv4paWlZGRkcPLkSaujWMrf35+oqCh8fV3v9nYR4anLu5FbUMw/Fu0gyN+H6/q3tTqWOofXV6Uxa0Uq18ZH8aT+VtZguV3Bz8jIICgoiOjoaI+9JtgYQ15eHhkZGcTEuOYdrt62a/QLi9czbd4WfL29mNA3yupYqgrv/7SXZ7/ZybiekTw7oSdeXp7578oTuN2HtidPnqR58+YeW+yhcgTdvHlzl/8tp5GPN6/f0o8L2jfn0c83sXDzQasjqTPM+zWDv3y1lYTO4bxwXW+8tdg3aG5X8AGPLva/cZfvgb+vN29NjKdfuxCmfrKRb7cdsjqSspm/8QCPfr6JCzs05+Wb+uLr7ZblQNWB/g0rpwvw8+F/t/WnR+tgpnz0K8u2Z1kdyeN9sT6Dhz7dSP/oUN68NR5/X72L1hPYVfBFJFRElolIiu05pIo2cSKy8ZTHMRF50J5+lfsJ8vdlzh0D6NoqmD9+sJ6vN+npHat8um4fj87dxIUdWvDu7QMIbOR2H+Wp82TvCH8akGiMiQUSba9PY4zZZYzpbYzpDfQDTgBf2tmvWxg6dCh1WbErPT2djz76yImJrBXc2JcP7hxA33YhTP1kA58l6dq49e2Dn/by+BdbGBIbxlsT42nspyN7T2Lvf+3jgaG27TnASuDxc7RPANKMMXvt7BeAZ77exvaDjp2hsWurpjx1eTeHvmdt/Vbwb7zxRkv6rw9B/r7MuX0Ak95P4rG5mzlZWs6tF0RbHavBM8bw6qo0nluyixFdwnn5pr46GZoHsneEH2GMyQSwPYfX0P564ONzNRCRSSKSJCJJOTk5dsZzvPT0dLp37/776xkzZvD0009X2/7zzz9nwIABdOrUidWrV//+HhdffDF9+/alb9++rF27FoBp06axevVqevfuzQsvvODU47BSY7/KD3JHdo3gr/O38eLyZF01y4kqKgzPfL2d55bsYnzvVrxyUz8t9h6qxhG+iCwHWlax68m6dCQifsAVwBPnameMeQN4AyoXMT9XW6tG4nVRVlbGL7/8wuLFi3nmmWdYvnw54eHhLFu2DH9/f1JSUrjhhhtISkpi+vTpzJgxg4ULF1od2+ka+Xjzyk19eWLeFl5cnsLBo0X886oeeqWIgxWXlfPIZ5tYuDmTOwfH8OSYLnqdvQerseAbY0ZUt09EskQk0hiTKSKRQPY53uoy4FdjjEddojFhwgQA+vXrR3p6OlB5t/CUKVPYuHEj3t7eJCcnW5jQOr7eXjx/dU9aN2vMzMQUDh0r5pWb+tJEP0R0iPyiUu77cD0/pObxxGWdmTSkvdtczqucw97h1AJgom17IjD/HG1voIbTOe7Ax8eHior/P/VvTTc/NWpUOY+4t7c3ZWWV88S/8MILREREsGnTJpKSkigpKXFeYBcnIjw0shP//kMPfkjN5brXf+RQvmvfUOYOducUcNUrP/Dz7sPMuKYX91zSQYu9srvgTwdGikgKMNL2GhFpJSKLf2skIgG2/fPs7M9yERERZGdnk5eXR3Fx8XmdfsnPzycyMhIvLy/ef/99ysvLAQgKCuL48eOOjuwWruvflrcmxpOeW8i4l9aQlH7Y6khua01KLle+/ANHCkv44K6BXN1Pp7RQlewq+MaYPGNMgjEm1vZ82Pb1g8aYMae0O2GMaW6Mybc3sNV8fX3561//ysCBAxk3bhydO3eu83vcd999zJkzh0GDBpGcnExgYCAAPXv2xMfHh169ejXoD22rMywunK8mX0STRt7c8OZPfPTzPqsjuRVjDO/9mM7Ed36hZbA/C6YM1vns1WnEla+OiI+PN2dex75jxw66dNHZ/KDhfi/yT5TywCcbWJWcww0D2vL0FV31qpIaFBaX8ecvtzB/40ESOofz4vW9CfJ3vZlUlfOJyHpjTHxV+/TTMeVyggN8+d9t/ZmxdBevrkxjc8ZRXrqhD+3DmlgdzSXtOnSc+z5cz57cQh4Z2YnJwzrqlTiqSnoNnANMnjyZ3r17n/Z45513rI7l1ry9hMdHd+bNW+M5cLSIcS+t4Yv1GVbHcinGGD5dt4/xL68hv6iMD+4ayP0JsVrsVbV0hO8AL7/8stURGqyRXSP4ZurFTP1kI498vonvU3J45opuNAvwszqapbKPn+SJL7aQuDObCzs058XrexMe5G91LOXitOArlxcZ3JiP7x7E7BWpvLQihbVpefzjyu6M6lbV/YAN36LNmfzfV1s4UVLOX8Z15fYLo3VUr2pFT+kot+DtJUwdEcv8KRcR1qQR97y/nvs/3kBuQbHV0erNgaNFTHovickf/Urb0AAWPXAxdw6O0WKvak1H+MqtdGsVzPwpF/HqyjReWpHCyl3ZPDyyEzcPatdgp2UoLa/g7TV7mLk8BYDHR3fm7otj8Gmgx6ucR39inOjdd9/l4EGd993RfL29eCAhlm+mXkzvNs145uvtjJ21mh9Sc62O5lDGGJZtz+KymauZ/s1OBse2YNnDQ/jj0A5a7NV50Z8aJzpXwf/t7lp1/jqGB/HeHQN445Z+FJWWc9NbP3PL2z+zYd8Rq6PZLSn9MNe89iN3v5dERYXhrVvjefPWeKJCAqyOptyY+5/SeWfs2V/rdiUMuBtKTsCH15y9v/eN0OcmKMyDz249fd/ti87ZXXp6OuPGjWPr1q1A5fTIBQUFZ02RPHfuXJKSkrjpppto3LgxP/74I126dOGOO+5g6dKlTJkyhddee40ZM2YQHx9Pbm4u8fHxpKenU15ezrRp01i5ciXFxcVMnjyZe+65pw7fFM8hIlzarSVDOoXx3o/pvLZqN1e9spaEzuE8OKITPaKCrY5Ya8YYftlzmFdXpbFyVw5hQY3451XduTa+TYM9XaXql/sXfBd19dVXM3v27N8L+m/8/f1Zs2YNAK+99lqVf/btt98mODiYdevWUVxczEUXXcSll15KTExMvWR3R/6+3kwa0oEbB7Zjztp0Xl+VxuWz1zAgJpQ7B8cwoksE3i764WZZeQWJO7N5fVUav+47SvNAP/40Ko7bL4omwE//iSrHcf+fpnONyP0Czr0/sHmNI3pHu+6662pss3TpUjZv3szcuXOBysnWUlJStODXQpNGPkwe1pFbLmjHp7/s59216dzz/nrahgZwbXwUV/WNonWzxlbHBCDjyAk+W7efT5P2k3WsmKiQxvxtfDeu6ddGlx5UTuH+Bb+e1XV65DP9NlHame916vsYY3jppZcYNWqUnWk9V1N/X+4e0p7bL4rm221ZvPdjOjOWJvOfZclc2KE543q2IqFzOOFN6/dmpUP5J1myNZPFWw+xzjYj6CWdwnjmiraM6BKuH8Yqp9KCX0enTo/cpEkTFi5cyOjRo6tsW9N0x9HR0axfv54BAwb8PpoHGDVqFK+++irDhw/H19eX5ORkWrdufdp/Fqp2fLy9GNszkrE9I9l/+ARf/JrBvF8P8MS8LQD0igpmaFw4/aND6dO2GYEOXnzlREkZSelHWJuWx9q0XDZnVE4YGxcRxIMJnfhDv9b6QayqN1rw6+jU6ZFjYmLOOT3ybbfdxr333vv7h7ZnevTRR7n22mt5//33GT58+O9fv+uuu0hPT6dv374YYwgLC+Orr75yxuF4lDahATw4ohNTE2LZlXWcxB3ZLNuexawVKRhTeXNX55ZBdG7ZlNiIJsSGNyEqJICwoEY0a+xb7Q1Oxhjyi0rJOFLEwaNF7M07wbaD+WzPPEZaTiHlFQZfb6FPmxAevbQTo7tH0jFcJ4JT9U+nR3Zj+r1wjGMnS9mw7yhJ6YfZsO8oyVnHyT5++h28Pl5CswBffL298PEWfLy8OFlaTmFxGSdKyimrOP3fUcum/nRr1ZSurZrSPzqU+OgQ/QBW1QunTY8sIqHAp0A0kA5ca4w56yJoEXkIuAswwBbgdmOMrmOnXEJTf18u6RTGJZ3Cfv9aflEpqdkFZOYXkXu8mJyCYo6cKKWsvILSckNpeQWNfb0JbORDgJ83oYF+RIU0plWzxkSFBBAa6NmTuynXZO+QYxqQaIyZLiLTbK8fP7WBiLQGHgC6GmOKROQz4HrgXTv7dhmTJ0/mhx9+OO1rU6dO5fbbb7cokbJXcGNf+rULAUKsjqKUw9hb8McDQ23bc4CVnFHwT+mnsYiUAgFAg5pvQKdHVkq5A3uvAYswxmQC2J7Dz2xgjDkAzAD2AZlAvjFmqT2duvLnDvVFvwdKqbqqseCLyHIR2VrFY3xtOhCRECp/E4gBWgGBInLzOdpPEpEkEUnKyck5a7+/vz95eXkeXfCMMeTl5eHvrwteKKVqr8ZTOsaYEdXtE5EsEYk0xmSKSCSQXUWzEcAeY0yO7c/MAy4EPqimvzeAN6DyKp0z90dFRZGRkUFV/xl4En9/f6KioqyOoZRyI/aew18ATASm257nV9FmHzBIRAKAIiABSKqiXa34+vrqFANKKXUe7D2HPx0YKSIpwEjba0SklYgsBjDG/AzMBX6l8pJML2wjeKWUUvXH7W68UkopVb1z3XilMzUppZSHcOkRvojkAHvP44+2ABrWenc102P2DHrMnsGeY25njAmraodLF/zzJSJJ1f1K01DpMXsGPWbP4Kxj1lM6SinlIbTgK6WUh2ioBd8TL/vUY/YMesyewSnH3CDP4SullDpbQx3hK6WUOoMWfKWU8hBuXfBFZLSI7BKRVNsCLGfuFxGZZdu/WUT6WpHTkWpxzDfZjnWziKwVkV5W5HSkmo75lHb9RaRcRK6uz3zOUJtjFpGhIrJRRLaJyKr6zuhotfjZDhaRr0Vkk+2Y3X6FIRH5n4hki8jWavY7toYZY9zyAXgDaUB7wA/YROWqWqe2GQN8AwgwCPjZ6tz1cMwXAiG27cs84ZhPabcCWAxcbXXuevh7bgZsB9raXodbnbsejvnPwL9t22HAYcDP6ux2HvcQoC+wtZr9Dq1h7jzCHwCkGmN2G2NKgE+onHf/VOOB90yln4Bmtmmc3VWNx2yMWWv+/7rCPwHuPodybf6eAe4HvqDqKbrdTW2O+UZgnjFmH4Axxt2PuzbHbIAgERGgCZUFv6x+YzqWMeZ7Ko+jOg6tYe5c8FsD+095nWH7Wl3buJO6Hs+dVI4O3FmNx2xbN/kq4LV6zOVMtfl77gSEiMhKEVkvIrfWWzrnqM0xzwa6ULlE6hZgqjGmon7iWcahNcze+fCtJFV87cxrTGvTxp3U+nhEZBiVBX+wUxM5X22O+UXgcWNMeeXgz+3V5ph9gH5Uri/RGPhRRH4yxiQ7O5yT1OaYRwEbgeFAB2CZiKw2xhxzcjYrObSGuXPBzwDanPI6irMXR69NG3dSq+MRkZ7AW8Blxpi8esrmLLU55njgE1uxbwGMEZEyY8xX9ZLQ8Wr7s51rjCkECkXke6AX4K4FvzbHfDsw3VSe3E4VkT1AZ+CX+oloCYfWMHc+pbMOiBWRGBHxA66ncgWuUy0AbrV90j2IygXUM+s7qAPVeMwi0haYB9zixqO9U9V4zMaYGGNMtDEmmsrFdu5z42IPtfvZng9cLCI+ttXkBgI76jmnI9XmmPdR+RsNIhIBxAG76zVl/XNoDXPbEb4xpkxEpgDfUvkJ//+MMdtE5F7b/teovGJjDJAKnKByhOC2annMfwWaA6/YRrxlxo1nGqzlMTcotTlmY8wOEVkCbAYqgLeMMVVe2ucOavn3/HfgXRHZQuWpjseNMW49bbKIfAwMBVqISAbwFOALzqlhOrWCUkp5CHc+paOUUqoOtOArpZSH0IKvlFIeQgu+Ukp5CC34SinlIbTgK6WUh9CCr5RSHuL/AWZsBblhSk1wAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "state_hat_split = np.hsplit(state_hat, 2)\n",
+    "\n",
+    "index = 100\n",
+    "\n",
+    "uhat = state_hat_split[0]\n",
+    "vhat = state_hat_split[1]\n",
+    "\n",
+    "err = np.linalg.norm(uhat[index]-u_true[index]) / np.linalg.norm(u_true[index])\n",
+    "print(err)\n",
+    "\n",
+    "plt.plot(x, state_hat_split[0][index], label='u_hat')\n",
+    "plt.plot(x, u_true[index],linestyle='dashed', label='u_true')\n",
+    "plt.legend()\n",
+    "plt.show()\n",
+    "\n",
+    "# print(u_true[500])\n",
+    "# print(t[500])\n",
+    "# print(c_squared)\n",
+    "\n",
+    "\n",
+    "\n",
+    "#print(uhat)"
+   ]
+  }
+ ],
+ "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.8.5"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "5f5b5ca3a926332263f32309627d0bc50222d2ae8afd567f76d7a8e89ef78637"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Archive/temp.py

+import numpy as np
+
+x = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
+
+y = np.hsplit(x, 2)
+print(y)
+
+
+# x = [np.array([0,1,2]), np.array([3,4,5])]
+# print(x)
+
+# y = np.array(x)
+# print(y)
+
+# x = np.array([0,1,2])
+# y = np.array([3,4,5])
+
+# z = np.concatenate((x,y), axis=0)
+# print(z)
+
+
+# a = np.array([[0,1], [4,5]])
+# b = np.array([[2,3], [6,7]])
+# c = np.array([[8,9], [12,13]])
+# d = np.array([[10,11], [14,15]])
+
+# x = np.bmat([[a, b], [c, d]])
+# print(x)
+
+# y = np.array([1,2,3,4])
+
+# print(x@y)
+
+
+# temp1 = np.concatenate((a,b), axis=1)
+# temp2 = np.concatenate((c,d), axis=1)
+# x = np.concatenate((temp1,temp2), axis=0)
+# print(x)
+# print(x@y)
+
+
+
+