Bai_12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tài liệu này mang giấy phép Creative Commons Attribution (CC BY).\n",
    "(c) Nguyễn Ngọc Sáng, Zhukovsky 12/2018.\n",
    "\n",
    "[@SangVn](https://github.com/SangVn) [@VnCFD](https://vncfdgroup.wordpress.com/)\n",
    "\n",
    "*Thực hành CFD với Python!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bài 12. Sơ đồ Godunov, bài toán Riemann\n",
    "\n",
    "### 1. Sơ đồ sai phân bảo toàn\n",
    "Xét phương trình Burgers không nhớt dạng không phân kì:\n",
    "$$\\frac {\\partial u}{\\partial t} + u\\frac{\\partial u}{\\partial x} = 0 \\qquad (1)$$\n",
    "Phương trình này có thể biểu diễn lại ở dạng phân kì:\n",
    "$$\\frac {\\partial u}{\\partial t} + \\frac{\\partial F}{\\partial x} = 0; \\ F=\\frac{u^2}{2} \\qquad (2)$$\n",
    "\n",
    "Sử dụng phương pháp sai phân hữu hạn với sơ đồ upwind ta có hai sơ đồ sai phân tương ứng `dạng không bảo toàn và bảo toàn`:\n",
    "$$\\frac{u_i^{n+1}-u_i^n}{\\Delta t} + u_i^n\\frac{u_i^n - u_{i-1}^n}{\\Delta x} = 0 \\qquad (3)\\\\\n",
    "\\frac{u_i^{n+1}-u_i^n}{\\Delta t} + \\frac{F_i^n - F_{i-1}^n}{\\Delta x} = 0 \\qquad (4)$$\n",
    "Xét thêm sơ đồ bậc hai MacCormack:\n",
    "$$u_i^* = u_i^n - \\frac{\\Delta t}{\\Delta x} [F_{i+1}^n - F_i^n]\\\\\n",
    "u_i^{n+1} = \\frac{1}{2}(u_i^n + u_i^*) - \\frac{\\Delta t}{2\\Delta x} [F_i^* - F_{i-1}^*] \\qquad (5)$$\n",
    "\n",
    "Tìm nghiệm bài toán `lan truyền sóng xung kích` đã được giới thiệu ở [bài 5, phần 1](https://nbviewer.jupyter.org/github/SangVn/CFD_Notebook_P1/blob/master/Bai_5.ipynb):\n",
    "$$u(x,0)=a,x≤x_0;u(x,0)=b,x>x_0;a>b$$\n",
    "Nghiệm chính xác:\n",
    "$$u(x,t)=a,x≤x_0+st;u(x,t)=b,x>x_0+st$$\n",
    "với s là vận tốc dịch chuyển gián đoạn và bằng (a+b)/2. Nhờ biết được nghiệm chính xác, ta có thể dễ dàng kiểm tra kết quả tính toán."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8VNXd+PHPmSWZmZA9YQlhCcoWyELYBQFLZRFkq0opLlRB0WprrTxgtdVHbR9r+zwqlIo8Kqj1UdD6U6wIFGQXyg7KKkuAsIZsZM9k5vz+mGRIQpIZINtMvu/XKzBz77n3nnsz+ebk3HO/R2mtEUII4V8MjV0BIYQQdU+CuxBC+CEJ7kII4YckuAshhB+S4C6EEH5IgrsQQvghCe5CCOGHJLgLIYQfkuAuhBB+yNRYB46KitIdO3ZsrMMLIYRP2rlz5yWtdbSnco0W3Dt27MiOHTsa6/BCCOGTlFInvSkn3TJCCOGHJLgLIYQfkuAuhBB+SIK7EEL4IQnuQgjhhyS4CyGEH5LgLoQQfkiCuxB+4Hz+eVakrmjsaogmpNEeYhJC1A2tNc9uepbt57czLHYYFpOlsaskmgBpuQvh4749+y3bzm9Do7lUeKmxqyOaCAnuQvgwp3by+q7XMSojgAR34SbBXQgf9vWJrzmUeYj74u8DJLiLKyS4C+Gj7A4783bPo2t4V+6Pvx+Q4C6ukOAuhI9aemQpZ/LO8GTvJ4mwRGBQBtIL0xu7WqKJkNEyQvigfHs+C/ctpF/rfgyKGYRSikhLpLTchZu03IXwQe/tf4/MokyeTHkSpRQAUdYoCe7CTYK7ED7mUuEl3tv/Hrd3uJ2E6AT38ihrFOkF0i0jXCS4C+FjFu5bSLGjmF/2+mWl5dG2aDIKMxqpVqKpkeAuhA85ffk0nxz5hEmdJ9ExtGOldZGWSDKKMnA4HY1TOdGkSHAXwofM2zMPs8HMo0mPXrUu2haNQzvIKs5qhJqJpkaCuxA+4mDGQb4+8TX3dr+XaFv0Veujra5l0jUjQIK7ED7j9V2vExoYys97/rza9VHWKAAZ6y4ACe5C+ITj2cf59uy3PNTzIYIDgqst4w7uMmJGIMFdCJ+w8cxGAEbHja6xTHlwzyiSbhkhwV0In7DxzEZuDruZ1kGtayxjMVkINgdLy10AEtyFaPIK7AXsurCLwW0HeywbZZOnVIWLBHchmrh/n/s3dqfdu+AuKQhEGQnuQjRxm85swmaykdIyxWPZKGuUjJYRgAR3IZo0rTWbzmyif5v+mI1mj+Wl5S7KSXAXogk7kXOCs/lnveqSAdeDTIWlheTb8+u5ZqKpk+AuRBNWPgTy1ra3elVexrqLchLchagjBfYC9lzcw+WSy3W2z01nNnFT6E20adHGq/LlwV26ZoTMxCTEddBacyr3FPvS97E3fS/70vdxJOsIDu2gb+u+vD3ibQzqxtpOBfYCdl7Yyc+6/czrbcrzy0hwFx6Du1LqXWAscFFr3bOa9VOB2WVv84BHtdZ767SWQjQBdqedFSdWsCJ1BfvS95FdnA1AkDmIhKgEHkp4CKd28vZ3b7Pk8BKmdJtyQ8fbdn6bawhkrHf97YA7oZgEd+FNy30x8Ffg/RrWnwCGaq2zlFKjgYVA/7qpnhCNr9hRzOc/fM6i/Ys4k3eGdsHtuK3dbSRFJ5EYnUin0E4YDUbA1aI/mHGQ13a+xuC2g2kX3O66j7vpzCasJqtXQyDLhQSEYDaYZTik8BzctdYblFIda1n/bYW3W4HYG69WzfZt/pLDqz6qz0MIL6my/3WFZbrsX33Vmpr2YcCIQqEwKAMGlOurrEvDqZ040Ti1EwcaJ04cZV+l2oFdl1KC63+7dlCi7di1AyeaFgYLwQYrIQYrIcpGiMFKkLJiUMpDrVxKdCkHSk6zrySVAorpZwgjJbAXHYzRqAP5KLZQoLawv+xMynf7cwIJyCnhg0P3Ma31bQSazZiMRlAKlBECgq58GQOvXMiKNORvX8nUoA4UrvgXhV7V2OX2H2wEpe8m5/xXlZYH9e+PKSrqGvYkfJnS2osfQFdw/2d13TJVyj0NdNNaT69h/cPAwwDt27fvffLkyWutL5++eD89/m/7NW8nRHMXNnkybf7zhcauhrhBSqmdWus+nsrV2Q1VpdRtwENAjR2EWuuFuLpt6NOnj+ffKtUY9MBvOdF7y3XVsTlT1TYPvaertsL1leUK5W59qrIjKWVw/V/LcTUah3a6WuXa4WqR61IcZe/RGqPBhFEZMSsjRuV6bVJGjMpIoCGAQGMgxlpuXDq0k2x7DpeKM7lUkklGSSbpJZlklWSTbc8h236ZUqqfli4pJJ5RLW8jLqj91TXXrkugy/5xvddorSmxOygsLmbx+Y84WXKKKbafEegIoqjEzoWcPC5lZNAqoJgBMSZ6RCoCS/OgKAeKL8OFA1BayBdRMbxj0bw94CVatk6q7VtzlT9s/QPnC84z70fz3MtOPzITR07ONe1H+LY6Ce5KqUTgbWC01rpe84226dCNNh261echRDOitSarOIuLBRe5kH+BCwUXyCnOYWi7oXQJ73JD+26ffxsTv5jImoiDvDPyHXdX0/bUTP629ihvHE6nxWUTUwe056HBcbQMtkBxHhz8ktV7/wdLSS7tVk2D9gMh6acQPwGsYR6Pa77QkcOnjhLYqZN7mTE0FGdhwQ2dj/AtNzzOXSnVHvgMuE9rfeTGqyREw1FKEWGJoFtEN4a2G8o9Xe9hRuKMGw7sAK2DWjOr7yx2XNjBx4c+di/v2zGCRT/vx1e/HMywrtH874bjDP7TWp77/DvOFRkp6DGeHQYHg7tPhuHPQ0EGfPkr+O+ucHiFx+NGWaPILMrE7rS7lxmsVnS+BPfmxJuhkB8Bw4AopVQa8DxgBtBaLwB+D0QCf1OuO0ql3vQHCdEcTLx5IqtOruL1Xa9za+ytlUbP9IgJ5a8/S+HEpXzeWn+MJdtPsyM1izk/cbqGQHYaBTEDYfCv4ewu+OIJV5DvsA0soTUeM8pWNmlHYYY7/7vBZqP0kgyPbE48tty11lO01m201matdazW+h2t9YKywI7WerrWOlxrnVz2JYFdiDJKKV4Y+AJGZeT3m3+PUzuvKhMXFcQrP0nk93f24ND5XL488g1Wk5XerXqX7wTa9obx8yDvAqx5sdZjVjdRtrJZcRZIy705kfQDQtSz1kGt+Y++/3FV90xVYxPaYDbC5rOb6d+6PwHGgMoF2vaGfg/D9nfgdM0jxqqbKNtgs+EsvJYBlcLXSXAXogFMuHkCg9oO4vVdr7Pm5Jpqy4QHBdC/i5N850VuiRlU/Y5+9BwEt3F1zzjs1RapNrhbbdJyb2YkuAvRAJRSvHjLi3QM6ciT657k2U3PkluSe1W5dm1PAWArreGREksI3PFnuLgftvy12iKRlkigcgoCabk3PxLchWggLW0t+fCOD3kk8RG+Ov4Vk5ZNYuu5rZXKXHLshZKWbDhY/dh7ALqPhW5jYd2fIPPEVavNRjPhgeFcKqgY3K1gt6NLSursfETTJsFdiAZkNpp5vNfjfDD6AyxGCzNWzeCVba9QWFpIYWkhuy7upKMthZX7z1NYUkuAH/0qGIzw1VNlT1JVVnWibIPNBiCt92ZEgrsQjSAhOoGldy5lavepfHjwQ+758h7e2/8eJc4SxncdTl5xKasPXqh5B6Ft4Ue/g2PfwPf/uGp1lKVycFdWK4D0uzcjEtyFaCRWk5U5/ebwvyP+l8LSQubvmY/VZOXepNtoE2rh891nat9BvxkQkwIr5kBBZqVV0bboq0bLgLTcmxMJ7kI0sgFtBvDZ+M+4p8s9/Lznz7GYAxmXFMP6I+lk5BXXvKHBCHe+4Qrsq5+vtCrSGsmlwkuUJwY0WMuCuzyl2mxIcBeiCQgJCOF3A3/Ho0mPAjChV1tKnZqvvjtX+4ZtEmHgY7DrfTh5Jft2tDUau9PunvLvSstdgntzIcFdiCaoe5sQurUO5v956poBGPYMhLaHf/7afXO1/CnV8omyDTbpc29uJLgL0URN6NWW3aeyOZmRX3vBgCAY9EtIPwTZrjkSIq1lY92LXDdVy1vuWvrcmw0J7kI0UeOSYlAKPt991nPh2LKUTmd3A9W03GW0TLMjwV2IJiomzMqAuEg+33MGjzOmtYwHY8CV4G6rnDxMlfe5F0jLvbmQ4C5EEzahVwwnLuWzN83DLEqmQGjVwx3cbSYbVpPVPRzSfUNVWu7NhgR3IZqwUT3bEGAyeB7zDhDTC87uBa1RShFljXIHdxUQAAaDjJZpRiS4C9GEhVrN/Lh7S77cexa74+pc8JXE9ILiHMg8DriyQ7q7ZZRyJQ+TlnuzIcFdiCZuQnJbMvJL2HTUw0xKMb1c/5d1zVRsuUPZVHsyWqbZkOAuRBM3rGtLwmxmz10z0d3AGFhpxMxVaX/lCdVmQ4K7EE1cgMnAHQltWLn/PHnFpTUXNJqhdQKc3QO4Wu65JbkUlRYBrhEzklum+ZDgLoQPmNirLUV2J2sPXay9YEwvOLcHnE73jEzlrXfpc29eJLgL4QMSY0NRCn64mFd7wZheUJIHGUevDu5Wq7TcmxEJ7kL4gECTkdYhFtIyPbS8K9xULX+QqXLL3UMqA+E3JLgL4SPahds4neUhuEd1AbMNzu6+aqJsg9WKlidUmw0J7kL4iNgIK6czPQRno8l1U/XcHsIDwzEow5WWe5D0uTcnEtyF8BHtwm1cyC2iuLSWuVWh7KbqXoxApCXSHdyV9Lk3KxLchfAR7SJsaA1nsjwE6JheYC+AS0eIskZV6nPXxcVoh4dfDsIvSHAXwke0j3Al/zrtTXAHd7/7lbS/Mo9qcyLBXQgf0S7ClZP9tKcRM5E3Q0AL94iZii13kHlUmwsJ7kL4iFbBFgKMBs8jZgxGaJMEZ3cTaYkksygTh9PhnmpPS2bIZkGCuxA+wmBQtA23kuZpxAy4umbOf0e0JQKHdpBVnFVhkmzplmkOJLgL4UNiw62eW+4AbZKhtIhoezHgepBJyVR7zYoEdyF8SLsIm+c+d3DfVI267MpFc6nwUoXZmKTl3hxIcBfCh7QLt5FVYCe3yF57wYhOEBhCVOYpwDVRtky117x4DO5KqXeVUheVUt/XsF4ppeYqpY4qpfYppVLqvppCCKg4YsZD69tggDZJRF04CEBGUUaFPncJ7s2BNy33xcCoWtaPBjqXfT0MvHnj1RJCVKddePlYd++6ZiwX9hNsbuFquUufe7PiMbhrrTcAmbUUGQ+8r122AmFKqTZ1VUEhxBXtyh9k8rbf3VFCVEAI6YVXumVkqr3moS763NsCpyu8TytbJoSoY+E2M0EBRtI8PaUKV26qagMZhRkoiwWQh5iai7oI7qqaZbragko9rJTaoZTakZ6eXl0RIUQtlFLej5gJ7wiWMKLsJaQXpqMMBplqrxmpi+CeBrSr8D4WOFtdQa31Qq11H611n+jo6Do4tBDNT6w3ed0BlIKYZKILcirPxiR97s1CXQT3ZcD9ZaNmBgA5WutzdbBfIUQ12kfYOJ1ZiNbV/oFcWUwvovIuUlhaSL493zUbk4yWaRZMngoopT4ChgFRSqk04HnADKC1XgAsB+4AjgIFwM/rq7JCCNdwyEK7g4z8EqJaBNZeOKYXUbtdY+LLR8xIy7158BjctdZTPKzXwC/qrEZCiFq5h0NmFngX3Mvyt18qvESUzSZT7TUT8oSqED6mnbd53QFC2xFtDgbKUxBIy725kOAuhI+JDfcyrzuAUkS3TABcE2XLaJnmQ4K7ED4mKNBEZFAAad6MmAFC2qRg1ppLeecxWGWS7OZCgrsQPii2bMSMN1TbFCIdDi5lHS0bLSMt9+ZAgrsQPqidt3ndAWJ6Ee5wkp17RkbLNCMS3IXwQe0ibJzNLsTh9GKse0gM4RjJLsrEYLOhCwvRTmf9V1I0KgnuQvigduE27A7N+ctFngsrRbg5iKzSQtc8qlqji7zYTvg0Ce5C+KAred2962IJDwgli1KUzKPabEhwF8IHlT/IdMrL4B5mjSJPgTa7nluUfnf/J8FdCB8UE2ZFKUjztuXeojUAhY4sQOZRbQ4kuAvhgwJMBmJCrd49pQqEh7QHIL80AwBnQX691U00DRLchfBRseFWr/vcw8I7AZBb4ppHQWZj8n8S3IXwUe0ivMzrDoRH3AxATqkrr7v0ufs/Ce5C+Kh24TYuXC6myO7wWDasrM89p9Q1HbKMlvF/EtyF8FHlwyHPZHsO1GGBYQBkOctuqMo8qn5PgrsQPsqd+teLfneTwUQIRjL0ZUBa7s2BBHchfJR70g5vR8wYrWQo15OpMlrG/0lwF8JHtQwOJMBk8Hqse1hAMJkmUAFmGS3TDEhwF8JHGQyK2DDvs0OGWyLINhgxWAJltEwzIMFdCB92LXndw4NakWk0oAKM8oRqMyDBXQgfdi153cOC27pa7mYlLfdmQIK7ED6sXYSN7AI7uUV2j2XDbdGUGBQYHTJaphmQ4C6ED3OPmPGia6Z8rLvDaJeWezMgwV0IH9a+fKy7F10zEZYIAEpVCc5CCe7+ToK7ED7sWibtCLO4Wu4lqhgtT6j6PQnuQviwUKuZ4EATaV48yBQeGA5AkdmJsyCvvqsmGpkEdyF8mFKK2AibVzMylbfcCwOU5JZpBiS4C+Hj2nmZ1z3YHIxJGSkIAGdRMVrrBqidaCwS3IXwce0ibKRlFXoM1kopwgLDyAtQ4HSiS0oaqIaiMUhwF8LHtQu3Umh3cCnPc7AOs4STazEDMmGHv5PgLoSPa3cNwyHDLeHkWE0AaAnufk2CuxA+7lryuocHhpMVaAQkp7u/k+AuhI+LDXeNdfdqOKQlnMwA12tnngyH9GcS3IXwcbYAE1EtArx7kCkwjCxTKQDOrHP1XTXRiLwK7kqpUUqpw0qpo0qpOdWsb6+UWquU2q2U2qeUuqPuqyqEqElsuM3rPveC8pZ7+ql6rpVoTB6Du1LKCMwHRgPxwBSlVHyVYs8BS7XWvYCfAn+r64oKIWrWzsu87uGB4RS7BsvgzDxTz7USjcmblns/4KjW+rjWugT4GBhfpYwGQspehwJn666KQghP2kdYOZtdSKnDWWu5MEsYReUt90zplvFn3gT3tsDpCu/TypZV9AJwr1IqDVgOPFHdjpRSDyuldiildqSnp19HdYUQ1WkfYaPUqTmXU1RrufDAcIrKWu46W34G/Zk3wV1Vs6zqo3BTgMVa61jgDuADpdRV+9ZaL9Ra99Fa94mOjr722gohqtUhMgiAkxm197uHW8KvtNxzMuq7WqIReRPc04B2Fd7HcnW3y0PAUgCt9RbAAkTVRQWFEJ51iHSNdU/NyK+1XFhgGA6jwmkAZ25OQ1RNNBJvgvt2oLNSKk4pFYDrhumyKmVOAcMBlFLdcQV3+ZtPiAbSKthCgMngMTukxWTBarJSGmDAmZ8LkjzMb3kM7lrrUuBxYCVwENeomP1KqReVUuPKiv0GmKGU2gt8BEzTknJOiAZjMCg6RNg46aHlDq5+d3uAAWeJA/IuNkDtRGMweVNIa70c143Sist+X+H1AWBQ3VZNCHEtOkTaPPa5g2vETEngeZylBsg+BcGtGqB2oqHJE6pC+In2EUGcyizwmPrXdVNVoR0Ksk82UO1EQ5PgLoSf6BBpo6DEQXpeca3lwgPDyTeD064g53StZYXvkuAuhJ8oHzHjqWsmLDCMApMDpza7umWEX5LgLoSfuJax7vkmBw6nBHd/JsFdCD/RNsyKQcEpDyNmwi2up1Qd5TdUhV+S4C6EnwgwGWgbbiXVU8u9LHmY065dwV1GLfslCe5C+JEOEUGc9PAgU1igK3mYtjuhtAjy5XlDfyTBXQg/0j7S5lW3TLFZoUqdaCeQLSNm/JEEdyH8SIcIG1kFdnIK7TWWqZQ8rFTGuvsrCe5C+JHyETOnaul3DwkIodjsSvbqCu5yU9UfSXAXwo+4x7pn1tw1YzKYUDbXpNpOY6gEdz8lwV0IP9I+wrsHmcxBLQBwWtpIcPdTEtyF8CNBgSaigwM9ZocMaBEKgA6MluDupyS4C+FnXKl/PeR1b+Ga8tgZEOXKLyNj3f2OBHch/Ez7SJvHSTuswREAOE3hYC+AAplyz99IcBfCz3SICOJcThFFdkeNZWzB4QA4DK6+dxkO6X8kuAvhZzpGuW6qnq6l9R4c2hKAYsoGvEu/u9+R4C6En/FmxExwiGv++vzSsgUS3P2OBHch/Iw79W8tLfew0JY4gaLCArCESgoCPyTBXQg/E24zE2wx1TocMswSTnEAFOXlQFh7abn7IQnuQvgZpZTHybLLc7rb83IhrIMEdz9kauwKVGS320lLS6OoqKixqyKaIYvFQmxsLGazubGrcsM6RARx4NzlGteHB4Zz1AymgjwI7QzH1rrGuivVgLUU9alJBfe0tDSCg4Pp2LEjSj5kogFprcnIyCAtLY24uLjGrs4Nax9pY+X+85Q6nJiMV/+BHmQOoihAYcvPd3XL2POhIBOCIhuhtqI+NKlumaKiIiIjIyWwiwanlCIyMtJv/mrsGGmj1Kk5l1P9+SilcASacBYWuYI7QFZqw1VQ1LsmFdwBCeyi0fjTZ699hOfJsp3WAFRhEbTs7lpw8UBDVE00kCYX3EXj+uMf/1jp/S233NKgx1+8eDGPP/54gx7TH5Wn/k2tLYGYJRBVXALhcWAOgvPfNVDtREOQ4O6ntNY4nc5r3q5qcP/222/rqkqiAbUOsRBgMtSeY8ZmxVRkB4MBWvWQ4O5nJLhXkJqaSvfu3ZkxYwY9evRgxIgRFBYWArBnzx4GDBhAYmIiEydOJCsrC4Bhw4Yxe/Zs+vXrR5cuXdi4cSMADoeDp59+moSEBBITE5k3bx4Aa9asoVevXiQkJPDggw9SXFwMQMeOHXn++edJSUkhISGBQ4cOAbB+/XqSk5NJTk6mV69e5ObmAvDnP/+Zvn37kpiYyPPPP1+p/o899hgpKSm89NJL/Md//If7/BYvXswTTzwBwIQJE+jduzc9evRg4cKFAMyZM4fCwkKSk5OZOnUqAC1auHKPTJ48meXLl7v3NW3aNP7xj3/gcDiYNWuWuy5vvfXWVdc1Pz+fMWPGkJSURM+ePVmyZAkA27dv55ZbbiEpKYl+/fq5z+3s2bOMGjWKzp07V6r/qlWrGDhwICkpKdx9993k5eW5r91vf/tbBg4cSJ8+fdi1axcjR47kpptuYsGCBe7tq7tm/spgULSPsNU61t1otWEqLss/0zoBLnwv2SH9ida6Ub569+6tqzpw4ID79QvLvtf3LPi2Tr9eWPb9Vces6MSJE9poNOrdu3drrbW+++679QcffKC11johIUGvW7dOa6317373O/2rX/1Ka6310KFD9VNPPaW11vqrr77Sw4cP11pr/be//U1PmjRJ2+12rbXWGRkZurCwUMfGxurDhw9rrbW+77779Guvvaa11rpDhw567ty5Wmut58+frx966CGttdZjx47VmzZt0lprnZubq+12u165cqWeMWOGdjqd2uFw6DFjxuj169frEydOaKWU3rJli9Za64sXL+qbbrrJfX6jRo3SGzdudNdHa60LCgp0jx499KVLl7TWWgcFBVW6JuXvP/vsM33//fdrrbUuLi7WsbGxuqCgQL/11lv6pZde0lprXVRUpHv37q2PHz9eaR+ffvqpnj59uvt9dna2Li4u1nFxcXrbtm1aa61zcnK03W7XixYt0nFxcTo7O1sXFhbq9u3b61OnTun09HR966236ry8PK211q+88or+z//8T/e1+9vf/qa11vrJJ5/UCQkJ+vLly/rixYs6Ojpaa61rvGZVVfwM+roHF23TI1+7+hzLff2bKXp3j2661FGq9fZ3tX4+ROvMEw1XQXFdgB3aixgrLfcq4uLiSE5OBqB3796kpqaSk5NDdnY2Q4cOBeCBBx5gw4YN7m0mTZpUqTzA6tWrmTlzJiaTa7RpREQEhw8fJi4uji5duni9n0GDBvHUU08xd+5csrOzMZlMrFq1ilWrVtGrVy9SUlI4dOgQP/zwAwAdOnRgwIABAERHR9OpUye2bt1KRkYGhw8fZtCgQQDMnTuXpKQkBgwYwOnTp93b12T06NF88803FBcX8/XXXzNkyBCsViurVq3i/fffJzk5mf79+5ORkXHVvhISEli9ejWzZ89m48aNhIaGcvjwYdq0aUPfvn0BCAkJcV+r4cOHExoaisViIT4+npMnT7J161YOHDjAoEGDSE5O5r333uPkySuZDMeNG+c+Vv/+/QkODiY6OhqLxUJ2dnat18xfdYgM4lRmAbqG1nhAUDCBpZBdmAmtE10LpWvGbzSpce4VPX9nj0Y5bmBgoPu10Wh0d8t4s43RaKS01JWJSWt91eiLmn7IatvPnDlzGDNmDMuXL2fAgAGsXr0arTXPPPMMjzzySKXtU1NTCQoKqrRs8uTJLF26lG7dujFx4kSUUqxbt47Vq1ezZcsWbDYbw4YN8zgE0GKxMGzYMFauXMmSJUuYMmWK+5zmzZvHyJEja9y2S5cu7Ny5k+XLl/PMM88wYsQIJkyYUOPolKrfg9LSUrTW3H777Xz00Ue1bmMwGCptbzAY3NtXd838WYdIGwUlDtLzimkZbLlqfWDZbEzZ2ReIbNkdlMEV3Lvf2dBVFfVAWu5eCA0NJTw83N2f/sEHH7hb8TUZMWIECxYscAfpzMxMunXrRmpqKkePHvV6P8eOHSMhIYHZs2fTp08fDh06xMiRI3n33Xfdfc5nzpzh4sWL1W4/adIkPv/8cz766CMmT54MQE5ODuHh4dhsNg4dOsTWrVvd5c1mM3a7vdp9/fSnP2XRokVs3LjRHcxHjhzJm2++6d7myJEj5OdX7uc9e/YsNpuNe++9l6effppdu3bRrVs3zp49y/bt2wHIzc11X6vqDBgwgM2bN7uvXUFBAUeOHKn12lV0LdfMX7QvGzFzqobhkJYWYQBk51yAABtEdpaWux9psi33pua9995j5syZFBQU0KlTJxYtWlRr+enTp3PkyBESExMxm83MmDGDxx9/nEWLFnH33XdTWlpK3759mTlzZq37ef3111mPnF04AAAc5ElEQVS7di1Go5H4+HhGjx5NYGAgBw8eZODAgYDrpuff//53jEbjVduHh4cTHx/PgQMH6NevHwCjRo1iwYIFJCYm0rVrV3c3DsDDDz9MYmIiKSkpfPjhh5X2NWLECO6//37GjRtHQECA+zxTU1NJSUlBa010dDSff/55pe2+++47Zs2ahcFgwGw28+abbxIQEMCSJUt44oknKCwsxGq1snr16hqvQ3R0NIsXL2bKlCnum9Avv/yyu4vLkxEjRlR7zVq2bOnV9r6oQ0T5cMgC+nSMuGq9LcS1LDc73bWgdQKc/neD1U/UM2865oFRwGHgKDCnhjL3AAeA/cD/edqnpxuqQjQGf/oMFtsdOm7OP/V/rzxU7fpTy5bqA1276WUrXDfy9cbXXDdV8zMasJbiWlFXN1SVUkZgPjAaiAemKKXiq5TpDDwDDNJa9wCerKtfPkKI6xNgMhATZq0xr3uL0LIJO3LK5k9t3dP1/4XvG6J6op550+feDziqtT6utS4BPgbGVykzA5ivtc4C0Fr7d2emED6iQ6SN1Br63MtvqBbmuZ7ZuDJiRoK7P/AmuLcFKk7Tkla2rKIuQBel1Gal1Fal1Ki6qqAQ4vp1iAziVA0PMhlsrj75otxs14IWLaFFK7mp6ie8uaFa3Xi1qmP6TEBnYBgQC2xUSvXUWmdX2pFSDwMPA7Rv3/6aKyuEuDYdImxkFdi5XGQnxFI5T73BagWgJL9C3vfWCRLc/YQ3Lfc0oF2F97HA2WrKfKG1tmutT+C6+dq56o601gu11n201n2io6Ovt85CCC91qGU4ZHnL3Z6Xe2Vh6wRIPwSlJQ1SP1F/vAnu24HOSqk4pVQA8FNgWZUynwO3ASilonB10xyvy4oKIa5deerf6rJDKqsruJcWVFjXOgGcdleAFz7NY3DXWpcCjwMrgYPAUq31fqXUi0qpcWXFVgIZSqkDwFpgltY6o74qLSpbvHgxZ89e+WNq+vTpHDjQcLm5U1NT6dmzZ4MdT3ivvOVeXV53g9X11KouqLBO0hD4Da8eYtJaLweWV1n2+wqvNfBU2Ze4AQ6Ho9qHkWqzePFievbsSUxMDABvv/12fVRN+KCgQBNRLQKr7ZZRRiOOABOGYjvFjmICjYEQ0QlMVgnufkDSD1Tx/vvvk5iYSFJSEvfddx8AJ0+eZPjw4SQmJjJ8+HBOnXLNFD9t2jR++ctfcsstt9CpUyc+/fRTAM6dO8eQIUNITk6mZ8+e7rQFtaWsffHFFxk8eDCvvvqq+0lScLWKExNdrakXX3yRvn370rNnTx5++GG01nz66afs2LGDqVOnkpycTGFhIcOGDWPHjh28+eabNab8/fvf/06/fv1ITk7mkUceweFwXHUt5syZQ3x8PImJiTz99NMAXLhwgYkTJ5KUlERSUpI737vD4ag2VfKxY8cYNWoUvXv35tZbb3WnMp42bRqPPvoot912G506dWL9+vU8+OCDdO/enWnTprnrUNM1E95zDYesfsSMtgZiKYGsorLhkAaj5Hb3F9486VQfXx6fUF0+W+t376jbr+Wza33y6/vvv9ddunTR6enpWusraXHHjh2rFy9erLXW+p133tHjx4/XWmv9wAMP6Lvuuks7HA69f/9+d3rdv/zlL/rll1/WWmtdWlqqL1++7DFl7Z/+9Cd3PZKSkvSxY8fc5cpT6pbXR2ut7733Xr1s2TKttSvt8Pbt293ryt/XlPL3wIEDeuzYsbqkpERrrfWjjz6q33vvvUrXIiMjQ3fp0kU7nU6ttdZZWVlaa63vueced5ri0tJSnZ2dXWuq5B/96Ef6yJEjWmutt27dqm+77Tb3tZs8ebJ2Op36888/18HBwXrfvn3a4XDolJQUvXv37lqvWX3xpydUy/16yW494I+rq123b8gg/eZPuuuDGQevLFz2K63/q53WZd970bTg5ROqklumgm+++Ya77rqLqCjXk3sREa7cG1u2bOGzzz4D4L777qvUGp4wYQIGg4H4+HguXLgAQN++fXnwwQex2+1MmDCB5ORk1q9f705ZC1BSUuLOcwK4k3oB3HPPPSxdupQ5c+awZMkS9+QWa9eu5dVXX6WgoIDMzEx69OjBnXfWnMGvYsrfzp07u1P+zp8/n507d7rT7RYWFl6VYyUkJASLxcL06dMZM2YMY8eOdV+j999/H3BlbAwNDSUrK6vaVMl5eXl8++233H333e79lueFAbjzzjtRSpGQkECrVq1ISEgAoEePHqSmppKWllbrNRPe6RARxGe7zlBkd2AxV+7yM9islVvu4LqpunMR5Jy+Mnm28DlNN7iPfqXBD6mrSdNbnYplKqaX1WUpfYcMGcKGDRv46quvuO+++5g1axbh4eG1pqytmKp38uTJ3H333UyaNAmlFJ07d6aoqIjHHnuMHTt20K5dO1544QWPaXrL91U15a/WmgceeID/+q//qnE7k8nEtm3bWLNmDR9//DF//etf+eabb2osX12qZKfTSVhYGHv27Kl1m5rS9BqNxlqvmfBO+U3V05kFdG4VXGmd0RZEYAFkF1d4JKXiTVUJ7j5L+twrGD58OEuXLiUjwzXQJzMzE3BNEv3xxx8D8OGHHzJ48OBa93Py5ElatmzJjBkzeOihh9i1a9c1pay96aabMBqNvPTSS+4WfXkgj4qKIi8vz92/DxAcHOyeoq6q6lL+Dh8+nE8//dSd8jYzM7PSxBcAeXl55OTkcMcdd/D666+7A/Tw4cN58803AVc/++XLl6lJSEgIcXFxfPLJJ4Drl9/evXtruXKV3WiaX+HSvpYRM+agFlhKNJlFmVcWtooHlPS7+zgJ7hX06NGDZ599lqFDh5KUlMRTT7kG/8ydO5dFixaRmJjIBx98wBtvvFHrftatW+ee8/Qf//gHv/rVryqlrE1MTGTAgAHum4vVmTx5Mn//+9+55557AAgLC2PGjBkkJCQwYcIEd5cKuG5Ozpw5031DtaLylL8nT55036iNj4/n5ZdfZsSIESQmJnL77bdz7ty5Stvl5uYyduxYEhMTGTp0KK+99hoAb7zxBmvXriUhIYHevXuzf//+Wq/Fhx9+yDvvvENSUhI9evTgiy++qLV8Rdd6zUT1Oka6/iqsLoGYOSgEi71Kyz0gCCJvluDu41R5V0JD69Onj96xY0elZQcPHqR79+6NUh8hwD8/g1prEl9YxbjkGP4wMaHSujO/eZqj337Nxtem8NyA566s+GQanNkJT0qAb2qUUju11n08lZOWuxB+TilF37gI1h66iNNZuTFnsNmwlKrKLXdw3VTNPgWFVZYLnyHBXYhmYGxiG87mFLH7dFal5QablcASXXm0DFy5qXqh9m430XRJcBeiGbg9vhUBJgNf7q18b0XZbASUOMmqeEMVXC13kH53HybBXYhmINhi5rau0Xz13TkcFbpmDFYbSkN+bpXulxatIChagrsPk+AuRDMxNjGG9Nxitp240kovT/tbmJdFpcEVSpXldt/X0NUUdUSCuxDNxPDuLbGajXy570oG0fIJO0zFpeTZq+TtkdzuPk2CexOzbt0696P+daEh0/HWdd1F3bIFmBjevSUrvj9PqcMJgCHI1XK32CG7qErXTKsEcJTAJXlwzBdJcPcz1WV3bMjtRdM2NjGGzPwSvj3megq7vOVuKYG0vLTKheWmqk+T4F5Bfn4+Y8aMISkpiZ49e7oTdu3cuZOhQ4fSu3dvRo4c6X6a8+jRo/z4xz8mKSmJlJQUjh07htaaWbNm0bNnTxISEtz7WLduHcOGDeOuu+6iW7duTJ061d3HuWLFCrp168bgwYPdCcqqs2bNGnr16kVCQgIPPvigOwlXxZTBn3zyCTt37iQpKYmBAwcyf/589/YOh4NZs2bRt29fEhMTeeutt9x1u+222/jZz37mTt5VcZtp06a5z6f8SdXqzh1caQuqO8earuGwYcP49a9/zZAhQ+jevTvbt29n0qRJdO7cmeeeu/JQjTcpioVnw7pG0yLQxJd7XV0z5X3uIc4AVqSuqFw48mYwWeDC9w1dTVEHmmzisD9t+xOHMuv2UfNuEd2Y3W92jetXrFhBTEwMX331FQA5OTnY7XaeeOIJvvjiC6Kjo1myZAnPPvss7777LlOnTmXOnDlMnDiRoqIinE4nn332GXv27GHv3r1cunSJvn37MmTIEAB2797N/v37iYmJYdCgQWzevJk+ffowY8YMvvnmG26++eZK2SErKioqYtq0aaxZs4YuXbpw//338+abb/Lkk08CYLFY2LRpEwCJiYnMmzePoUOHMmvWLPc+3nnnHUJDQ9m+fTvFxcUMGjSIESNGALBt2za+//574uLiKh13z549nDlzhu+/d/2AZ2e7/nSv7txPnz5d7Tn279+/xmsIEBAQwIYNG3jjjTcYP348O3fuJCIigptuuolf//rXXLx4kSVLlrB582bMZjOPPfYYH374Iffff/+1fQAEFrOREfGtWLn/PH+YmOCeaq9/aBLvpa7kmX7PYDG5ZmjCaIKW8XJT1UdJy72ChIQEVq9ezezZs9m4cSOhoaEcPnyY77//nttvv53k5GRefvll0tLSyM3N5cyZM0ycOBFwBVebzcamTZuYMmUKRqORVq1aMXToULZv3w5Av379iI2NxWAwkJycTGpqKocOHSIuLo7OnTujlOLee++ttm6HDx8mLi6OLl26APDAAw+wYcMG9/ryXwo5OTlkZ2czdOhQAPeEI+Ca+OL9998nOTmZ/v37k5GRwQ8//OCuW9XADtCpUyeOHz/OE088wYoVKwgJCanx3Gs6x5quYblx48a5r3+PHj1o06YNgYGBdOrUidOnT7NmzRp3iuLk5GTWrFnD8eMyRe/1GpvUhstFpWz8Id3dcu8flky+PZ9vTlXJ/Nk6wdUt00hpSsT1a7It99pa2PWlS5cu7Ny5k+XLl/PMM88wYsQIJk6cSI8ePdiyZUulsjVlQ6wtV0/VtLilpaUANaYZHjlyJBcuXKBPnz48/vjjtda9PGVwbWmLtdbMmzePkSNHVlq+bt26SimHKwoPD2fv3r2sXLmS+fPns3TpUl5//fUa61HdOWqtq72GVbepKfWvNymKhfcG3xxNqNXMl3vPMmRYGwA6BrSijbENy44t445Od1wp3DoBdr0Hl89AaGwj1VhcD2m5V3D27FlsNhv33nsvTz/9NLt27aJr166kp6e7A5Pdbmf//v2EhIQQGxvL559/DrgmoSgoKGDIkCEsWbIEh8NBeno6GzZsqDRtXlXdunXjxIkT7j7rirnLV65cyZ49e3j77bfp1q0bqamp7vS3H3zwgbt1XlFYWBihoaHuLpoPP/zQvW7kyJG8+eab2O12AI4cOUJ+fvXTr5W7dOkSTqeTn/zkJ7z00kvs2rWrxnOvSU3X0FvepCgW3gswGRjVozX/OnCBErPrl6kuLGJsp7FsObeFiwUXrxSWCbN9lgT3Cr777jv3Tbs//OEPPPfccwQEBPDpp58ye/ZskpKSSE5Ods8b+sEHHzB37lwSExO55ZZbOH/+PBMnTnTPwfqjH/2IV199ldatW9d4TIvFwsKFCxkzZgyDBw+mQ4cONZZbtGgRd999NwkJCRgMBmbOnFlt2UWLFvGLX/yCgQMHYi0bDQEwffp04uPjSUlJoWfPnjzyyCPuvx5qcubMGYYNG0ZycjLTpk1zt56rO/ea1HYNveFNimJxbcYmtSG/xMHG066x7c6CAsbdNA6ndvLV8a+uFGwV7/pfgrvPkZS/QlTQXD6DpQ4n/f+4hgGdInnilWlETHuAlr/5Dfcuv5d8ez6fjfvsSvfe3F6uSbMn/71xKy0ASfkrhKiFyWhgdEJr1hy6gLJacea7utXG3TSOo9lHOZh58Erhtn3g+HrIrfmvM9H0SHAXopkamxhDkd1JsdmCs2wGr5EdRxJgCGDZsWVXCg6dDaXFsHxWDXsSTZEEdyGaqb4dI2gVEshlTDjLboiHBoYyrN0wlh9fjt3huvFO1M0wbDYcXAYH/9mINRbXQoK7EM2U0aC4I6ENmU4jJXlXRk2Nv3k8WcVZbDqz6UrhW34JrXrC8qehKKcRaiuulQR3IZqxsYkxFBrMZF66ErAHxgwkwhJRuWvGaIY750LeBVj9QsNXVFwzCe5CNGMp7cPQVhu52VceyjMbzIzpNIZ1aesqZ4qM7Q39H4Ud78JJ74eyisYhwb0BrVu37prGd9fE2zS+f/zjHyu9v+WWW2742MK/KKWIjAqhNL+AtKwrD6KNv2k8pc5Svk79uvIGP3oWwtrDl78Ce1ED11ZcCwnuDaiugru3qgb3hjy28B3t20ZhdZQwYf5m/n3clQq4a0RXuoZ35ctjX1YuHBAEY19z5Xjf+N+NUFvhLQnuVVSXWvbkyZN07tzZ/Sj+rbfeyqpVqwCYMGECvXv3pkePHixcuNC9nxUrVpCSkkJSUhLDhw8nNTWVBQsW8Nprr5GcnMzGjRsrHXf9+vUkJyeTnJxMr169yM3NrTF9cEWLFy+ulHdm7NixrFu3jjlz5lBYWEhycjJTp04FoEWLFgDXlZZY+K/Q8BAijQ5CLGamvv1vFm8+gdaaO2+6k+8ufcfxnCpJ2m7+MSROhk3/AxcONE6lhUdNNnHY+T/+keKDdZvyN7B7N1r/9rc1rj948GCNqWVnz57NzJkz6d+/P/Hx8e5Uue+++y4REREUFhbSt29ffvKTn+B0OpkxYwYbNmwgLi6OzMxMIiIimDlzJi1atODpp5++6th/+ctfmD9/PoMGDSIvLw+LxVJr+mBPXnnlFf7617+yZ8+eq9Zda1riwYMHe3VM4ZsMNhsUFfH/fnELv1m6lxe+PMC+Mzn8ZtRoXtv5GsuOLuPJ3k9W3mjkf8HR1bDsCXhoFRiMjVN5USNpuVdQW2rZ6dOnk5uby4IFC/jLX/7i3mbu3LkkJSUxYMAATp8+zQ8//MDWrVsZMmSIO4VuRESEx2MPGjSIp556irlz55KdnY3JZKo1ffCNuNa0xMK/GWxWKC0l2AgL7+vDr3/chc92neGR9w7Tu+UA/nn8nzicVSZHCYp0BfgzO2D7241TcVErr1ruSqlRwBuAEXhba/1KDeXuAj4B+mqtd1RXxlu1tbDrS22pZQsKCtw5yPPy8ggODmbdunWsXr2aLVu2YLPZGDZsGEVFRbWm3a3JnDlzGDNmDMuXL2fAgAGsXr3aqy4Rk8mE0+l0vy8q8nyT63rSEgv/VZ7TXRcUYAwL41c/7kzPtiE8+fEeThd3whm9mW3ntzEwZmDlDRPvgX1LYPV/QnRX6DDINWRSNAkeW+5KKSMwHxgNxANTlFLx1ZQLBn4J/LuuK9lQakstO3v2bKZOncqLL77IjBkzANfEGOHh4dhsNg4dOsTWrVsBGDhwIOvXr+fEiRPu/QAEBweTm5tb7bGPHTtGQkICs2fPpk+fPhw6dMir9MEdO3Zkz5497pmQtm3b5l5nNpvd6X0ruta0xMK/qbLMoc4KaZuHd2/FF48PIsqQgnZY+O3a1/jzxs/YlXYKp7OscaCU6+aqKQDeHw+v3gRLH4A9/wd5F6s7lGhA3rTc+wFHtdbHAZRSHwPjgap3Ul4CXgWu7lD2ERVTyzqdTsxmM/Pnzyc1NZXt27ezefNmjEYj//jHP1i0aBE/+9nPWLBgAYmJiXTt2pUBAwYAEB0dzcKFC5k0aRJOp5OWLVvyr3/9izvvvJO77rqLL774gnnz5nHrrbe6j/3666+zdu1ajEYj8fHxjB49moCAALZs2UJSUhJKKXf64IpdJYMGDSIuLo6EhAR69uxJSkqKe93DDz9MYmIiKSkplfK6T5w4sdr9HjpUt/c4hG8ob7mX55cp1ym6BV/8YhhTPxnPMfsnvH/8ed4/DtoeTgs60aFFd3q1TGTIuFXEZe8hJO0brKlrMBxw5fknJgU6j4COg8AaDmYbBLRwjbgx28AgvcL1yWPK37KullFa6+ll7+8D+mutH69QphfwnNb6J0qpdcDTnrplJOWvaIqa42cw95u1pD32GObYWJQlsNoyTqcmr6SQgtIiShzFlOpiNOX98Ap02f9l/yo0BjTKtYLqOil1tUubh/O92jLlr6uva1tvU/5603Kv/vty5UAG4DVgmheVehh4GKB9+/ZeHFoIUd9sKb0InTDhqpZ7VdYq74tLi0gvyCazKBu7wzUdohONU7u+tHaCdmDQdtDlgV6jqryuEE4qUYCuYZ03mvKvDlNYZP0fw4syaUC7Cu9jgbMV3gcDPYF1ZTcRWwPLlFLjqrbetdYLgYXgarnfQL2FEHXEGBZGzCvXNz/tTXVcF1F3vOn02g50VkrFKaUCgJ8C7oxCWuscrXWU1rqj1rojsBW4KrALIYRoOB6Du9a6FHgcWAkcBJZqrfcrpV5USo2r6wrJE5GischnT/gTr8a5a62XA8urLPt9DWWHXW9lLBYLGRkZREZGXvM4cSFuhNaajIwMLBZLY1dFiDrRpNIPxMbGkpaWRnp6emNXRTRDFouF2NjYxq6GEHWiSQV3s9nsfmRfCCHE9ZOnCIQQwg9JcBdCCD8kwV0IIfyQx/QD9XZgpdKBk9e5eRRwqQ6r4wvknJsHOefm4UbOuYPWOtpToUYL7jdCKbXDm9wK/kTOuXmQc24eGuKcpVtGCCH8kAR3IYTwQ74a3Bd6LuJ35JybBznn5qHez9kn+9yFEELUzldb7kIIIWrRpIO7UmqUUuqwUuqoUmpONesDlVJLytb/WynVseFrWbe8OOenlFIHlFL7lFJrlFIdGqOedcnTOVcod5dSSiulfH5khTfnrJS6p+x7vV8p9X8NXce65sVnu71Saq1SanfZ5/uOxqhnXVFKvauUuqiU+r6G9UopNbfseuxTSqVUV+66aa2b5BdgBI4BnYAAYC8QX6XMY8CCstc/BZY0dr0b4JxvA2xlrx9tDudcVi4Y2IBrvoA+jV3vBvg+dwZ2A+Fl71s2dr0b4JwXAo+WvY4HUhu73jd4zkOAFOD7GtbfAXyNa9KoAcC/6/L4Tbnl7p6YW2tdApRPzF3ReOC9stefAsOVb+cK9njOWuu1Wuvyaeq34poZy5d5832GKxOwFzVk5eqJN+c8A5ivtc4C0FpfbOA61jVvzlkDIWWvQ6k845vP0VpvADJrKTIeeF+7bAXClFJt6ur4TTm4twVOV3ifVras2jLaNalIDlD/kxPWH2/OuaKHcP3m92Uez7lsAvZ2Wut/NmTF6pE33+cuQBel1Gal1Fal1KgGq1398OacXwDuVUql4Zo/4omGqVqjudaf92vSpFL+VlHrxNzXUMaXeH0+Sql7gT7A0HqtUf2rswnYfYg332cTrq6ZYbj+OtuolOqptc6u57rVF2/OeQqwWGv930qpgcAHZefsrP/qNYp6jV9NueXuaWLuSmWUUiZcf8rV9mdQU+fNOaOU+jHwLK65aosbqG715VomYE/F1Te5zMdvqnr72f5Ca23XWp8ADuMK9r7Km3N+CFgKoLXeAlhw5WDxV179vF+vphzca52Yu8wy4IGy13cB3+iyOxU+yuM5l3VRvIUrsPt6Pyw0zwnYvflsf47r5jlKqShc3TTHG7SWdcubcz4FDAdQSnXHFdz9eVq2ZcD9ZaNmBgA5Wutzdbb3xr6j7OFu8x3AEVx32Z8tW/Yirh9ucH3zPwGOAtuATo1d5wY459XABWBP2deyxq5zfZ9zlbLr8PHRMl5+nxXwP8AB4Dvgp41d5wY453hgM66RNHuAEY1d5xs834+Ac4AdVyv9IWAmMLPC93h+2fX4rq4/1/KEqhBC+KGm3C0jhBDiOklwF0IIPyTBXQgh/JAEdyGE8EMS3IUQwg9JcBdCCD8kwV0IIfyQBHchhPBD/x99lWciNnSjKwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f353a48e6d0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#chia lưới trên đoạn [0, 1]\n",
    "nx = 41\n",
    "x, dx = np.linspace(0, 1, nx, retstep=True)\n",
    "\n",
    "#bước thời gian\n",
    "r = 0.4 #điều kiện ổn định: u_max dt/dx <= 1\n",
    "dt = r*dx\n",
    "\n",
    "#điều kiện ban đầu\n",
    "uL = 1.2\n",
    "uR = 0.4\n",
    "\n",
    "u0 = np.zeros(nx)\n",
    "i = int(0.5/dx)\n",
    "u0[:i] = uL\n",
    "u0[i:] = uR\n",
    "\n",
    "#sơ đồ sai phân 3\n",
    "def nonconservative_scheme(nt):\n",
    "    u = u0.copy()\n",
    "    un = np.zeros(nx)\n",
    "    for n in range(nt):\n",
    "        un = u.copy()\n",
    "        u[1:] = un[1:] - dt/dx*un[1:]*(un[1:] - un[:-1])\n",
    "    return u\n",
    "\n",
    "#sơ đồ sai phân 4 (F=u^2/2)\n",
    "def conservative_scheme(nt):\n",
    "    u = u0.copy()\n",
    "    un = np.zeros(nx)\n",
    "    for n in range(nt):\n",
    "        un = u.copy()\n",
    "        u[1:] = un[1:] - dt/(2*dx)*(un[1:]**2 - un[:-1]**2)\n",
    "    return u\n",
    "\n",
    "#sơ đồ bậc hai MacCormack\n",
    "def second_order_scheme(nt):\n",
    "    u = u0.copy()\n",
    "    u1 = u0.copy()\n",
    "    un = np.zeros(nx)\n",
    "    for n in range(nt):\n",
    "        un = u.copy()\n",
    "        u1[1:-1] = un[1:-1] - dt/(2*dx)*(un[2:]**2 - un[1:-1]**2)\n",
    "        u[1:-1] = (un[1:-1] + u1[1:-1])/2 - dt/(4*dx)*(u1[1:-1]**2 - u1[:-2]**2)\n",
    "    return u    \n",
    "\n",
    "#so sánh nghiệm tìm được\n",
    "nt = 25\n",
    "u_1 = nonconservative_scheme(nt)\n",
    "u_2 = conservative_scheme(nt)\n",
    "u_3 = second_order_scheme(nt)\n",
    "#nghiệm chính xác: vận tốc lan truyền sóng xung kích s=(uL+uR)/2. = 0.8\n",
    "u_exact = np.zeros(nx)\n",
    "i = int((0.5 + 0.8*nt*dt)/dx)\n",
    "u_exact[:i] = uL\n",
    "u_exact[i:] = uR\n",
    "\n",
    "plt.plot(x, u_1, x, u_2, x, u_3, x, u_exact)\n",
    "plt.legend(['nonconservative scheme', 'conservative scheme', 'second-order scheme', 'exact solution'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta thấy, với sơ đồ không bảo toàn, `sóng xung kích chuyển động sai vận tốc`, khác với trường hợp sơ đồ bảo toàn. Sử dụng sơ đồ bảo toàn là cần thiết. Ta cũng thấy rằng sơ đồ bậc hai cho gián đoạn có độ dốc lớn nhưng lại làm xuất hiện dao động.\n",
    "\n",
    "Như đã nói ở bài trước, ngoài sử dụng phương pháp sai phân hữu hạn biến đổi phương trình vi phân ở dạng phân kì để thu được sơ đồ sai phân bảo toàn, ta có thể sử dụng phương pháp `thể tích hữu hạn`.\n",
    "\n",
    "### 2. Phương pháp Godunov\n",
    "\n",
    "Sử dụng phương pháp `thể tích hữu hạn` ta thu được sơ đồ sai phân:\n",
    "$$\\frac{du}{dt} + \\frac{F_{i+1/2} - F_{i-1/2}}{\\Delta x} = 0$$\n",
    "với xấp xỉ thời gian upwind bậc 1 ta có:\n",
    "$$\\frac{u_i^{n+1} - u_i^n}{\\Delta t} + \\frac{F_{i+1/2} - F_{i-1/2}}{\\Delta x} = 0$$\n",
    "ở đây $u_i$ là giá trị trung bình trên đoạn [i-1/2, i+1/2]: \n",
    "$$u_i^n = \\frac{1}{\\Delta x}\\int_{i-1/2}^{i+1/2} u(x, t_n) dx$$\n",
    "\n",
    "Trong phương pháp Godunov hàm dòng F xác định như sau:\n",
    "$$F_{i-1/2}^* = F(u_{i-1/2}^*(u_{i-1/2}^L, u_{i-1/2}^R)); \\quad F_{i+1/2}^* = F(u_{i+1/2}^*(u_{i+1/2}^L, u_{i+1/2}^R))$$\n",
    "\n",
    "với $u^*$ là nghiệm của `bài toán Riemann - phân rã gián đoạn` trên bề mặt thể tích hữu hạn: cho giá trị u tại hai bên trái phải của bề mặt là uL, uR, tìm nghiệm $u^*(u_L, u_R)$ ở trên bề mặt. Godunov lấy uL, uR chính là giá trị trung bình u ở hai thể tích hữu hạn hai bên bề mặt đang xét $u_{i+1/2}^*(u_i, u_{i+1})$.\n",
    "\n",
    "\n",
    "**Bài toán Riemann**\n",
    "\n",
    "Chúng ta đã xem xét qua bài toán Riemann ở phần một, bây giờ ta sẽ nghiên cứu cụ thể hơn. Giả thiết bài toán:\n",
    "$$u_t + F(u)_x = 0 \\qquad (6)\\\\\n",
    "u(t=0)=u_0 = \n",
    "\\begin{cases}\n",
    "u_L; \\quad x<0 \\\\ \n",
    "u_R; \\quad x>0\n",
    "\\end{cases}\n",
    "$$\n",
    "<img src='img/Bai_12.png' width=400>\n",
    "\n",
    "Để giải bài toán này ta cần biết về `đường đặc trưng và sóng`.\n",
    "\n",
    "**Đường đặc trưng**\n",
    "\n",
    "Xét phương trình $u_t + au_x = 0 \\qquad (7)$. \n",
    "\n",
    "*Đường đặc trưng được định nghĩa là đường trong không gian mà dọc theo nó phương trình đạo hàm riêng (PDE) trở thành phương trình vi phân thường (ODE).*\n",
    "\n",
    "Xem $x=x(t)$, ta có:\n",
    "$$\\frac{du}{dt} = \\frac{\\partial u}{\\partial t} + \\frac{\\partial u}{\\partial x} \\frac{dx}{dt}$$\n",
    "Nếu $x'=a$ ta có:\n",
    "$$\\frac{du}{dt} = \\frac{\\partial u}{\\partial t} + a\\frac{\\partial u}{\\partial x} = 0$$\n",
    "Nghĩa là $u$ là hằng số dọc theo đường $\\frac{dx}{dt} = a$. Ở đây $a$ được gọi là vận tốc đặc trưng và bằng độ dốc của đường $x(t)$. Trong trường hợp tổng quát $a=a(x,t)$.\n",
    "\n",
    "Giả thiết tại thời điểm ban đầu $x(t=0)=x_0$, ta có nghiệm phương trình $x'=a$: $x=x_0+at$. Dọc theo đường đặc trưng u là hằng số nên ta có:\n",
    "$$u(x,t) = u(x(0),0)=u_0(x_0)=u_0(x-at)$$\n",
    "Nghiệm u chính là $u_0$ ban đầu được dịch chuyển với vận tốc $a$ sang trái hoặc sang phải tùy theo dấu của $a$. Quá trình này ta gọi là `dịch chuyển sóng`.\n",
    "\n",
    "<img src='img/Bai_12_2.png' width=600>\n",
    "\n",
    "Quay lại với bài toán Riemann, phương trình (6) có thể viết lại ở dạng: \n",
    "$$u_t + Au_x = 0$$ \n",
    "với $A=\\frac{\\partial F}{\\partial u}$ gọi là ma trận Jacobian. `Giá trị riêng` của ma trận A chính là vận tốc đường đặc trưng hay vận tốc dịch chuyển sóng. Đối với phương trình trình (7): $A=|a|, \\lambda = a$. Với a là hằng số, sóng dịch chuyển với vận tốc không đổi nên hình dạng không đổi. Đối với phương trình (2): $A=|u|, \\lambda=u$, sóng dịch chuyển với vận tốc phụ thuộc thời gian và không gian điều này có thể dẫn đến sự giao cắt, phân tách đường đặc trưng làm xuất hiện sóng va đập và sóng giãn như đã giới thiệu ở phần 1.\n",
    "\n",
    "Trong trường hợp bài toán Riemann cho phương trình (2), ta thấy sóng ở hai bên trái phải chuyển động với hai vận tốc khác nhau uL và uR.\n",
    "\n",
    "Tại $x_{i+i/2}$ 'cấu hình' đường đặc trưng có nhiều dạng phụ thuộc mối liên hệ uL, uR, ví dụ:\n",
    "\n",
    "<img src='img/Bai_12_3.png'>\n",
    "\n",
    "Hình 1 ngoài cùng bên trái: $u_L<0<u_R$, hình 2:$u_R>u_L>0$, hình 3: $u_L<u_R<0$, hình 4 ngoài cùng bên phải: $u_L>0>u_R, (u_L+u_R)/2<0$. **Chú ý:$u=dx/dt$.**  \n",
    "\n",
    "Căn cứ vào vị trí $x_{i+1/2}$ thuộc vùng 1, 2 hay 3 và mối liên hệ uL, uR - sóng nào dịch chuyển nhanh hơn, về phía nào để ta xác định $u^*$. Ví dụ như trên hình, hình 1 - $x_{i+1/2}$ nằm giữa (1) và (2) - lấy bằng 0, hình 2 - nằm trong (1) - lấy uL, hình 3 - nằm trong (2) - lấy uR, hình 4 - nằm trong vùng giao thoa, ở đây ta lấy sóng nào chuyển động nhanh hơn - ảnh hưởng mạnh hơn, chạy tới ${i+1/2}$ nhanh hơn. Tổng hợp tất cả các trường hợp ta có:\n",
    "\n",
    "Nếu $u_L > u_R$ - sóng xung kích:\n",
    "$$\n",
    "u* =\n",
    "\\begin{cases}\n",
    "u_L; \\ if\\quad \\frac{u_L+u_R}{2}>0 \\\\ \n",
    "u_R; \\ in\\ other\\ case\n",
    "\\end{cases}\n",
    "$$\n",
    "Nếu $u_L \\leq u_R$ - sóng giãn:\n",
    "\n",
    "$$\n",
    "u* =\n",
    "\\begin{cases}\n",
    "u_L; \\ if\\quad u_L>0 \\\\ \n",
    "u_R; \\ if\\quad u_R<0 \\\\\n",
    "0; \\quad if\\quad u_L \\leq 0 \\leq u_R\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.86 ms, sys: 0 ns, total: 3.86 ms\n",
      "Wall time: 2.73 ms\n",
      "CPU times: user 1.42 ms, sys: 186 µs, total: 1.61 ms\n",
      "Wall time: 1.15 ms\n"
     ]
    }
   ],
   "source": [
    "#tìm nghiệm bài toán Riemann\n",
    "def decay_godunov(uL, uR):\n",
    "    ustar = 0\n",
    "    if uL > uR:\n",
    "        if (uL + uR)/2. > 0: ustar = uL\n",
    "        else: ustar = uR\n",
    "    else:\n",
    "        if uL > 0: ustar = uL\n",
    "        elif uR < 0: ustar = uR\n",
    "        else: ustar = 0\n",
    "    return ustar\n",
    "\n",
    "#tính hàm dòng F(u*)\n",
    "def flux(u):\n",
    "    return u**2/2.\n",
    "\n",
    "#sơ đồ Godunov\n",
    "def godunov_method(nt):\n",
    "    u = u0.copy()\n",
    "    un = np.zeros(nx)\n",
    "    for n in range(nt):\n",
    "        un = u.copy()\n",
    "        for i in range(1, nx-1):\n",
    "            u_im05 = decay_godunov(un[i-1], un[i]) #u_{i-1/2}\n",
    "            u_ip05 = decay_godunov(un[i], un[i+1]) #u_{i+1/2}\n",
    "            u[i] = un[i] - dt/dx*(flux(u_ip05) - flux(u_im05))\n",
    "    return u\n",
    "\n",
    "#có thể viết lại code như sau\n",
    "def godunov_scheme(nt):\n",
    "    u = u0.copy()\n",
    "    un = np.zeros(nx)\n",
    "    for n in range(nt):\n",
    "        un = u.copy()\n",
    "        ustar = np.array([decay_godunov(ul, ur) for ul, ur in zip(un[:-1], un[1:])])\n",
    "        u[1:-1] = un[1:-1] - dt/dx*(flux(ustar[1:]) - flux(ustar[:-1]))\n",
    "            \n",
    "    return u\n",
    "\n",
    "%time u_g = godunov_method(nt)\n",
    "%time u_g_2 = godunov_scheme(nt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So sánh thời gian chạy hai đoạn code khác nhau, rõ ràng sử dụng tính toán mảng numpy nhanh hơn nhiều."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8VNXd+PHPN8mEJJCEJUBAkEXZwiJLgkprJGjZtAjqE8DaihuPKGKtivq0QkCtwmOV8pNHRKqobSERq3VBbZVQQMEmFEQWWQSKKQgBZAmQ/fz+mEnIMslMYCYzc+/3/XrNi7n3nnvmexn45uTcc88RYwxKKaWsJSzQASillPI9Te5KKWVBmtyVUsqCNLkrpZQFaXJXSikL0uSulFIWpMldKaUsSJO7UkpZkCZ3pZSyoIhAfXBCQoLp3LlzoD5eKaVC0oYNG44YY1p7Khew5N65c2dyc3MD9fFKKRWSROTf3pTTbhmllLIgTe5KKWVBmtyVUsqCNLkrpZQFaXJXSikL8pjcReRVETksIlvqOP4zEdnsen0hIpf5PkyYOxeys10bGRmAc3vu3As/HqhzQ7Vujcu35yrlF8aYel9AKjAQ2FLH8SFAC9f7UcCXnuo0xjBo0CDTECtXGpOQ4PzTQPXtCzweqHNDtW6Ny7fnzplz7r2ZObOyvjlzvPqvoWwGyDVe5FgxXiyzJyKdgQ+MMX08lGvh+iFwkac6k5OTTUPHuWdnw803lXBDxyW8s2ciT/33Gwzqubvy+IZvLuU3L/+CcV2XNvh4oM4N1bo1rrrPvbHLUv6ydyLPTHmTlKRvEdexR//vVtI7vkHWd79g4cOZXDtwP9GOcNZt6cL4J28k64l3SPvgLrKvX0z6k+Oc2wO8GtLsnT43QZuevqtPBYSIbDDGJHsq5+uHmO4EPqrroIhMBiYDXHzxxQ2rOSODtFmzuLHLSyzefA+/uWouU5r+Fr47V2RwU/h+QAFPrZne4OOBOjdU69a4PJ87Ofq3sNd5rH8T2N//ZOWxm8t+CznOY0OBZTdkkT7jNe4ZNIuFM0aSefME0k6tgdX4iIGTB2DsAl9VqIKdN817oDN1dMtUKZMGbAdaeVNnQ7tljDn36+0TzKr2a60vjgfq3FCtW+Oq69xy82tmm4SEcvPRJ6XmVGGJOX6m2Lz7YbFp1arcPBD5nGnRssy88PoP5t2NeeYPa/aY//2fReaxEfeZAT1WGjAmfshOM+bnz5u//Ga+KSwp9U23zfxBxmTd1oATVLDCy24ZnyR3oB/wLdDdm/qM9rmHdN0al2/PrXr+o2FPm9jmpab/PRtNp0c/MIOe/LuZMne/adWqvN7zPXrpx8b8Kb0BJ6hg1WjJHbgY2A0M8aauildDk7un1suFHA/UuaFat8bl23PdJ/9y88LrP5jbX/un6fToB6b9LetMVGyxeShiTsMTuzHGLB5uzJLrG3iSCkbeJnePN1RFZCnObsEE4BAwE3C4unQWishi4Cag4s5PqfGis/98bqgqZUVz50JKCqSl4RwqmZFBdjbk5MD06bB3xjO88cVe5p/+FcfW9+QJZjObmTBzZuXQSo/eGAvFBXDXp368EtUYvL2h6tVoGX/Q5K6U97Kz4adjS4novYewHZfy9vIw5w8Dby29BY7/G6Z87rcYVePwNrnrE6pKBbnsbEhPh6WZ0CnlX6TcuYP09CoPRnnDEQ0lZ/wWowo+mtyVCnI5OZCVBT8dGcF9scfZEbaHGS+cICenAZU4oqGk0G8xquCjyV2pIDd9OpVdML+YeTft46P49PjXPPJIA7pUteVuO5rclQohUY5wHvxJd77KO8FHW773/kRHNJSc9V9gKuhoclcqxNw4sAPd2zbjfz/ZQUlZuXcnOWKgrAjKy/wbnAoamtyVCjHhYcL0ET3Ze+Q0mTnfeT4BnC130Na7jWhyVyoEXdOrDSmdWzDv012cLir1fIIjxvmnJnfb0OSuVAgSER4b1ZMjBUW8unav5xMqWu6lmtztQpO7UiFqUKeWDE9qy8ur93DsdHH9hbVbxnY0uSsVwprv7svRXfG8uHJ3/as8RVQkdx0OaRea3JUKYaPSmnDiw2ReXnaS7+YtrHyaNSWlRkFtuduOJnelQlhaGrz5x3IOvTuQ/2r9BunpzqdZa807U3lDVVvudqHJXalQlpHBTWOa8OM2n/Dl7uHceey3pA2T2rNFasvddjS5KxXKMjLIXmnYdOI64ofs4uXmj5G90rhJ7joU0m40uSsVwir62F96tYjmV+3k/qeOuJ8xUlvutuMxuYvIqyJyWES21HG8p4isE5EiEXnY9yEqpepSMWPk+DFRRFFGRIcjZGVRe8ZITe62E+FFmSXAi8AbdRw/BkwDxvooJqWUl6ZPr3gndO/Qkp2HTjHjp+5uqOpQSLvx2HI3xqzGmcDrOn7YGJMDlPgyMKVUw3RvG8s3359yfzAiyvmnttxtQ/vclbKInomxHCko4mhBUe2DIs6bqtpyt41GTe4iMllEckUkNz8/vzE/WinL6942FoAdh+poveuc7rbSqMndGLPIGJNsjElu3bp1Y360UpbXM9GZ3HfW1TXjiNHkbiPaLaOURbSObULzGEf9LXedFdI2PI6WEZGlwFAgQUTygJmAA8AYs1BEEoFcIA4oF5FfAknGmJN+i1opVYuI0KNtLDvqbLlrt4ydeEzuxpiJHo5/D3TwWURKqfPWIzGWv/zrPxhjEJHqByN0kWw70W4ZpSykR2IsBUWl/Oe4mxa6ttxtRZO7UhbSo2LEjLuuGb2haiua3JWykO6J9QyHdGi3jJ1oclfKQuKiHLSPj3I/HFK7ZWxFk7tSFtMjsY5pCLRbxlY0uStlMd0TY9mTf5qSsvLqB7Tlbiua3JWymJ6JsRSXlbPvyOnqBxzRUFYE5WWBCUw1Kk3uSllMnXPM6JzutqLJXSmLuaR1M8LDpPZwSF1qz1Y0uStlMVGOcDq3inGT3HXBDjvR5K6UBfVMjNNuGZvT5K6UBXVvG8v+Y2c4U1x6bmdFt4zODGkLmtyVsqAeibEYA7sOFZzbqUvt2Yomd6UsqIe7aQgqb6hqn7sdaHJXyoIubhlDlCOs+k1V7XO3FU3uSllQeJjQrU0sO9223DW524HH5C4ir4rIYRHZUsdxEZH5IrJbRDaLyEDfh6mUaqhac8zoUEhb8ablvgQYWc/xUUA312sy8NKFh6WUulA92saSf6qIY6eLnTu0W8ZWPCZ3Y8xq4Fg9RW4A3jBO64HmItLOVwEqpc5P5U3Vita7dsvYii/63C8Cvquynefap5QKoIrkXtnvHtEEEE3uNuGL5C5u9hm3BUUmi0iuiOTm5+f74KOVUnVpE9uE5jGOc/3uIroak434IrnnAR2rbHcADrgraIxZZIxJNsYkt27d2gcfrZSqi4jQvW3NETM6p7td+CK5vwf8wjVq5grghDHmoA/qVUpdoJ6Jsez8/hTGuH6Z1tWYbCPCUwERWQoMBRJEJA+YCTgAjDELgRXAaGA3cAa43V/BKqUapnvbWE4VlXLgRCEXNY/Wbhkb8ZjcjTETPRw3wH0+i0gp5TM9K26qfn+qSnLXlrsd6BOqSllYN9eqTN9UHQ6ps0LagiZ3pSwsPtpB+/ioKsMho7TlbhOa3JWyuO5VpyHQG6q2ocldKYvrkRjLt4cLKC0r1xuqNqLJXSkLmzsXivcnUFxWzr5Zc8ERTfb2fsydG+jIlL95HC2jlApdKSlw082tiBzeih0bP+C76J+T/ubNZL0f6MiUv2nLXSkLS0uDpcsM+X8dwO/PPET6nFvJ+q87SUsLdGTK3zS5K2VlGRmMGB5Bh15b+ft345gSP4+0TtmQMSPQkSk/0+SulJVlZJC90nBoZ1969fmcl04+SPbeq+Dx6YGOTPmZJnelLCw7G9LTYfQv95EwZDtZT68kffkSsj8rDXRoys80uStlYTk5kJUFKVeW8n2Ltgy94geybp5ETo7bWbmVhehoGaUsbLqr92Xf51EUE8ZpE0lalzWkTTkEtAhobMq/tOWulA0kxkUBcKzY1Z7TB5ksT5O7UjbQNt6Z3I8Wuf7L6xQElqfJXSkbaOdK7vmF4c4dOjOk5WlyV8oGWjdrQpjAobOuJY+15W55XiV3ERkpIjtEZLeIPObmeCcR+UxENovIKhHp4PtQlVLnKyI8jIRmTfj+jCZ3u/CY3EUkHFgAjAKSgIkiklSj2HPAG8aYfsBs4BlfB6qUujCJ8VHkna5I7npD1eq8abkPBnYbY/YYY4qBZcANNcokAZ+53me7Oa6UCrDEuCj+45rWXVvu1udNcr8I+K7Kdp5rX1VfATe53o8DYkWkVc2KRGSyiOSKSG5+fv75xKuUOk+J8VH8+5Tr4SVtuVueN8ld3Oyr+Xjbw8DVIrIRuBr4D1Dr+WZjzCJjTLIxJrl169YNDlYpdf7axkWRXwgG0Za7DXjzhGoe0LHKdgfgQNUCxpgDwI0AItIMuMkYc8JXQSqlLpxzOKRgHNGIJnfL86blngN0E5EuIhIJTADeq1pARBJEpKKux4FXfRumUupCVTylWhami2TbgcfkbowpBaYCnwDbgSxjzFYRmS0iY1zFhgI7RGQn0BZ42k/xKqXOU8VTqiVhTTS524BXE4cZY1YAK2rsm1Hl/XJguW9DU0r5UkXLvUiaEKM3VC1Pn1BVyiaaNokgNiqCs0Zb7nagyV0pG0mMi+K0idShkDagyV0pG0mMj6KgzKEtdxvQ5K6UjSTGRXGiNAJKCwMdivIzTe5K2UhifBQnSyMw2i1jeZrclbKRtnFRnDGRlBdpcrc6Te5K2Ui7+CjO0kRb7jagyV0pG2kbF0UhkYTpSkyWp8ldKRtJjI/irGlCWHkJlJcFOhzlR5rclbKRljGRFIc1cW7ocEhL0+SulI2EhQmOqKbODU3ulqbJXSmbiYquSO56U9XKNLkrZTNRTZs532jL3dI0uStlM82axgLocEiL0+SulM3ExsYBcPr0KQ8lVSjzKrmLyEgR2SEiu0XkMTfHLxaRbBHZKCKbRWS070NVSvlCnCu5Hz+hK2FamcfkLiLhwAJgFJAETBSRpBrFfoNzhaYBOJfh+z9fB6qU8o3mzeMBOH5SW+5W5k3LfTCw2xizxxhTDCwDbqhRxgBxrvfx1FhAWykVPFo1bw5AwamTAY5E+ZM3y+xdBHxXZTsPuLxGmQzgbyJyP9AUuNYn0SmlfK6Vq+V+ukBb7lbmTctd3OwzNbYnAkuMMR2A0cCbIlKrbhGZLCK5IpKbn5/f8GiVUhcs0vUQ05kzmtytzJvkngd0rLLdgdrdLncCWQDGmHVAFJBQsyJjzCJjTLIxJrl169bnF7FS6sI4YgAoOlMQ4ECUP3mT3HOAbiLSRUQicd4wfa9Gmf3ANQAi0gtnctemuVLBKKIJ5QglhacDHYnyI4/J3RhTCkwFPgG24xwVs1VEZovIGFexh4C7ReQrYCkwyRhTs+tGKRUMRCgJi6JMF+ywNG9uqGKMWQGsqLFvRpX324Af+TY0pZS/lIVHYc6epai0jCYR4YEOR/mBPqGqlA2ZiGiipZjDJ4sCHYryE03uStmQOKKJoojvTxYGOhTlJ5rclbKhsMgYoinm4AlN7lalyV0pG4qIako0RRzS5G5ZmtyVsqHwJjE0DSvWbhkL0+SulA2JI4Zm4SV8ry13y9LkrpQdOaJpKiXacrcwTe5K2VFEFNFSrC13C9PkrpQdOWJoYoo4fKqQ8nJ9mNyKNLkrZUeOaBzlhZSUGY6eLg50NMoPNLkrZUeOGMJNCeGUcUj73S1Jk7tSduSIBiAK7Xe3Kk3uStmRK7lHU8xBbblbkiZ3pezItWBHTFiJPqVqUZrclbIjRxQAFzU1OtbdojS5K2VHrpZ7+6ZGb6halFfJXURGisgOEdktIo+5Of6CiGxyvXaKyHHfh6qU8hlXn3v7pkZnhrQojysxiUg4sAD4Cc7FsnNE5D3X6ksAGGMerFL+fmCAH2JVSvmKq+XeNspwKE+TuxV503IfDOw2xuwxxhQDy4Ab6ik/Eec6qkqpYOVqubeOLuNUUSkFRaUBDkj5mjfJ/SLguyrbea59tYhIJ6ALsPLCQ1NK+Y2r5d4qsgxAx7pbkDfJXdzsq2syignAcmNMmduKRCaLSK6I5Obn53sbo1LK11wt95aR5QB6U9WCvEnueUDHKtsdgAN1lJ1APV0yxphFxphkY0xy69atvY9SKeVbEc6hkPEOZ3eMttytx5vkngN0E5EuIhKJM4G/V7OQiPQAWgDrfBuiUsrnXN0yceElADrW3YI8JndjTCkwFfgE2A5kGWO2ishsERlTpehEYJkxRucPVSrYRTQBBEd5Ic1jHNpytyCPQyEBjDErgBU19s2osZ3hu7CUUn4l4my9l5wlMS5KW+4WpE+oKmVXjmgoOUPbuCi9oWpBmtyVsitXy71dfJQ+pWpBmtyVsqm5/5hM9sZOtI2L4sipQkrKysnOhrlzAx2Z8gVN7krZVErX3aTPu4+ju+IxCO+uKCE9HVJSAh2Z8gVN7krZVFqfHWRNeYYXf9Oa42u6M3mSg6wsSEsLdGTKFzS5K2VX+w+QdnQ+t52Zx4kvujEs4s+kDRPIyAh0ZMoHvBoKqZSyoJ59yT6eyp9jH6TlgG/4eNsEslfeqi13i9CWu1I2lb2jP+mLZ5OVJaT0+pSr/nsX6emQnR3oyJQvaHJXyqZy9nYj69ZfkpYG3RLjOJ1wkKwsyMkJdGTKF7RbRimbmj5+HXy1CoBLh13JB5/t5Mofl5GWFh7YwJRPBFVyLykpIS8vj8JCfaDCiqKioujQoQMOhyPQoShwzgxZchaAS9s0wxj4Nr+A3u3jAxyY8oWgSu55eXnExsbSuXNnRNxNI69ClTGGo0ePkpeXR5cuXQIdjgLnE6plxVBWyqVtmgGw+7Amd6sIqj73wsJCWrVqpYndgkSEVq1a6W9lwcS1YAelZ+mcEEOYwLeHCwIbk/KZoErugCZ2C9PvNshUJPeSszSJCKdTq6bsztfkbhVBl9wDLTw8nP79+1e+9u3bx6pVq4iPj2fAgAH06NGD1NRUPvjgg8pzMjIyEBF2795due+FF15ARMjNzQ3EZSjlmWvBDkrOAHBJ62bs1pa7ZVgjufvwibro6Gg2bdpU+ercuTMAV111FRs3bmTHjh3Mnz+fqVOn8tlnn1We17dvX5YtW1a5vXz5cpKSknwWl1I+V6XlDs6bqnuPnKa0rDyAQSlf8Sq5i8hIEdkhIrtF5LE6yqSLyDYR2Soif/ZtmB7MmtWoH9e/f39mzJjBiy++WLlv7Nix/PWvfwVgz549xMfHo+vEqqBW2XI/l9xLygz7j50JYFDKVzwmdxEJBxYAo4AkYKKIJNUo0w14HPiRMaY38Es/xNoozp49W9klM27cuDrLDRw4kG+++aZyOy4ujo4dO7JlyxaWLl3K+PHjGyNcpc6fw7lIdtXkDmjXjEV403IfDOw2xuwxxhQDy4AbapS5G1hgjPkBwBhz2LdhupGR4VwqrOImXcX7C+yiqdot884779RZzt1SsRMmTGDZsmW8++679f5gUCoo1Gi5X9K6KYDeVLUIb5L7RcB3VbbzXPuq6g50F5HPRWS9iIx0V5GITBaRXBHJzc/PP7+IK2RkgDHOF5x730gz2m3cuJFevXpV2/fTn/6UN998k4svvpi4uLhGiUOp81bZ5+7shomNcpAYF6Utd4vw5iEmd+PXajZbI4BuwFCgA7BGRPoYY45XO8mYRcAigOTk5NpN3xCxefNmnnzySRYvXlxtf3R0NHPmzKF79+4BikypBqjRcgdn14yOdbcGb5J7HtCxynYH4ICbMuuNMSXAXhHZgTPZN84URDNn+v0j1qxZw4ABAzhz5gxt2rRh/vz5XHPNNbXKTZgwwe+xKOUTNVru4EzuyzfkYYzR5xJCnDfJPQfoJiJdgP8AE4BbapR5F5gILBGRBJzdNHt8GWi9fNgVU1BQu9UydOhQTpw4Uc/Hu//8VatW+SgqpfygxlBIgEvaNKOgqJTvTxbSLj46QIEpX/DY526MKQWmAp8A24EsY8xWEZktImNcxT4BjorINiAbeMQYc9RfQSulfKCiW6a0SrdMax0xYxVeTRxmjFkBrKixb0aV9wb4leullAoF4ZGA1OpzB2dyv6qbPqcRyqzxhKpSquFEnK33Ksk9oVkk8dEOdmnLPeRpclfKzhzR1W6oigiXttE5ZqxAk7tSdlaj5Q7OfncdDhn6NLkrZWc1Wu7g7Hc/erqYH04XBygo5Qua3Gs4dOgQt9xyC127dmXQoEFceeWV9U5D4M7QoUMDMtXvvn37EBGeeOKJyn1HjhzB4XAwderUes9dtWoVX3zxReX2pEmTWL58+XnHcqHnq0biiIaS6guoVN5U1WkIQpom9yqMMYwdO5bU1FT27NnDhg0bWLZsGXl5eYEOzWtdu3atNtf8W2+9Re/evT2eVzO5K5twxLhtuYMOhwx1mtyrWLlyJZGRkdxzzz2V+zp16sT9998POJcBvP322+nbty8DBgwgOzsbcM4kOWHCBPr168f48eM5e/ZcH2azZs0q3y9fvpxJkyYBzpbttGnTGDJkCF27dq1s5RpjeOSRR+jTpw99+/YlMzMTgPHjx7NixbnRqJMmTeLtt9+udQ3R0dH06tWr8jeHzMxM0tPTK4/n5+dz0003kZKSQkpKCp9//jn79u1j4cKFvPDCC/Tv3581a9YAsHr1aq/jM8YwdepUkpKSuO666zh82P9zxykfcETV6nO/qHk0UY4wTe4hLqgWyK5q1vtb2XbgpE/rTGofx8yf1t2K3bp1KwMHDqzz+IIFCwD4+uuv+eabbxg+fDg7d+7kpZdeIiYmhs2bN7N58+Z666jq4MGDrF27lm+++YYxY8Zw880385e//IVNmzbx1VdfceTIEVJSUkhNTWXChAlkZmYyevRoiouL+eyzz3jppZfc1lsxO2ViYiLh4eG0b9+eAwecM0Y88MADPPjgg/z4xz9m//79jBgxgu3bt3PPPffQrFkzHn74YQD+8Ic/NCi+devWsWPHDr7++msOHTpEUlISd9xxh1d/DyqAHDFw6lC1XWFhQtcEHTET6oI2uQeD++67j7Vr1xIZGUlOTg5r166tbMX37NmTTp06sXPnTlavXs20adMA6NevH/369fOq/rFjxxIWFkZSUhKHDjn/g61du5aJEycSHh5O27Ztufrqq8nJyWHUqFFMmzaNoqIiPv74Y1JTU4mOdv94+MiRI3niiSdo27ZtrXnlP/30U7Zt21a5ffLkSU6dOnXB8a1evbpyf/v27Rk2bJhXfwcqwNzcUAVn18yGf/8QgICUrwRtcq+vhe0vvXv3rtbVsWDBAo4cOUJycjLgfg73CnVNslR1f2Fh9RtXTZo0qXxfUXddnxEVFcXQoUP55JNPyMzMZOLEiXXGEhkZyaBBg/jd737H1q1bef/99yuPlZeXs27dujp/MJxvfKALYIckR3StbhlwJvf3vjrAmeJSYiKDNk2oemifexXDhg2jsLCwWnfHmTPnWjWpqan86U9/AmDnzp3s37+/csHsiv1btmxh8+bNlee0bduW7du3U15e7tWom9TUVDIzMykrKyM/P5/Vq1czePBgwNnd8tprr7FmzRpGjBhRbz0PPfQQc+bMoVWrVtX2Dx8+vNrygJs2bQIgNja2zha8N/GlpqaybNkyysrKOHjwYOX9CBXk3Ixzh3M3Vffkn27siJSPaHKvQkR49913+cc//kGXLl0YPHgwt912G3PmzAHg3nvvpaysjL59+zJ+/HiWLFlCkyZNmDJlCgUFBfTr14+5c+dWJmOAZ599luuvv55hw4bRrl07jzGMGzeOfv36cdlllzFs2DDmzp1LYmIi4EzMq1ev5tprryUyMrLeenr37s1tt91Wa//8+fPJzc2lX79+JCUlsXDhQsC50Mg777xT7YZqQ+IbN24c3bp1o2/fvkyZMoWrr77a47WqIOCIrjZxWAUdMRP6pL5fs/0pOTnZ1BwLvn379lqrGylr0e84yPxjLmQ/DU8chfBz3S/FpeX0mvExU66+hIdH9AhggKomEdlgjEn2VE5b7krZWYRrkewarffIiDA6tYzRlnsI0+SulJ25WbCjwiVtmulTqiFMk7tSdla5jqr74ZD7jpympKy8kYNSvuBVcheRkSKyQ0R2i8hjbo5PEpF8Ednket3l+1CVUj5XT8v90tbNKC03/Pto7cSvgp/H5C4i4cACYBSQBEwUkSQ3RTONMf1dr8U+jlMp5Q8eWu6gI2ZClTct98HAbmPMHmNMMbAMuMG/YSmlGkVly72w1qFLXMn9W+13D0neJPeLgO+qbOe59tV0k4hsFpHlItLRXUUiMllEckUkNz8//zzC9b/w8HD69+9f+Xr22Wd9VvemTZuqTf51IapOSFaXefPmVXsIa/To0Rw/ftwnn68sorLlXrtbplmTCNrHR2nLPUR581yxu2fKaw6Ofx9YaowpEpF7gNeBWpOLGGMWAYvAOc69gbE2iujo6MqnNn1t06ZN5ObmMnr0aL/UX9O8efO49dZbiYlx/gf21Q8WZSEO11BIN90y4Boxo8k9JHnTcs8DqrbEOwAHqhYwxhw1xhS5Nl8BBvkmvOBw4sQJevTowY4dOwCYOHEir7zyCgBTpkwhOTmZ3r17M3PmzMpzcnJyGDJkCJdddhmDBw/mxIkTzJgxg8zMTPr37185VW6FrVu3MnjwYPr370+/fv3YtWsXAM8//zx9+vShT58+zJs3r1Zsq1at4vrrr6/cnjp1KkuWLGH+/PkcOHCAtLQ00tLSAOjcuTNHjhyps959+/bRq1cv7r77bnr37s3w4cOrTV+sLKieljs4+92/zS+gvDwo22KqHt603HOAbiLSBfgPMAG4pWoBEWlnjDno2hwDbL/gyD56DL7/+oKrqSaxL4yqv5vl7Nmz9O/fv3L78ccfZ/z48bz44otMmjSJBx54gB9++IG7774bgKeffpqWLVtSVlbGNddcw+bNm+nZsyfjx48nMzOTlJQUTp48SUxMDLNnzyY3N7fa3C4VFi5cyAMPPMDPfvYziouLKSuVpzjrAAAOIklEQVQrY8OGDbz22mt8+eWXGGO4/PLLufrqqxkwYIDHS502bRrPP/882dnZJCQkVDtWV70tWrRg165dLF26lFdeeYX09HTefvttbr31Vm/+dlUoquxzd99yv7RNM84Ul3HgxFk6tIhpxMDUhfKY3I0xpSIyFfgECAdeNcZsFZHZQK4x5j1gmoiMAUqBY8AkP8bsV3V1y/zkJz/hrbfe4r777uOrr76q3J+VlcWiRYsoLS3l4MGDbNu2DRGhXbt2pKSkABAXF+fxc6+88kqefvpp8vLyuPHGG+nWrRtr165l3LhxNG3aFIAbb7yRNWvWeJXc61NXvWPGjKFLly6VP9wGDRrEvn37LuizVJCrZygkOIdDgnPEjCb30OLVXJ7GmBXAihr7ZlR5/zjwuE8j89DCbmzl5eVs376d6Ohojh07RocOHdi7dy/PPfccOTk5tGjRgkmTJlFYWIgxpsHT395yyy1cfvnlfPjhh4wYMYLFixfXO71uhYiICMrLzz1kUnNaYXfqq7fqNL/h4eHaLWN19QyFnDsXuveJBWD3m28z9KkpZGdDTg5Mn96YQarzoU+oeumFF16gV69eLF26lDvuuIOSkhJOnjxJ06ZNiY+P59ChQ3z00UeAcyGPAwcOkJOTA8CpU6coLS2td1rdPXv20LVrV6ZNm8aYMWPYvHkzqampvPvuu5w5c4bTp0/zzjvvcNVVV1U7r1OnTmzbto2ioiJOnDjBZ599Vnmsrs/zpl5lE+GRIGFQWrtRkJICd98WScShNnybu43sbEhPd+5XwU9n4a+hZp/7yJEjueOOO1i8eDH//Oc/iY2NJTU1laeeeopZs2YxYMAAevfuTdeuXfnRj34EOBfLyMzM5P777+fs2bNER0fz6aefkpaWxrPPPkv//v0r+/IrZGZm8sc//hGHw0FiYiIzZsygZcuWTJo0qXIK4bvuuqtWl0zHjh1JT0+nX79+dOvWrdrxyZMnM2rUKNq1a1dtfvWBAwe6rVe7YGxIpM453dPSICsLRo7pzztdSlmc7tx23Z9Xwc4YE5DXoEGDTE3btm2rtU9Zi37HQWhOV2Pe/2Xt/TNnGgNmRIflBoy5u9mLxoBzvwoYnPc6PeZY7ZZRyu7qaLmTkUH2SkNu4Y20S97M60xm5WflkJHR6CGqhtPkrpTd1bFIdkUf+1tZwv+WPEqL6zdw480GXUExNGhyV8ru6lgkOyfnXB/7+Bsup3dyIV3Gb+bLL3UK4FCgyV0pu6ujW2b69HM3TyNmZfDYqF78EP8fOgz9rlZZFXw0uStld3W03Gu6tlcbBnduybxPd1JQVNoIgakLocldKbvzMrmLCI+P7smRgmJeWb2nEQJTF0KTu1J2V8cNVXcGXNyC6/q1Y9HqPRw+6flpaBU4IZvc586l1l377Gzn/mD097//nUGDBtG3b18GDRrEypUrAx2SUk5ettwrTB/Rg9Lycl74dJcfg1IXKmSTe0qKc5hWRYIP9kejExISeP/99/n66695/fXX+fnPfx7okJRyqmucex06tWrKrVd0IjNnP7sOuZ9OQwVeyCb3ikej09Nhxgznn754NHrfvn306dOncvu5554jo46HNnJycujXrx9XXnkljzzySLXzahowYADt27cHoHfv3hQWFlJUVFRneaUaTQO6ZSrcP6wbTSMjmPPxN34KSl2okE3u4EzkU6bAk086/2zsOS9uv/12Fi5cyLp16wgPD/f6vLfffpsBAwZUm4FRqYBxxEB5CZR5PwKmZdNILv1+EB98XMb6PUcrn1oN5q5Ruwnp5J6dDS+9BE884fyzMZ+cO378OKdOnWLIkCGAc8peb2zdupVHH32Ul19+2Z/hKeW9ijndSxs2vfM9N7fg6HsDeWjeAcpnzQ76rlG78Sq5i8hIEdkhIrtF5LF6yt0sIkZEkn0XonsV/5CysmD27HNdNBea4L2dH914Mdd6TXl5eYwbN4433niDSy655LxjVMqnIirWUW1Ych/xk3B+/fwJ1i/uTkqPTxl3Uzl//HMZaWk1Bjxoqz4gPCZ3EQkHFgCjgCRgoogkuSkXC0wDvvR1kO5UfTQazvXBu6ZQP29t27bl8OHDHD16lKKiIj744AO35Vq0aEFsbCzr168HYNmyZfXWe/z4ca677jqeeeaZyqmBlQoK9SzYUa+MDJ6Ykkhqm7/xrx3DoNe3PLzibX71P68T3vYY6emueWhmzarVqveU/Os7Hqhz/V23z3maNhK4EvikyvbjwONuys0DrgdWAcme6g3mKX9///vfm0suucRce+215rbbbjMz65jidP369aZv377miiuuMI899pgZMmRInXU++eSTJiYmxlx22WWVr0OHDvnpCoJXsHzHqoqv3zZmZpwxhxr+3axcaUxCgjH/w5MmvkWZSZ+xy/Sd+bHp9OgHpssv/mmiY0vMxJZ/MM1blplXs06ZfUcKzJFThebjv5WahIRys3KlMQYq61m5snq97o7Xd8yf5/q7bm/h5ZS/Yjx0L4jIzcBIY8xdru2fA5cbY6ZWKTMA+I0x5iYRWQU8bIzJra/e5ORkk5tbvcj27dvp1auXtz+XAq6goIBmzZxrTD777LMcPHiQ3//+9wGOKriF2ndsCzs+hqXjoXmnc/3vXsjeOYj0V58m645fk1byBtmOX5D+6tMsu/3XXN7ibxQUlvDUF4+yYO293HvV//HA0OoLw3+5bzAPvP08Pxu4lD//ayLz/+shhnQ596v3+n0pTM36HT8b8Gf+tPEWXkx/iCtdx9ftrfuYp+MXcq6v6r6h89u8f/ju8xrhJyIbjDEeu769WYnJ3WKglT8RRCQMeAEvFsUWkcnAZICLL77Yi48Obh9++CHPPPMMpaWldOrUiSVLlgQ6JKUaruNguOwWKDndoNNyvkgl68GFpPUpgG09SUsqICt+ITnfpnLNjwvI2dKLzK9u5ZG+v+PVjXczLCWfy7tvobTcUHrwe1KafUl6UiYvrr2PySkL6RP3L06WtcLEx8OJE/SM28iNfZbz/z6fyl2DF9EjbhPHShMA6BG3qfqxeNexuDg4eZIe8TWPf3Xu3Piv3B+Li4WTp87/eAPqXvzPyTzBbNKGzYSZM/0zR76npj0eumWAeOAIsM/1KgQO4KFrJpi7ZWq69957q3WnXHbZZebVV191W/bjjz+uVXbs2LGNHHHwCtbvWPmWN10QFfueYJbb7on6jgfqXH/X7Q287JbxJrlHAHuALkAk8BXQu57yqzwldhNiyV35jn7H9jBnTpXE5bpntXKlc3/F+0D0bdupz93jaBljTCkwFfgE2A5kGWO2ishsERnjg18ean6er6tUQUK/W/uoOhd8RZdDWppzP9QY7TZzZq3RbvUdD9S5/q7b1zzeUPUXdzdU9+7dS2xsLK1atULEXVe/ClXGGI4ePcqpU6fo0qVLoMNRKmT58oZqo+nQoQN5eXnk5+cHOhTlB1FRUXTo0CHQYShlC0GV3B0Oh7bqlFLKB0J6bhmllFLuaXJXSikL0uSulFIWFLDRMiKSD/z7PE9PwPnglJ3oNduDXrM9XMg1dzLGtPZUKGDJ/UKISK43Q4GsRK/ZHvSa7aExrlm7ZZRSyoI0uSullAWFanJfFOgAAkCv2R70mu3B79cckn3uSiml6heqLXellFL1COrk7mlhbhFpIiKZruNfikjnxo/St7y45l+JyDYR2Swin4lIp0DE6UvBuAC7v3lzzSKS7vqut4rInxs7Rl/z4t/2xSKSLSIbXf++RwciTl8RkVdF5LCIbKnjuIjIfNffx2YRGejTALyZFzgQLyAc+Bboyrl55JNqlLkXWOh6PwHIDHTcjXDNaUCM6/0UO1yzq1wssBpYjxfrBQTzy8vvuRuwEWjh2m4T6Lgb4ZoXAVNc75OAfYGO+wKvORUYCGyp4/ho4COcq91dAXzpy88P5pb7YGC3MWaPMaYYWAbcUKPMDcDrrvfLgWsktOcK9njNxphsY0zFMvXrgVCfZtGb7xngSWAuzpW+Qp0313w3sMAY8wOAMeZwI8foa95cswHiXO/jca7oFrKMMauBY/UUuQF4wzitB5qLSDtffX4wJ/eLgO+qbOe59rktY5yLipwAWjVKdP7hzTVXdSfOn/yhzOM1uxZg72iM+aAxA/Mjb77n7kB3EflcRNaLyMhGi84/vLnmDOBWEckDVgD3N05oAdPQ/+8NElRT/tZQ78LcDSgTSry+HhG5FUgGrvZrRP7nswXYQ4g333MEzq6ZoTh/O1sjIn2MMcf9HJu/eHPNE4ElxpjficiVwJuuay73f3gB4df8Fcwt9zygY5XtDtT+Na2yjIhE4PxVrr5fg4KdN9eMiFwL/BoYY4wpaqTY/MXTNccCfYBVIrIPZ9/keyF+U9Xbf9t/NcaUGGP2AjtwJvtQ5c013wlkARhj1gFROOdgsSqv/r+fr2BO7jlANxHpIiKROG+YvlejzHvAba73NwMrjetORYjyeM2uLoqXcSb2UO+HBQ/XbIw5YYxJMMZ0NsZ0xnmfYYwxJtd9dSHBm3/b7+K8eY6IJODsptnTqFH6ljfXvB+4BkBEeuFM7lZelu094BeuUTNXACeMMQd9Vnug7yh7uNs8GtiJ8y77r137ZuP8zw3OL/8tYDfwT6BroGNuhGv+FDgEbHK93gt0zP6+5hplVxHio2W8/J4FeB7YBnwNTAh0zI1wzUnA5zhH0mwChgc65gu83qXAQaAEZyv9TuAe4J4q3/EC19/H177+d61PqCqllAUFc7eMUkqp86TJXSmlLEiTu1JKWZAmd6WUsiBN7kopZUGa3JVSyoI0uSullAVpcldKKQv6/5LQ261P2PtTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8e4c1c9310>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#lưu giá trị u_g vào file để so sánh ở bài 13\n",
    "np.savetxt('data/burgers_godunov_test1.dat', u_g)\n",
    "\n",
    "plt.plot(x, u_2, 'r+')\n",
    "plt.plot(x, u_g, x, u_exact)\n",
    "plt.plot(x, u_g_2, 'bx')\n",
    "\n",
    "plt.legend(['FDM', 'Godunov Method', 'Exact solution', 'u_g_2'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Với bài toán chúng ta đang xét có uL = 1.2 > uR = 0.4, sơ đồ Godunov suy biến về dạng sơ đồ sai phân bảo toàn (4). Nghiệm thu được từ hai sơ đồ là như nhau."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#xét bài toán sóng giãn\n",
    "uL = -0.4 # <0 - dịch chuyển sang trái\n",
    "uR = 1.2  # >0 - dịch chuyển sang phải \n",
    "\n",
    "u0 = np.zeros(nx)\n",
    "i = int(0.5/dx)\n",
    "u0[:i] = uL\n",
    "u0[i:] = uR\n",
    "\n",
    "u_4 = nonconservative_scheme(nt)\n",
    "u_5 = conservative_scheme(nt)\n",
    "u_6 = godunov_scheme(nt)\n",
    "\n",
    "#lưu giá trị u_g vào file để so sánh ở bài 13\n",
    "np.savetxt('data/burgers_godunov_test2.dat', u_6)\n",
    "\n",
    "u_exact = u0.copy()\n",
    "t = nt*dt\n",
    "i = int((0.5 + uL*t)/dx)\n",
    "j = int((0.5 + uR*t)/dx)\n",
    "for k in range(i+1, j): u_exact[k] = (x[k]-0.5)/t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3Xd8zVcfwPHPyU5kk9jEiBFZQoi9Wpuithq12xqtUlqtTWvUVmqrrbSqrQelRm0h9h5BrESW7HHvef4It4lMzZXBeT8vr+bee37nnPt7ku899/zO73uElBJFURTl7WKQ2x1QFEVRcp4K/oqiKG8hFfwVRVHeQir4K4qivIVU8FcURXkLqeCvKIryFlLBX1EU5S2kgr+iKMpbSAV/RVGUt5BRbncgPYUKFZJOTk653Q1FUZR85fTp00+llA6Zlcuzwd/JyQlfX9/c7oaiKEq+IoS4m5VyatpHURTlLaSCv6IoyltIBX9FUZS3UJ6d809LQkICAQEBxMbG5nZXlLeQmZkZJUqUwNjYOLe7oijZlq+Cf0BAAFZWVjg5OSGEyO3uKG8RKSXBwcEEBARQpkyZ3O6OomSbXqZ9hBArhRCBQoiL6bzeQwhx/vm/o0IIj//STmxsLAULFlSBX8lxQggKFiyovnUqbwx9zfmvBppn8PodoIGU0h2YDCz9rw2pwK/kFvW7p7xJ9DLtI6U8JIRwyuD1o8keHgdK6KNdRVHebPH+/oTv+B3Ig9vNSpDaRBITo4lPiCEhIZqEhBgSNLHEJcaSkJhAolZDolaD5sU/qUWj1aKVWqSUJP2Pf39+sa2uvT1dZu14rd3PjTn/fsD/0npBCDEQGAhQqlSpnOyTkkXTpk3jq6++0j2uXbs2R48ezeAI/Vq9ejW+vr4sXLgwx9pUcoc2Lo57gwaRcPce5OK3Lpn8g0em+E+ajMh+YA0oGpLNGjKXo8FfCNGIpOBfN63XpZRLeT4lVL169Tz4Uf/mkDJplGFg8Gozfy8H/5wM/MrbJXjpMhLu3qPUyhUUqF37tbalffaIO7d2c/7+Ic6FXuNcQii3jAyQyT50hJQU1ApsEg0xjTfDMNEctKZIrSkCc8yMrShgaoOVmS22BewpYG6FuYk5ZiZmmJuYU8DUHAszcyxNLbAws8DI2BBjA0OMDI0wMjBKeiyMMDIypOIr/l3+FzkW/IUQ7sByoIWUMjin2tUnf39/WrRoQd26dTl69CjFixfnt99+w9zcnLNnzzJ48GCio6MpV64cK1euxM7OjoYNG1KzZk32799PWFgYK1asoF69emg0GkaPHs3u3bsRQjBgwACGDh3Kvn37GDlyJImJiXh7e7N48WJMTU1xcnKid+/e/P777yQkJPDzzz9TqVIlDh48yPDhw4GkOelDhw5hZWXFzJkz2bJlC3FxcbRv356JEyfq+t+oUSOOHTtGu3btiIqKYsaMGUDSqPr06dMsWLCAdu3acf/+fWJjYxk+fDgDBw5kzJgxxMTE4OnpSZUqVVi/fj2WlpZERkbSpUsXevfuTcuWLQHo06cPbdq0oV27dowZM4YDBw4QFxfHJ598wqBBg1Kc16ioKDp37kxAQAAajYZvvvmGLl26cOrUKYYPH05UVBSmpqbs27cPgIcPH9K8eXNu3bpF+/btdf3fs2cP48ePJy4ujnLlyrFq1SosLS1xcnKie/fu7N+/n4SEBJYuXcqXX37JzZs3GTVqFIMHDwZI85wpuSPe35/gpUuxbtXqtQR+7bNH+J1ewvEHRzgf85gLhpIIw6SAay3B1awQlUVJImOLcCfEnjvhjkQkFkZjYkKJ0nZUq2SHewkbittaUMTGDGszo3x3TShHgr8QohTwC9BTSnldH3VO/P0Slx8+00dVOi7FrBnfpkqGZW7cuMHGjRtZtmwZnTt3Ztu2bXzwwQf06tWLBQsW0KBBA8aNG8fEiROZO3cuAImJiZw8eZKdO3cyceJE9u7dy9KlS7lz5w5+fn4YGRkREhJCbGwsffr0Yd++fVSoUIFevXqxePFiPv30UwAKFSrEmTNn+OGHH5g1axbLly9n1qxZLFq0iDp16hAZGYmZmRl79uzhxo0bnDx5Eiklbdu25dChQ5QqVYpr166xatUqfvjhB4KCgqhVq5YueG7evJmxY8cCsHLlSuzt7YmJicHb25v333+f7777joULF3L27NlU56Vr165s3ryZli1bEh8fz759+1i8eDErVqzAxsaGU6dOERcXR506dWjatGmK5ZK7du2iWLFi/PnnnwCEh4cTHx9Ply5d2Lx5M97e3jx79gxzc3MAzp49i5+fH6amplSsWJGhQ4dibm7OlClT2Lt3LwUKFGD69OnMnj2bcePGAVCyZEmOHTvGZ599Rp8+fThy5AixsbFUqVKFwYMHp3vO6tevn51fKeU/kFLyeNIkhJkZhceM1l/F0SHc9lvF7zd/5c/EYB4ZGWEgwdnMguZWZXAvWgMD09rsvWLCrouBxGu0FLUxo7qTPR3r2FGttB2VilhhZPhm3Burl+AvhNgINAQKCSECgPGAMYCUcgkwDigI/PD80zFRSlldH23ntDJlyuDp6QlAtWrV8Pf3Jzw8nLCwMBo0aABA79696dSpk+6YDh06pCgPsHfvXgYPHoyRUdL/Bfb29pw7d44yZcpQoUIFXT2LFi3SBf/k9fzyyy8A1KlThxEjRtCjRw86dOhAiRIl2LNnD3v27KFq1aoAREZGcuPGDUqVKkXp0qXx8fEBwMHBgbJly3L8+HGcnZ25du0aderUAWD+/Pn8+uuvANy/f58bN25QsGDBdM9LixYtGDZsGHFxcezatYv69etjbm7Onj17OH/+PFu3bgWSAvuNGzdSBH83NzdGjhzJ6NGjad26NfXq1ePChQsULVoUb29vAKytrXXlmzRpgo2NDQAuLi7cvXuXsLAwLl++rOt/fHw8tWrV0h3Ttm1bXVuRkZFYWVlhZWWFmZkZYWFh6Z4zFfxz3rM/dxJ19BiFv/kaI4dMk1NmLC6Cpxe2sPPKBv6ICeCKqQmGEmpZlmR4xU40rNKdhARjtp4OYMHhe9wKCsHKzIhuNUrSrWYpKhWxzryNfEpfq326ZfJ6f6C/Ptp6IbMR+utiamqq+9nQ0JCYmJgsH2NoaEhiYiKQNLp5+Wui7kr/K9QzZswYWrVqxc6dO/Hx8WHv3r1IKfnyyy9TTa/4+/tToECBFM916dKFLVu2UKlSJdq3b48QggMHDrB3716OHTuGhYUFDRs2zHR9u5mZGQ0bNmT37t1s3ryZbt266d7TggULaNasWbrHVqhQgdOnT7Nz506+/PJLmjZtSrt27dL9Gv3y/weJiYlIKXn33XfZuHFjhscYGBikON7AwEB3fFrnTMlZmmfPeDL9O8xcXbHr2vW/VxT+gBN/f82qJ0c4ZmaCVghcLB0ZXbY1zd37UciiEBcCwvn6lxv8eeER8YlaqpayZWZHd1q7F8PcxFB/byqPejO+v+QyGxsb7Ozs+OeffwBYu3at7ltAepo2bcqSJUt0QTwkJIRKlSrh7+/PzZs3s1zPrVu3cHNzY/To0VSvXp2rV6/SrFkzVq5cSWRkJAAPHjwgMDAwzeM7dOjA9u3b2bhxI126dAGSRud2dnZYWFhw9epVjh8/ritvbGxMQkJCmnV17dqVVatW8c8//+iCfbNmzVi8eLHumOvXrxMVFZXiuIcPH2JhYcEHH3zAyJEjOXPmDJUqVeLhw4ecOnUKgIiICN25SouPjw9HjhzRnbvo6GiuX8/6DOOrnDPl9QmaOw9NcAhFJkxAGP6HABz+gHPb+9N/UyP6P/PlhoUV/Zxa8VvbX9nc/R8+8BmNgbRmzLbztF10mL8uP6FL9ZLsHFaPXz+uQ6fqJd+KwA/5LL1DXrZmzRrdBd+yZcuyatWqDMv379+f69ev4+7ujrGxMQMGDGDIkCGsWrWKTp066S74vrgYmZ65c+eyf/9+DA0NcXFxoUWLFpiamnLlyhXdtIelpSXr1q3DMI0/Jjs7O1xcXLh8+TI1atQAoHnz5ixZsgR3d3cqVqyomyYCGDhwIO7u7nh5ebF+/foUdTVt2pRevXrRtm1bTExMdO/T398fLy8vpJQ4ODiwffv2FMdduHCBUaNGYWBggLGxMYsXL8bExITNmzczdOhQYmJiMDc3Z+/evemeBwcHB1avXk23bt2Ii4sDYMqUKboptMw0bdo0zXPm6OiYpeOV7Iu5cIHQjRux69EDc9dX/GYf/oDrByax4OHfHLAww97ckjFVPqSj50BMDZO+6Wm0ko0n7zFz9zUi4xLpV6cMw99xxsrs7czVJDKbasgt1atXly9v5nLlyhUqV66cSz1SFPU7+LpIjQb/Tp1JDAqi7M4/MbSyytqB4Q+4d3AKiwL+4n8WZlgaGPNhxe708PoEC2MLXTG/e6GM++0SFx6E41PWnknvuVKhcBbbyGeEEKezck1VjfwVRcl1oRs2Env5MsXnzM5a4NdqiDy6gLnnFrK1gBkmBQrQz/l9+lQbjo2pja5YSFQ8M3ZdZdOp+zhamTK/W1XauBfNd8syXwcV/BVFyVUJTwIJmjuXAnXqYNU8oxRhzwVd5+TvA/haBvLE0pzOTq0YVGMkhcwLpSj2vwuPGPPLBaLiEhlYvyzDmjhjaapC3gvqTCiKkquefPctMiGBIuO+yXhErkkk9sgc5p1bzDrrApQ2deSnxvPwcPRMWUwrmbn7GksO3sKjZNIKnjd1iic7VPBXFCXXRP5zmIj/7aLQ0CGYlC6dfsEnl7iwYyBfEYy/dQG6l2vHpz5fYW5knqJYaFQ8Qzf6cfjmU7rXLMX4Ni6YGr0dq3delQr+iqLkCm1sLI8nT8bEyYmCAwakXUiTQMKhWfx4/keW21hSyKQgyxp9j09Rn1RFLz4IZ9Da0wRFxDH9fTe6eKvkkBlRwV9RlFwRvHQpCffuUWrVSgyeLw1O4dkjbm/pxhjtI67YWtG2dDNG1x6HtUnqu25/ORPAl79cwM7ChC2Da+FZ0jYH3kH+pm7yegusXr2ahw8f6h7379+fy5cv51j7/v7+uLq65lh7St4Xd/sOwcuWY92mDQWSpeHQuXecQ6sb0d0wiMcFbJnbcC5TG85KFfgTNFom7LjEiC3n8Chpy+9D66rAn0Vq5J/PaDSaNG/Wysjq1atxdXWlWLFiACxfvvx1dE1RsiRF4rbRX7z8IvLUCn46MpnZdtZUsHZiQbNlFClQJFU9z2ITGLDGlxN3QuhbpwxftqyE8RuSdC0nqDP1in766Sfc3d3x8PCgZ8+eANy9e5cmTZrg7u5OkyZNuHfvHpCU1njYsGHUrl2bsmXL6pKbPXr0iPr16+Pp6Ymrq6suLcSePXuoVasWXl5edOrUSZdqwMnJiUmTJlG3bl1mzJihuxMXkkbV7u7uAEyaNAlvb29cXV0ZOHAgUkq2bt2Kr68vPXr0wNPTk5iYGBo2bIivry+LFy/miy/+/eNbvXo1Q4cOBWDdunXUqFEDT09PBg0ahEajSXUuxowZg4uLC+7u7owcORKAJ0+e0L59ezw8PPDw8NDl+9doNAwYMIAqVarQtGlTXU6kW7du0bx5c6pVq0a9evW4evWq7tx99NFHNGrUiLJly3Lw4EH69u1L5cqV6dOnj64P6Z0zJe969scfRB8/juOIzzAqlGx5ZmIcCTuGMP7kVGbZW9O4RAPWtNmcZuAPjYqnx7ITnL4byuzOHoxr46IC/6t6salHXvtXrVo1+bLLly//+2DnaClXttTvv52jU7WZ3MWLF2WFChVkUFCQlFLK4OBgKaWUrVu3lqtXr5ZSSrlixQr53nvvSSml7N27t+zYsaPUaDTy0qVLsly5clJKKWfNmiWnTJkipZQyMTFRPnv2TAYFBcl69erJyMhIKaWU3333nZw4caKUUsrSpUvL6dOn6/rh4eEhb926pSs3efLkFP2RUsoPPvhA7tixQ0opZYMGDeSpU6d0r714HBgYqOuTlFI2b95c/vPPP/Ly5cuydevWMj4+Xkop5UcffSTXrFmT4lwEBwfLChUqSK1WK6WUMjQ0VEopZefOneWcOXN07y0sLEzeuXNHGhoaSj8/PymllJ06dZJr166VUkrZuHFjef36dSmllMePH5eNGjXSnbsuXbpIrVYrt2/fLq2srOT58+elRqORXl5e0s/PL8Nz9rqk+B1UXlliWJi8VruOvN2ps9QmJv77QvgDGby0oey1pLx0Xe0q55+eJzVaTZp1BD6LlU1nH5TOY3fKfVce51DP8w/AV2Yhxqppn1fw999/07FjRwo9H63Y29sDcOzYMV2K5Z49e6YYTbdr1w4DAwNcXFx48uQJAN7e3vTt25eEhATatWuHp6cnBw8ezDAl8YukawCdO3dmy5YtjBkzhs2bN7N582YA9u/fz4wZM4iOjiYkJIQqVarQpk2bdN9PeimdFy1axOnTp3XplGNiYlLluLG2tsbMzIz+/fvTqlUrWrdurTtHP/30E5CUcdPGxobQ0NA0U2FHRkZy9OjRFOmvX+TlAWjTpg1CCNzc3ChcuDBubm4AVKlSBX9/fwICAjI8Z0reEzh3LprQUEotW/pv4ra7x7ixrTdDbYwIMrXku7pTaFW2VZrHPwqPoceyEzwKj2VVH2/qlC+UZjklc/k3+Lf4LseblGmkYU5L8jLJ0wfL53mU6tevz6FDh/jzzz/p2bMno0aNws7OLsOUxMlTMXfp0oVOnTrRoUMHhBA4OzsTGxvLxx9/jK+vLyVLlmTChAmZpmF+UdfLKZ2llPTu3Ztvv/023eOMjIw4efIk+/btY9OmTSxcuJC///473fJppcLWarXY2tqmuTlM8mPSS8NsaGiY4TlT8paY8+cJ27QZu54fYObikvTk+Z85tPtTvnCwx9zMllVNFuHu4J7m8feCo+m+/Djh0Qn81K8G3k72Odj7N49eJsmEECuFEIFCiIvpvC6EEPOFEDeFEOeFEF76aDenNWnShC1bthAcnLQLZUhI0ibLtWvXZtOmTQCsX7+eunXT3KJY5+7duzg6OjJgwAD69evHmTNnXiklcbly5TA0NGTy5Mm6bwQvAn2hQoWIjIzUXV8AsLKyIiIiIs260krp3KRJE7Zu3apLaRwSEsLdu3dTHBcZGUl4eDgtW7Zk7ty5ugDepEkTFi9eDCTN8z97lv5ua9bW1pQpU4aff/4ZSPpwPHfuXAZnLqXspnFWco5MTOTRhAkYOTjgMGxY0pMnfuTXPcMZ6mBPKTtnNrb5Od3AfzMwks4/HiMyLpH1A2qqwK8H+rpCshrIKClHC8D5+b+BwGI9tZujqlSpwtixY2nQoAEeHh6MGDECSNr1atWqVbi7u7N27VrmzZuXYT0HDhzA09OTqlWrsm3bNoYPH54iJbG7uzs+Pj66i59p6dKlC+vWraNz584A2NraMmDAANzc3GjXrp1uygaSLp4OHjxYd8E3uRcpne/evau7kOzi4sKUKVNo2rQp7u7uvPvuuzx69CjFcREREbRu3Rp3d3caNGjAnDlzAJg3bx779+/Hzc2NatWqcenSpQzPxfr161mxYgUeHh5UqVKF3377LcPyyb3qOVNyT+iGDcRdvkLhr77EsEAB+Hsq649MYpxDQWoV82F1y3VpXtgFuPLoGV2XHiNRq2XTQB/cS6ilnPqgt5TOQggn4A8pZaoF3UKIH4EDUsqNzx9fAxpKKR+9XPYFldJZyYvU7+CrS3jyhNstWmJerRoll/yA2DmSFTe3MtfejsYlGzGzwSxMDNO4yQu4EBDOBytOYG5syPoBNSnnYJnDvc9/8lpK5+LA/WSPA54/l27wVxTlzfDk2++QGg1FvvwCtn7IwkcH+dHejhZOLZhabyrGBmlvpnIzMIJeK09gZWbExgE+lLS3SLOc8t/k1MLYtK6SpvrKIYQYKITwFUL4BgUF5UC3FEV5nSIPHSJi1y4K9e+D8T+fMevJIX60s6GDcwe+rfdtuoH/QVgMPVecxNDAgPX9a6rA/xrkVPAPAEome1wCePhyISnlUilldSlldQcHhxzqmqIor0NS4rYpmJQuha3hL0x5dpGfbKzpXqk742uNx9Ag7TvVgyPj6LniBJFxifzUtwalCxZIs5ySPTkV/HcAvZ6v+vEBwjOa71cUJf97+uOPJNy/j0PVQMbxkC3WlvRz7ceYGmMwEGmHnsi4RD5cfYoHoTGs6O2NS7HUSdwU/dDLnL8QYiPQECgkhAgAxgPGAFLKJcBOoCVwE4gGPtRHu4qi5E1xt28TvGw5VhWNmFQ6it0W5gzxHMJA94Hp3isTl6hh4E++XHr4jKU9q1GjjFrO+TrpJfhLKbtl8roEPtFHW4qi5G1SSh5/8zUGhomsqpfIbgszRlQbwYeu6Y/5NFrJ8I1nOXormNmdPWhSuXAO9vjtpDIh5aD8ktrYycmJp0+f5nY3lHzq2dYNRJ/240ytRLY4mPGxx8cZBn4pJWN/vcCuS48Z19qFDl4lcrC3by8V/BVF0RvNozs8mTaVsMJavvUxpZdLLwZ7DM7wmJm7r7Hp1H2GNCpP37plcqinigr+r2jy5MlUqlSJd999l27dujFr1iwAzp49i4+PD+7u7rRv357Q0FAATp8+jYeHB7Vq1WLRokW6elavXs2QIUN0j1u3bs2BAwcAsLS0ZOzYsXh4eODj46NLCJdW6ujw8HCcnJzQarVAUoqDkiVLkpCQkKLfP//8M66urnh4eFC/fn0gKf3CyJEjcXNzw93dnQULFujKL1iwAC8vL9zc3HR3zUZFRdG3b1+8vb2pWrWq7m7c1atX065dO9q0aUOZMmVYuHAhs2fPpmrVqvj4+OjSYKSXvll5Q0SHEDi0PYmxWqa2MqZDhY6MrD4yw3xYa4/f5YcDt+hesxSfN62Qg51V8m1it+knp3M1RL/Bo5J9JUbXGJ3u676+vmzbtg0/Pz8SExPx8vKiWrVqAPTq1YsFCxbQoEEDxo0bx8SJE5k7dy4ffvih7vlRo0ZlqR9RUVH4+PgwdepUvvjiC5YtW8bXX3/NkCFD6NWrF71792blypUMGzaM7du34+HhwcGDB2nUqBG///47zZo1w9g45frpSZMmsXv3booXL05YWBgAS5cu5c6dO/j5+WFkZKQL0pCUI+jMmTP88MMPzJo1i+XLlzN16lQaN27MypUrCQsLo0aNGrzzzjsAXLx4ET8/P2JjYylfvjzTp0/Hz8+Pzz77jJ9++olPP/2UgQMHsmTJEpydnTlx4gQff/xxhsnglHwkOoSY6S0IvRjLn94GuNRsyTc+32QY+A9eD2LCjks0qeTI5Pdcs5Q0UdEfNfJ/BYcPH+a9997D3NwcKysrXbrk8PBwwsLCaNCgAQC9e/fm0KFDqZ5/sflLZkxMTHQpkl+kP4ak1NHdu3fX1XX48GEgKc/Pi7TOmzZtSpH++YU6derQp08fli1bptuYZe/evQwePBgjo6QxwIsU1ZCU8O3l9vfs2cN3332Hp6cnDRs2JDY2VrdxTaNGjbCyssLBwQEbGxvduXFzc0uVvvnFBjEv5wtS8qmYUOTqtlz7K5QQK8Gjrg2YWm9quuv4Aa4/iWDI+jM4O1oyr1tVDA1U4M9p+Xbkn9EI/XV51TxIGaWANjIy0k3VACnSLxsbG+uOMzQ0JDExMc06XpRp27YtX375JSEhIZw+fZrGjRunKrtkyRJOnDjBn3/+iaenJ2fPns2wfy9SKCdvX0rJtm3bqFixYoqyJ06cSJVyOXk65sTExEzTNyv5VFwkrO/EiRMPsAkxY1ffSnzbfG66d+4CPI2Mo+/qU5iZGLKyjzeWpvk2DOVrauT/CurWrcvvv/9ObGwskZGR/PnnnwDY2NhgZ2en245x7dq1NGjQAFtbW2xsbHQj9PXr1+vqcnJy4uzZs2i1Wu7fv8/JkyczbT+91NGWlpbUqFGD4cOH07p16zT3+L116xY1a9Zk0qRJFCpUiPv379O0aVOWLFmiC+7Jp33S0qxZMxYsWKD7EPTz88u0zy9kN32zkgclxsHmDzhx7xImfmbcqGzN0E/XYWpomu4hsQkaBq09TVBEHMt7VaeYrXkOdlhJTn3kvgJvb2/atm2Lh4cHpUuXpnr16tjY2ACwZs0aBg8eTHR0NGXLlmXVqlUArFq1ir59+2JhYUGzZs10ddWpU4cyZcrg5uaGq6srXl6Zb3Ewf/58+vbty8yZM3FwcNC1Af9u8PLiovHLRo0axY0bN5BS0qRJEzw8PHB1deX69eu4u7tjbGzMgAEDUlyEftk333zDp59+iru7O1JKnJyc+OOPP7Jy6oCkD6yPPvqIKVOmkJCQQNeuXfHw8Mjy8UoeotXAtv5cCzjMtfNF8NRqqTNrNQVM0k/FIKVk9LbznL4byg89vPAoqVIz5ya9pXTWt7ya0jkyMhJLS0uio6OpX78+S5cuzVLgVt4MeeF3MNdJCTuG8vj8BqZqnPh4Uwxmn/SnzNDPMzxs3t4bzNl7nVHNKvJJo/I51Nm3T15L6fzGGDhwIJcvXyY2NpbevXurwK+8XaSEv77h2bl1DCtdhSGLw8CpBE6DhmZ42I5zD5mz9zrve5Xg44blcqizSkZU8H9FGzZsyO0uKEruOTyHhKML+KyCB1X3huAQpqXU/KkIk7Q3YwE4cy+UkT+fo4aTPdM6qCWdeYW64KsoStb4rkTum8g456o8eBTCeyfBpl07Cjzf/jMtj8NjGfjTaYramLGkZzVMjdJf/qnkLDXyVxQlcxe3wR8jWFC2Kn8kPGXpkaIYFojA8Yv0b1yMTdAwaN1pYuIT2TigJvYF0v92oOQ8NfJXFCVjt/6GXwaxpbQby2QwnwdVw/ZyAI6fj8DIPu20y1JKxv12kXP3w/i+syfOha1yuNNKZlTwVxQlfY/OweaeHCxclqmGz3jH1ofav97E3NMT244d0z1s3fG7bPENYGjj8jR3LZKDHVaySgX/POTAgQMcPXo02/VkNXX0tGnTUjyuXbt2tttW3iCh/rCuI1cs7RhVQEsl+0qMOOWIJjycIhMnIAzSDh8n74Qw8ffLNK7kyGfvqGRteZVegr8QorkQ4poQ4qYQYkwar5cSQuwXQvgJIc7rsxDkAAAgAElEQVQLIVrqo903jb6Cf1a9HPxzsm0lj4sKhrUdeCoTGFq4ENamtswpOJjIrb9g36sXZi+l+HjhUXgMH68/TUl7C+Z08cRA5ezJs7Id/IUQhsAioAXgAnQTQri8VOxrYIuUsirQFfghu+3mlnXr1lGjRg1dcjKNRsPdu3dxdnbm6dOnaLVa6tWrx549ewBo164d1apVo0qVKixdulRXz65du/Dy8sLDw4MmTZrg7+/PkiVLmDNnDp6enrpUES8cPHgQT09PPD09qVq1KhEREUgpGTVqFK6urri5uemSuyWXXuroMWPGEBMTg6enJz169ACS0kQA6dZ74MABGjZsSMeOHalUqRI9evR45XxHSj4QHwUbOhP37AHDy7vzLDGaBfXnEjt9PkZFiuAwJO1N+WITNAxee5qYeA1Le1bDxjz9/D5K7tPHap8awE0p5W0AIcQm4D3gcrIyEnixE7MN8DC7jT6eNo24K/pN6WxauRJFvvoq3devXLnC5s2bOXLkCMbGxnz88cesX7+eXr16MXr0aAYPHkzNmjVxcXGhadOmAKxcuRJ7e3tiYmLw9vbm/fffR6vVMmDAAA4dOkSZMmUICQnB3t6ewYMHY2lpyciRI1O1PWvWLBYtWkSdOnWIjIzEzMyMX375hbNnz3Lu3DmePn2Kt7e3Lld/Zr777jsWLlyYZqK1jOr18/Pj0qVLFCtWjDp16nDkyBFdjiHlDaBJhK19kQ/PMK5aK84Hn2NOwzk4/s+XwKtXKb5gPgYFUqdwkFLyzfaLnAsI58ee1dQF3nxAH9M+xYH7yR4HPH8uuQnAB883d98JZHw7YB61b98+Tp8+jbe3N56enuzbt4/bt28D0L9/fyIiIliyZIlugxdIysfzYlOW+/fvc+PGDY4fP079+vUpUyZp1yL7dFZMJFenTh1GjBjB/PnzCQsLw8jIiMOHD9OtWzcMDQ0pXLgwDRo04NSpU9l+nxnVW6NGDUqUKIGBgQGenp66dM/KG0BK+ONTuL6L5d6d2Bl8jqFVh9LA1JWn8xdg2bAhVs/3b3jZ2uN3+fl0AMMal6dZFXWBNz/Qx8g/rUm9l+cCugGrpZTfCyFqAWuFEK5SSm3yQkKIgcBAgFKlSmXYaEYj9NdFSknv3r359ttvU70WHR1NQEAAkJT/x8rKigMHDrB3716OHTuGhYWFLgd+RqmU0zNmzBhatWrFzp078fHxYe/evVmacskodXR6Mqo3eermjNJNK/nQgW/Bby17q3dlftBRWpZpyQC3ATwYNhyp1VL467Fp/t6evBPCpN8v06SSI5+qC7z5hj5G/gFAyWSPS5B6WqcfsAVASnkMMAMKvVyRlHKplLK6lLK6g4ODHrqmX02aNGHr1q0EBgYCSSmQ7969C8Do0aPp0aMHkyZNYsCAAUDSJi92dnZYWFhw9epVjh8/DkCtWrU4ePAgd+7c0dUDYGVlRURERJpt37p1Czc3N0aPHk316tW5evUq9evXZ/PmzWg0GoKCgjh06BA1XrrbMqPU0cbGxqm2ewSyVK/yhvFdBQenc8WtHV+Fn8G9kDsTa08k8uBBIv76i0Iff4xJidQbqwdGxPLJhjNJF3i7qgu8+Yk+Rv6nAGchRBngAUkXdLu/VOYe0ARYLYSoTFLwD9JD2znKxcWFKVOm0LRpU7RaLcbGxixatAh/f39OnTrFkSNHMDQ0ZNu2baxatYru3buzZMkS3N3dqVixIj4+PgA4ODiwdOlSOnTogFarxdHRkb/++os2bdrQsWNHfvvtNxYsWEC9evV0bc+dO5f9+/djaGiIi4sLLVq0wMTEhGPHjuHh4YEQghkzZlCkSJEUUzEZpY4eOHAg7u7ueHl5pdhroH379mnWq/bcfUNd3wN/juBpuYYM1T7A2sSGeY3nYZIgCZg8BZPy5SjYp3eqwxI1WoZs8CMiNoG1/WpgbaYu8OYneknp/Hzp5lzAEFgppZwqhJgE+Eopdzxf/bMMsCRpSugLKeWejOrMqymdlbfbG/c7+NAPVrUirlA5+hYvzo3w26xpvobKBSsTOHsOwUuXUnrtT1h4e6c69Nv/XeHHg7eZ3dmDDl6pvxUouSNHUzpLKXeSdCE3+XPjkv18Gaijj7YURdGTUH9Y3xlpUZAJztU4f38fcxrOoXLBysTdvEnwqlXYtG+fZuDffekxPx68TfeapVTgz6fUHb6K8jaKDoH1nUATx7o6ffjj/j4+9vyYd0q/g5SSxxMmYmBhgeOo1MuO/Z9GMXLLOdxL2DCu9cu39Cj5Rb4L/uqmIiW3vDG/ewmxsKkHhPpzstkEvr+6lsYlGzPIfRAA4dt/I9rXF8eRn6dK3BYTr2HwutMYGAgWdffCzFilaM6v8lXwNzMzIzg4+M35I1TyDSklwcHBmJmZ5XZXskerhe2D4d5RHrb4lpHXVlPKuhRT607FQBigCQsjcMYMzKtWxfb991McKqXkm98ucu1JBHO7elLS3iKX3oSiD/kqn3+JEiUICAggKCjfLRRS3gBmZmaUSGO5Y76ydzxc+pXYxt/w6aM9JGgTmNdoHpYmSak9Ar+fjebZM4pMGJ8qcdvmU/fZ+vxGrkYVHXOj94oe5avgb2xsrLsrVlGUV3RyGRydj6zej4kimKshV1nQeAFlbJL+pqLP+BH288/Y9+2bKnHbxQfhjNtxiXrOhRiubuR6I+SraR9FUf6jq3/C/76Aii1ZX8aTP27/wceeH9OgZAMAZEICjydMwKhoURw++TjFoeExCXy0/jQFC5gwr2tVDNWNXG+EfDXyVxTlP7h/Crb2g2JVOVn3I2btH0bjko0Z6D5QVyRk7Trirl+nxMIFKRK3SSn5Yus5HoXFsmVwLbUV4xtEBX9FeZMF34KNXcCqCA/bzmPk/k9SXOAFSHj0iKCFC7Fs2BDLJk1SHL7qiD+7Lz3h61aV8SpllxvvQHlN1LSPorypIoNgXdKKndhuG/j05ORUF3gBnkybBlothb/+OkXitrP3w/j2f1d4p3Jh+tVV19reNGrkryhvovjopBF/xCPo/QfTbmziSsiVFBd4ASL+3k/EX3tx+HwEJiX+zcQeHp3AkA1ncLQyY1Yn91fOQqvkfSr4K8qbRquBbf3gwRnoso5fYu/z681fGeQ+iIYlG/5bLDqaJ1OeJ27r/W/iNiklo7ae43F4LD8ProWthZrnfxOp4K8obxIpYecouLYTWs7ismNZpu7sSa2itfjI46MURZ8uXkLCw4eUXvsTwuTfAL/qiD97LifN81dV8/xvLBX8FeVNcmQu+K6A2sMI9+jMiD+6YG9uz/T60zE0+DcVQ9yNG0mJ2zp0SJG4Tc3zvz1U8FeUN8X5n2HvBHB9H22T8Yw98ClPop+wpvka7Mz+HcFLrZZHEyZiWKBAisRt4dEJfLI+aZ7/+04eap7/DadW+yjKm+DmPtj+EZSuC+0Ws+LSKg4GHOQL7y9wd3BPUTT81+3EnD6N46iRGNklfShIKRm59RyBEbEs6uGFjYXamOVNp4K/ouR3Aadhc09wqAhd13Ms8AwLzy6kZZmWdK3YNUXRxNBQAmfOxNzLC5sOHXTPrzziz1+XnzCmRWU8S9rm9DtQcoFegr8QorkQ4poQ4qYQYkw6ZToLIS4LIS4JITboo11FeesFXYf1HaFAIfhgG4+1sYw+NJqyNmUZX2t8qqmbwO+/RxMRQZHx/yZuOx8Qxnf/u8K7LoXpW8cpF96EkhuyPecvhDAEFgHvkrSZ+ykhxI7nu3e9KOMMfAnUkVKGCiFUSkBFya7wB7CuAxgYQs9fSbAoyOe7+xCniWN2w9lYGKdMuRx95gzhW7dh368vZhWTkrNFxCYwdKMfDpamzOyo1vO/TfQx8q8B3JRS3pZSxgObgPdeKjMAWCSlDAWQUgbqoV1FeXtFhyQF/pgw+GAbFCzHLN9ZnA86z+Q6k1PcyAXPE7eNn4BRsaI4fPJJ0nNSMvbXiwSExjC/W1W1nv8to4/gXxy4n+xxwPPnkqsAVBBCHBFCHBdCNE+rIiHEQCGErxDCV+XsV5R0xEfDhi4Qchu6bYCiHuzx38OGqxvo6dKTpk5NUx0S8tNPxN24QZGxYzGwSPpG8LNvADvOPeSzd5yp7mSf6hjlzaaP4J/W98SXt9oyApyBhkA3YLkQItVVJSnlUilldSlldQcHBz10TVHeMJoE+Lk3PPCF91dAmfrcf3af8UfH4+7gzmfVPkt1SMKDBwQtXIRl48ZYPU/cduNJBON2XKRO+YJ81LB8Tr8LJQ/QR/APAEome1wCeJhGmd+klAlSyjvANZI+DBRFySqtFn4bAjf2QKvZ4NKWeE08nx/8HANhwMz6MzE2SL1E8/HUaQAUGfsVALEJGoZs8KOAiRFzOnuq/PxvKX0E/1OAsxCijBDCBOgK7HipzHagEYAQohBJ00C39dC2orwdpIRdY+D8Jmj0NVT/EICZp2ZyJeQKU+tOpZhlsVSHRezbR+Tff+Mw5BOMiyfNxk764zLXnkQwu4snjtb5fE9i5T/LdvCXUiYCQ4DdwBVgi5TykhBikhCi7fNiu4FgIcRlYD8wSkoZnN22FeWtICXs+hJO/gi1hkD9pLtyd/vvZtO1TfR26Z0iYdsL2qgoHk+ZiqmzM/a9egHw5/lHbDhxj0ENytKggppafZvpJb2DlHInsPOl58Yl+1kCI57/UxQlq6SE3V/BicXg8wk0nQJCcO/ZPd08//Bqw9M8NOiHH0h89Iji69chjI25HxLNmG3n8Sxpy8imFdM8Rnl7qDt8FSWvkhJ2j4XjP0DNj6DZVBCCOE0cIw+OxFAYMqv+rDTn+WOvXSdkzU/YdHwfi2rVSNBoGbLRDwQs6FYVY0P1p/+2U4ndFCUvkhL2fA3HF0GNQdD8W3h+A9aLef4FjRdQ1LJo6kO1Wh5PnIihpSWOn38OwPd7rnPufhiLuntR0t4i1THK20d9/CtKXiMl/DUOji0E7wHQYrou8O/y38Xma5vpU6VPmvP8AOG//krMmTM4jhqFkZ0dh288ZcnBW3SrUYpW7qk/LJS3kwr+ipKXSAl7x8PR+VC9H7ScqQv8957dY8LRCbg7uDPMa1iahyeGhhI4Yybm1ath074dwZFxfLblLOUdLRnX2iUn34mSx6lpH0XJK6SEfRPhyDyo3hdaztIF/nhNfKbz/ACBM2ehiYqi6PjxIASjtp4nPCaBNR/WwNzEMM1jlLeTCv6KkhckxsHvw+HcRqjWB1p+Dwb/fjGfe2YuV0KuMK/RvDTn+QGifX0J/+UXCg7oj6mzM6uO3OHvq4FMaOOCSzHrHHojSn6hgr+i5LaoYNjcA+4dg4ZfQoPRuhE/wKGAQ6y9vJauFbvSuFTjNKuQ8fE8njgR42LFKPTRR1x6GM63O6/SpJIjvWs75dAbUfITFfwVJTcFXYMNneHZo6RcPW4dU7wcGB3I14e/poJdBUZ6j0ynEghes4a4Gzcp8cMPxBqZMGzjSWwtjJmptmNU0qGCv6Lkllt/w5Y+YGQCff6Ekt4pXtZoNXz1z1fEamKZWX8mpoamaVYTH/CAp4t+wPKdJlg1bsSYbee5/TSKdf1qYl9ApWlW0qZW+yhKbji1HNZ1BJsSMODvVIEfYOXFlZx4fIIxNcZQ1rZsmtVIKXkyZQoYGFDkq6/YeeERm07dZ3CDctQpX+h1vwslH1Mjf0XJSVpN0l27JxaDczPouAJMrVIVOxt4lkVnF9HcqTnty7dPt7rIffuIPHAAxy++4Im5LWO2/YNHSVtGvFvhdb4L5Q2ggr+i5JTHF2HnKLh39HmenslJWzC+5Fn8M0YfGk2RAkUYV2tcunP22qgoHk+dhmmFClh3787g1afRSljQVaVvUDKngr+ivG4xYbB/GpxaBma20G4xeHZPs6iUkglHJxAYHciaFmuwMkn9reCFoEXPE7d9/z2Lj9zjlH8oc7t4UqqgSt+gZE4Ff0V5XbRaOLsO9k6EmJCkG7cajQWL9LdM3HpjK3/d/YtPvT7F3cE93XKx164RsmYNtp06caWgE/O2HqOdZzHaVX15B1VFSZsK/oryOgSchp0j4eEZKOmTlKahaPrBHOB22G1mnJxBraK1+ND1w3TLSa2WxxMmYmhtjcWQYXy6xo+iNmZMaueq73ehvMH0MjEohGguhLgmhLgphBiTQbmOQggphKiuj3YVJU/RaiHAF377BJY3hmcPocMy6Lsr08CfoElgzD9jMDMyY2rdqRiI9P80w7ZtI8bPD8cvvmDiwQAehMYwt4sn1mZpp3xQlLRke+QvhDAEFgHvkrRX7ykhxA4p5eWXylkBw4AT2W1TUfKM2PCk9frX98DNvyAqCAyMofYwaPBFmit50rLg7AJd+gYHi/R32EoMCSFw1vdYeHtzqIw3v2w6y/AmzlR3Sn8qSVHSoo9pnxrATSnlbQAhxCbgPeDyS+UmAzOA9G9TVJT8IOw+XN4O13cnpWTQJiZdyHV+N2n5ZvkmGc7rv+zko5OsvriajhU6ppu+4YXAmbPQRkXBZ18wdvtFvErZMrRx+ey+I+UtpI/gXxy4n+xxAFAzeQEhRFWgpJTyDyGECv5K/ra6FYTdBUcXqD00KeCX8AbDV/9zCo8L56vDX1HaujSjqo/KsGz0qVOE//or9gMGMPRkBFLC3C5VMVLLOpX/QB/BP61FyFL3ohAGwBygT6YVCTEQGAhQqlQpPXRNUfRMSnj2IGkj9WZTs1mVZPLxyQTHBLOu5TosjNNfoinj43k0YSLGxYuz1aUpJw/eZXZnD7WsU/nP9DFkCABKJntcAniY7LEV4AocEEL4Az7AjrQu+kopl0opq0spqzs4pD/vqSi5Ju5Z0jSPVZFsV/X77d/Z7b+bT6p+QpVCVTIsG7x6DfG3bhE1+FNm/3OPth7FaK+WdSrZoI+R/ynAWQhRBngAdAV0d7BIKcMBXZIRIcQBYKSU0lcPbStKzooOSfqvefYusAZEBDDtxDSqFa7Gh1XSX9YJEB8QwNMffsC8yTsMulOAItaSye1cVbZOJVuyPfKXUiYCQ4DdwBVgi5TykhBikhCibXbrV5Q85UXwtyj4n6tI1Cby5T9fYoAB0+pOwzCNFA8vSCl5MnkKwsCAVVXbExAazdyuntiYq2WdSvbo5SYvKeVOYOdLz41Lp2xDfbSpKLkiJvvBf/mF5ZwNOst39b6jmGWxDMtG7N1L5MGDhPT6iDU3YxjauDzealmnogfqDl9FeRXRwUn/fYWlnMmdDzrPknNLaFmmJa3KtsqwrCYyiidTp2HoXIGh8c64l7BiWBPn/9SuorxMrRFTlFehm/Z59eAfkxjDV4e/wtHCkbE+YzMt/3ThQhKfPGFFjc5Ea2BOF0+VrVPRGzXyV5RXER0MwgBMbV750Pln5nP32V1WNF2BtUnGG6rHXr1KyNq1PKnfnM1Rtkxp50I5B8v/2mtFSUUNIxTlVUQHJ630MXi1P52Tj06y7so6elTuQY2iNTIsK7VaHo+fAFbWfG5Tm8aVHOlRU933ouiXCv6K8ipiQl75Ym9UQhTjjo6jtHVphnsNz7R82M9biTl3jg1e7cDKhunvu6tlnYreqWkfRXkV0SGvPN8/y3cWj6Iesab5GsyNzDMsmxgcTODs2QSVd2WtlQvL33fHwSrtjdsVJTvUyF9RXkX0q438Dz84zNbrW+lTpQ+ejp6Zlg+cMRNNVBRfl25Bt5qlecelcHZ6qyjpUsFfUV5FdDCY22WpaHhcOOOPjKe8bXk+8fwk0/JRJ04S/ttv7HRpgmGZsnzTunJ2e6so6VLTPoqSVVImBf8sjvynn5xOSGwIC5oswMTQJOOq4+N5PHEiz2wdWFGmERu7eGJhov48lddHjfwVJaviI0GbkKXgv+/uPn6//TsD3QfiUtAl0/LBK1cRf/s2syq1ZfC7LniWtNVHjxUlXSr4K0pWZfHu3pDYECYdn0Rl+8r0d++fabXx9+8TtHgxx0t6Eu9di08aldNHbxUlQ+p7paJkVRaSukkpmXJ8ChHxEaxougJjg4wTsEkpeTxpMnFSsMyjLWs7e6jNWZQcoX7LFCWrspDOec/dPfx19y+GVB1CebvMt1eM2PMXUf/8w+oKTRnUoZa6i1fJMSr4K0pW6aZ90h75P4t/xncnv8OloAu9XXpnWp0mMoqHk6dw27Y4wU3fo6dPaX32VlEypIK/omRVTMZJ3eadnkdIbAjja43PMEf/C4Hz56N5+pTl3p2Y3rkqBgbqLl4l56jgryhZ9SKpm1nqlThnA8+y5foWPqj8QZZW98Revkzo2nX8z8mHXn1bUcw24zt/FUXf9BL8hRDNhRDXhBA3hRBj0nh9hBDishDivBBinxBCfb9V8p/okKQbvF5K6pagSWDisYkULVA0SzdzSY2GO1+NI9zEgoD3+9DOU+3Fq+S8bAd/IYQhsAhoAbgA3YQQLw99/IDqUkp3YCswI7vtKkqOe5HR8yWrL63mZthNxtYci4WxRabVPN20Ga5eYmO1DnzTraZK2qbkCn2M/GsAN6WUt6WU8cAm4L3kBaSU+6WU0c8fHgdK6KFdRclZadzde+/ZPX48/yPvln6XBiUbZFpF4tOnPJr5PWcLlaPt5x9S0FIlbVNyhz6Cf3HgfrLHAc+fS08/4H96aFdRclZMaIrgL6Vk8vHJGBsYM6ZGqtnONF0ZPwXi47jdYwjvVCnyunqqKJnSR/BP6zurTLOgEB8A1YGZ6bw+UAjhK4TwDQoK0kPXFEWPooPB4t+kbn/c/oPjj44z3Gs4jhaOmR4e8s9RjPbtZpfruwz7sMnr7KmiZEofwT8AKJnscQng4cuFhBDvAGOBtlLKuLQqklIulVJWl1JWd3Bw0EPXFEVPpEyRzjksNoyZp2biXsidzhU7Z3q4Nj6em199wyOLgtQdPxIrs4zv/FWU100fwf8U4CyEKCOEMAG6AjuSFxBCVAV+JCnwB+qhTUXJWfFRoInTXfCdfXo2EfERjKs1DgOR+Z/RhVkLsQp6yJUug6hdRa3uUXJftoO/lDIRGALsBq4AW6SUl4QQk4QQbZ8XmwlYAj8LIc4KIXakU52i5E3J7u499fgUv978lV5VelHRvmKmh0bcuoNYv4pTTlXp/Wm319xRRckavSR2k1LuBHa+9Ny4ZD+/o492FCXXxPyb1G36yekUtyzOYI/BmR4mpcR3xFdYC0OcJ3yDpanKpajkDeoOX0XJiucj/1AjY66FXqNzxc6Z7scLcG7dNopcO8ul5t2o5aN25lLyDjUMUZSsiA4F4EJ80oeAWyG3LBwSTsycWTwuWJL2E4e/1u4pyqtSwV9RsuL5yP9iZAAGwoAqBatkesi+LyZTNvoZYspMrAqYve4eKsorUdM+ipIV0cGA4Hz4TcrZlss0jcOZvccp889Obvi8i0/LejnTR0V5BSr4K0pWxIQgze24GHwR90LuGReNjefJxAlEmlnSaMa4DMsqSm5RwV9RsiI6mPsF7AiPC8e1kGuGRX+bvBCnoLsYfPIZNo6Zb/auKLlBBX9FyYroEC5YFAAyvth79twtyu1Yy6Nybnj375pTvVOUV6aCv6JkRXQIF4wNMDcyp5xtuTSLJGq0XPp6EqaaRDxnT1WpmpU8TQV/RcmK6GAuiARcCrpgZJD2IrntK7bjdeMkUR26Y1/ROYc7qCivRgV/RcmMlCTEhHBFE5XulM/DwDDsls4l1NaRGl9/lsMdVJRXp9b5K0pmEqK5ZqAlAW26wX/X17OoFRmE2ZyFGJqpNf1K3qdG/oqSmegQLpiaAGlf7D3w92mqHf6NJ9XrUaaFytOv5A8q+CtKZqKDuWBqSiFjK4oUSLn7VlRcAkFTpqIxMqLmzEm51EFFeXUq+CtKZmKSRv5u1mVTreDZ9v1qXB9eQfT7CPOialtGJf9QwV9RMhH+7AH+Jsa4vZTP58rNh5TbvIynxcriOaRfLvVOUf4bFfwVJROXQq8C4FbYS/ecVis5NnYatnGRVJw+BWFomFvdU5T/RC/BXwjRXAhxTQhxUwgxJo3XTYUQm5+/fkII4aSPdhUlJ1x45o+QkipFa+ie+33LPmqc+5vQd9vi6F01F3unKP9NtoO/EMIQWAS0AFyAbkIIl5eK9QNCpZTlgTnA9Oy2qyg55UL0I8okarEytwPgaXg0hvNnEG1hjc/Ur3K5d4ry3+hj5F8DuCmlvC2ljAc2Ae+9VOY9YM3zn7cCTYS6913JB6SUXEgIxVX+e0vMb5MXUi7kPrajRmFkbZ2LvVOU/04fwb84cD/Z44Dnz6VZ5vmG7+FAqnSHQoiBQghfIYRvUFCQHrqmKNnzMOohISTibmgFQPyTQDx2beRheXecu3bI5d4pyn+nj+Cf1ghe/ocySCmXSimrSymrOzg46KFripI9F55eAMDNNGmscnfKNIy1iYQOGK4Styn5mj6CfwBQMtnjEsDD9MoIIYwAGyBED20rymt1IegCplLibFGMyCNHSPxrN5srNKa4a8Xc7pqiZIs+gv8pwFkIUUYIYQJ0BXa8VGYH0Pv5zx2Bv6WUqUb+ipLXXHx6kcrxiRga2/J40iRiCxfnZ+dGlHMskNtdU5RsyXbwfz6HPwTYDVwBtkgpLwkhJgkh2j4vtgIoKIS4CYwAUi0HVZS8JkGbwOXgS7jGxhB88D4Jd+9xrE1fzCzMcbA0ze3uKUq26CWrp5RyJ7DzpefGJfs5Fuikj7YUJafcCrtFrCYOr8BEgv93DutWrThmV46yZglqvl/J99QdvoqSjvNB50FKSh82RZgYU3jMaG4FRlLWwTK3u6Yo2aaCv6Kk48LTCzS9ZoLmsSkOH75PvI0dD8NjKVtIzfcr+Z8K/oqSjpv3z9Ftbzxm9vHYdWzPnadRwP/bu/foqOprgePfnceQByGQhBBF3g+V64NHRFjXolShXJcVXVVK121RW69dWl2XolVbrQi4hPqgWi+tCEVFqz2BF6AAAA5jSURBVBXbVQ2I9XVFvVXECNHrgwCJQnkmJCEhmZBMZvb9Yw7cGJLMJDN5zJn9WWuYM5nfnPPbmcnmd/Y5c37YyN+4giV/Y1pR21jL5IIS0mr95OVXIxm5lJYfT/428jexz5K/Ma3Y/o8NzNwaoOGCIaRm+SB1AKXldYjACCv7GBew5G9MC+r3ow/8gSN9YfjM0ZCSCYlJlB6uZXD/VFKS7fLNJvZZ8jemharnnqfvV2W88t1BDJAGSAte2qG0vM7q/cY1LPkb04zvUBnljzzCl6P6EJh+PngrIDULVaW0vNbO9DGuYcnfmGbKfrOMgK+R31/SxNkDz4H6SkjL5lBNA3WNfkbZwV7jEpb8jXHUvvc/1Gx8lZq5MzmUJZydczZ4g8m/tLwWsNM8jXtY8jcGCBw7xsElS0gePoxV55STmpTKGVlnBMs+aVmUHLbTPI27WPI3Bqh4YhW+PXt47XvD+LBiK/dMvQdPwA8+L6RlUVpeS5onkbx+KT3dVWOiwpK/iXsNpV9RsWoV5d8ax+Oe9/npOT/lspGXBUs+AKlZzpk+6XZBN+MalvxNXFNVDi5ejN+TxF0TdvCd4d/hpvE3BZ+sd5J/Wjalh2sZmWP1fuMelvxNXKvZsAHv5s08My3A0GHncN+/3keCOH8W3goAGvoMYG9VvdX7jatE5Xr+xsQif3U1B5Yu5evByWybOpA/fftRUpKa1fSd5L+vIQXVGjvTx7hKRCN/EckSkTdEZKdzP6CVNuNF5AMR+VxEPhWR70eyTWOiZf/yh/BXVbHm3zw8NmMFOak532zg1PxL64L/IdgXvIybRFr2uRN4S1XHAG/R+vSMXmCeqv4LMAt4RET6R7hdYyJS90kRR9f9hb/nJ3DL3OWMHTD25EZO8t9RHdxBtrKPcZNIk/9s4Gln+WngipYNVHWHqu50lvcDZcDACLdrTKdpUxP/e/vNVKXDqfNvZdpp01pvWF8JfTLZVdHAqZkppHmsSmrcI9LkP0hVDwA497ntNRaRyYAHKGnj+RtEpFBECsvLyyPsmjGt2/Dgz8jcXcGOeRcwd+J1bTds9gUvq/cbtwmZ/EXkTRH5rJXb7I5sSEROAZ4BrlPVQGttVPUJVc1X1fyBA23nwESXqrLyzfsZ/Py77DtrED/62R/aP2/fW4GmZVFaVmslH+M6IfdjVfWStp4TkUMicoqqHnCSe1kb7foBrwB3q+rmTvfWmE5SVR4ufJis/3qWZE3kgoefIikxxMffW0ljSg5HG5rsYK9xnUjLPgXANc7yNcDLLRuIiAf4G7BWVV+McHvGdFhAA9z/4f18sv5Jpm5X8m68iZRhw0O/0FtJjWQAdkE34z6RJv9lwAwR2QnMcB4jIvkistppMweYBlwrIkXObXyE2zUmLAENsPiDxfz1s+f5z03peEaMIPv668N7cX0lFYFg0reyj3GbiE5fUNUK4OJWfl4IXO8sPws8G8l2jOkMf8DPPe/fQ0FJAct2jSe97GPynnqMBI8n9It9x6CxlgO+NFKSEzg1M7XrO2xMN7LLOxhX8gV8/PK9X1JQUsBtA3/AyFc+pd/l3yV9yvnhrcC5rs+eY6mMyOlLQoJd0M24iyV/4zo1jTXMf3s+r379Kj+fOJ+L/lxMQmoqg26/PfyVOF/wKqn1WMnHuJIlf+MqxZXFzN0wl/f3vc+vp/ya7301EO+WLeQuWEBSTk7oFRznXNdn51EPo+xMH+NC9pVF4xobSjew6P1FZHgyeHLWk5zdZwQl111Kyrnn0H/O1R1bmVP2qdAMO9PHuJIlfxPzfH4fDxU+xHPbn2PSoEk8dOFD5KTmcGDhvfiPHGHoH1cjCR3cyXVG/lWaYWUf40qW/E1MK/OWceumWykqL2LeuHnMnzSf5IRk6ouKOLJuHVnz5pFy5pkdX7G3CoAj9GWElX2MC1nyNzGr8GAht71zG94mLw9Oe5BZI2YBwQu3Hbh3EUm5ueTcckvnVu6toD4hnax+6WSkJEex18b0Dpb8Tcyp89Wx8tOVrP18LUMyhrB65mpGDxh94vnKZ5+lYft2Bj/6KIl9Ozlq91ZQTYZN3Whcy5K/iRmqysavNrK8cDll9WVcMfoK7jjvDvp6/j9B+w4e5PDvHiP9wmlkzJzR+W3VV3I4kG71fuNalvxNTCiuLGbplqV8fOhjxmWPY/n05Zw78NyT2h26fynq95N3993tX7EzBP/Rw5T7+9qZPsa1LPmbXq26oZoVRSt4ofgF+nn6sXDqQq4cfSWJCYknta195x2Ovv46A+fPxzNkSETbbaqroIphNvI3rmXJ3/RKXp+X9SXrWVG0gurGauaMncPNE24ms09mq+0D9fUcXHIfnlGjyP5xOxO0hCnxWCVVehb5VvM3LmXJ3/Qqu6p28eKOF1lfsp6jvqNMzJ3Ir87/Fadnnd7u6w4/vhLf3r0MXfs0Es6F29rT1EhyUx01ksHgAXZBN+NOlvxNj2v0N/LG7jdYV7yOrWVbSU5IZubwmcwZO4cJuRNC1u4bSkqoWLOGzNmzSZ88OfIOOd/uTUzPIdEu6GZcypK/6RFNgSa2V27n9d2v89LOl6hqqGJIxhAWTFrA7NGzyUrJCms9qsrBRYtJSEsj9/ZfRKdzzrd7U/rbVKLGvSJK/iKSBbwADAe+BuaoalUbbfsBXwJ/U9WbI9muiT2+gI8vKr6g8GAhHx36iKKyIup8dSRKItOHTOfq069myilTSJCOXYahpqAA75Yt5C1eRFJ2dlT62lRbQRLQLys3KuszpjeKdOR/J/CWqi4TkTudx3e00XYJ8E6E2zO9nKpScayCfbX72Hd0H7uP7mbboW0UlRdR31QPwKjMUVw28jLyB+VzXt55ZKd2Lmn7jxzh0G8eIHX8ePpfdVXUYjhctp88IHvgqVFbpzG9TaTJfzZwkbP8NLCJVpK/iEwCBgF/B/Ij3KbpBgEN4Av4aAo00eBvoK6xjlpfbfDWWPuN5fL68hPJfn/d/hNJ/rgxA8Zw5egryc/LZ9KgSWGXdEIp++0j+Kurybt3Yccv3NaOqvKD5AF5p1jyN+4VafIfpKoHAFT1gIictJ8sIgnAw8CPaGXKx2jbU7yVkh//e1dvJmZomD/XZvcaxjHOPs4tGxihwQ9SokKSs5zUbFn4kmDFby1lQFmHo2g9AK2BpDOEg+uiN+oHyA3UADB08GlRXa8xvUnI5C8ibwJ5rTx1V5jbuAnYqKr/DHXWhojcANwAMHTo0DBX/02JnhSqs+w4dnMtf+vi/CstfiYanN1HnOeO3xJUSAQSERLVuW+23N6Y2x+9ME6ipyVRf142JEd3TqJKYEf/sUzNsHP8jXuJaltjwzBeLFIMXOSM+k8BNqnq6S3a/An4FhAA+gIe4Peqemd7687Pz9fCwsJO980YY+KRiHysqiHL65EOkQuAa4Blzv3LLRuo6okajIhcC+SHSvzGGGO6VqT7y8uAGSKyE5jhPEZE8kVkdaSdM8YY0zUiKvt0JSv7GGNMx4Vb9onukTJjjDExwZK/McbEIUv+xhgThyz5G2NMHLLkb4wxcajXnu0jIuXA7ghWkQMcjlJ3YkW8xRxv8YLFHC8iiXmYqoa8HnmvTf6REpHCcE53cpN4izne4gWLOV50R8xW9jHGmDhkyd8YY+KQm5P/Ez3dgR4QbzHHW7xgMceLLo/ZtTV/Y4wxbXPzyN8YY0wbYjr5i8gsESkWkV3OHMItn+8jIi84z38oIsO7v5fRFUbMC0TkCxH5VETeEpFhPdHPaAoVc7N2V4mIikjMnxkSTswiMsd5rz8Xkee6u4/RFsZne6iIvC0i25zP96U90c9oEZE1IlImIp+18byIyO+c38enIjIxqh1Q1Zi8AYlACTCS4AQxnwDjWrS5CXjcWZ4LvNDT/e6GmKcDac7yjfEQs9MuA3gX2Exwzoge73sXv89jgG3AAOdxbk/3uxtifgK40VkeB3zd0/2OMOZpwETgszaevxR4leCEelOAD6O5/Vge+U8Gdqlqqao2An8mOKF8c7MJTiwP8BfgYgk1l2TvFjJmVX1bVb3Ow81ArE9EG877DLAEeAA41p2d6yLhxPwfwApVrQJQ1ahMjdyDwolZgX7Ociawvxv7F3Wq+i7BWUPbMhtYq0Gbgf7OjIlREcvJfzDwz2aP9zo/a7WNqjYB1QTnHI9V4cTc3E8IjhxiWciYRWQCMERVN3Rnx7pQOO/zWGCsiPxDRDaLyKxu613XCCfme4EfisheYCNwS/d0rcd09O+9Q2J5pvPWRvAtT10Kp00sCTseEfkhkA9c2KU96nrtxiwiCcBvgWu7q0PdIJz3OYlg6ecignt374nIWap6pIv71lXCifkHwFOq+rCITAWecWIOdH33ekSX5q9YHvnvBYY0e3waJ+8GnmgjIkkEdxXb283q7cKJGRG5BLgLuFxVG7qpb10lVMwZwFnAJhH5mmBttCDGD/qG+9l+WVV9qvoVUEzwP4NYFU7MPwHWAajqB0AKwWvguFVYf++dFcvJ/yNgjIiMEBEPwQO6BS3aHJ9gHuAq4L/VOZISo0LG7JRAVhJM/LFeB4YQMatqtarmqOpwVR1O8DjH5aoay3OAhvPZfongwX1EJIdgGai0W3sZXeHEvAe4GEBEziSY/Mu7tZfdqwCY55z1MwWoVtUD0Vp5zJZ9VLVJRG4GXiN4psAaVf1cRBYDhapaAPyR4K7hLoIj/rk91+PIhRnzg0Bf4EXn2PYeVb28xzodoTBjdpUwY34NmCkiXwB+4BeqWtFzvY5MmDHfCqwSkZ8TLH9cG8uDORF5nmDZLsc5jrEQSAZQ1ccJHte4FNgFeIHrorr9GP7dGWOM6aRYLvsYY4zpJEv+xhgThyz5G2NMHLLkb4wxcciSvzHGxCFL/sYYE4cs+RtjTByy5G+MMXHo/wBHmqkBW0TEsgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8e4c1aa190>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, u_4, x, u_5, x, u_6, x, u_exact)\n",
    "plt.legend(['nonconservative scheme','conservative scheme', 'godunov scheme', 'exact solution'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kết quả ba phương pháp là hoàn toàn khác biệt. Phương pháp Godunov cho kết quả chính xác hơn.\n",
    "\n",
    "Sơ đồ Godunov ta đang dùng là sơ đồ bậc 1, nó có tính chất tiêu tán, tính nhớt như đã phân tích ở phần 1.\n",
    "\n",
    "### [Bài 13. Định lý Godunov, sơ đồ TVD, sơ đồ MUSCL](Bai_13.ipynb)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}