{ "cells": [ { "cell_type": "markdown", "id": "03617dd0-ce8d-4c5f-a9e5-edb8395c21b2", "metadata": {}, "source": [ "# Nonlinear heat PDE\n", "\n", "Diffrax can also be used to solve some PDEs.\n", "\n", "(Specifically, the scope of Diffrax is \"any numerical method which iterates over timesteps\". This means that e.g. semidiscretised evolution equations are in-scope, but e.g. finite volume methods for elliptic equations are out-of-scope.)\n", "\n", "---\n", "\n", "In this example, we solve the nonlinear heat equation\n", "\n", "$$ \\frac{\\partial y}{\\partial t}(t, x) = (1 - y(t, x)) \\Delta y(t, x) \\qquad\\text{in}\\qquad t \\in [0, 40], x \\in [-1, 1]$$\n", "\n", "subject to the initial condition\n", "\n", "$$ y(0, x) = x^2, $$\n", "\n", "and Dirichlet boundary conditions\n", "\n", "$$ y(t, -1) = 1,\\qquad y(t, 1) = 1. $$\n", "\n", "---\n", "\n", "We spatially discretise $x \\in [-1, 1]$ into points $-1 = x_0 < x_1 < \\cdots < x_{n-1} = 1$, with equal spacing $\\delta x = x_{i+1} - x_i$. The solution is then discretised into $y(t, x_i) \\approx y_i(t)$, and the Laplacian discretised into $\\Delta y(t,x_i) \\approx \\frac{y_{i+1}(t) - 2y_{i}(t) + y_{i-1}(t)}{\\delta x^2}$.\n", "\n", "In doing so we reduce to a system of ODEs\n", "\n", "$$ \\frac{\\mathrm{d}y_i}{\\mathrm{d}t}(t) = (1 - y_i(t)) \\frac{y_{i+1}(t) - 2y_{i}(t) + y_{i-1}(t)}{\\delta x^2} \\qquad\\text{for}\\qquad i \\in \\{1, ..., n-2\\},$$\n", "\n", "subject to the initial condition\n", "\n", "$$ y_i(0) = {x_i}^2, $$\n", "\n", "for which the Dirichlet boundary conditions become\n", "\n", "$$ \\frac{\\mathrm{d}y_0}{\\mathrm{d}t}(t) = 0,\\qquad \\frac{\\mathrm{d}y_{n-1}}{\\mathrm{d}t}(t) = 0. $$\n", "\n", "---\n", "\n", "This example is available as a Jupyter notebook [here](https://github.com/patrick-kidger/diffrax/blob/main/examples/nonlinear_heat_pde.ipynb).\n", "\n", "\n", "!!! danger \"Advanced example\"\n", "\n", " This is an advanced example, as it involves defining a custom solver." ] }, { "cell_type": "code", "execution_count": 1, "id": "0a89f429-bab4-4a0f-800c-a0c8e1c7bf9b", "metadata": { "tags": [] }, "outputs": [], "source": [ "from typing import Callable\n", "\n", "import diffrax\n", "import equinox as eqx # https://github.com/patrick-kidger/equinox\n", "import jax\n", "import jax.lax as lax\n", "import jax.numpy as jnp\n", "import matplotlib.pyplot as plt\n", "from jaxtyping import Array, Float # https://github.com/google/jaxtyping\n", "\n", "\n", "jax.config.update(\"jax_enable_x64\", True)" ] }, { "cell_type": "code", "execution_count": 2, "id": "16da14af-420a-4d25-aa06-515a9baa50c2", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Represents the interval [x0, x_final] discretised into n equally-spaced points.\n", "class SpatialDiscretisation(eqx.Module):\n", " x0: float = eqx.field(static=True)\n", " x_final: float = eqx.field(static=True)\n", " vals: Float[Array, \"n\"]\n", "\n", " @classmethod\n", " def discretise_fn(cls, x0: float, x_final: float, n: int, fn: Callable):\n", " if n < 2:\n", " raise ValueError(\"Must discretise [x0, x_final] into at least two points\")\n", " vals = jax.vmap(fn)(jnp.linspace(x0, x_final, n))\n", " return cls(x0, x_final, vals)\n", "\n", " @property\n", " def δx(self):\n", " return (self.x_final - self.x0) / (len(self.vals) - 1)\n", "\n", " def binop(self, other, fn):\n", " if isinstance(other, SpatialDiscretisation):\n", " if self.x0 != other.x0 or self.x_final != other.x_final:\n", " raise ValueError(\"Mismatched spatial discretisations\")\n", " other = other.vals\n", " return SpatialDiscretisation(self.x0, self.x_final, fn(self.vals, other))\n", "\n", " def __add__(self, other):\n", " return self.binop(other, lambda x, y: x + y)\n", "\n", " def __mul__(self, other):\n", " return self.binop(other, lambda x, y: x * y)\n", "\n", " def __radd__(self, other):\n", " return self.binop(other, lambda x, y: y + x)\n", "\n", " def __rmul__(self, other):\n", " return self.binop(other, lambda x, y: y * x)\n", "\n", " def __sub__(self, other):\n", " return self.binop(other, lambda x, y: x - y)\n", "\n", " def __rsub__(self, other):\n", " return self.binop(other, lambda x, y: y - x)\n", "\n", "\n", "def laplacian(y: SpatialDiscretisation) -> SpatialDiscretisation:\n", " y_next = jnp.roll(y.vals, shift=1)\n", " y_prev = jnp.roll(y.vals, shift=-1)\n", " Δy = (y_next - 2 * y.vals + y_prev) / (y.δx**2)\n", " # Dirichlet boundary condition\n", " Δy = Δy.at[0].set(0)\n", " Δy = Δy.at[-1].set(0)\n", " return SpatialDiscretisation(y.x0, y.x_final, Δy)" ] }, { "cell_type": "markdown", "id": "7482e079-5ed1-4bc7-85a5-dcee3717ce7f", "metadata": {}, "source": [ "First let's try solving this semidiscretisation directly, as a system of ODEs." ] }, { "cell_type": "code", "execution_count": 3, "id": "d304a9d0-58c7-4d29-91e6-10bc65406b73", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Problem\n", "def vector_field(t, y, args):\n", " return (1 - y) * laplacian(y)\n", "\n", "\n", "term = diffrax.ODETerm(vector_field)\n", "ic = lambda x: x**2\n", "\n", "# Spatial discretisation\n", "x0 = -1\n", "x_final = 1\n", "n = 50\n", "y0 = SpatialDiscretisation.discretise_fn(x0, x_final, n, ic)\n", "\n", "# Temporal discretisation\n", "t0 = 0\n", "t_final = 1\n", "δt = 0.0001\n", "saveat = diffrax.SaveAt(ts=jnp.linspace(t0, t_final, 50))\n", "\n", "# Tolerances\n", "rtol = 1e-10\n", "atol = 1e-10\n", "stepsize_controller = diffrax.PIDController(\n", " pcoeff=0.3, icoeff=0.4, rtol=rtol, atol=atol, dtmax=0.001\n", ")" ] }, { "cell_type": "code", "execution_count": 4, "id": "d1ad0404-5a13-4506-bdab-cdcfaf5be609", "metadata": { "tags": [] }, "outputs": [], "source": [ "solver = diffrax.Tsit5()\n", "sol = diffrax.diffeqsolve(\n", " term,\n", " solver,\n", " t0,\n", " t_final,\n", " δt,\n", " y0,\n", " saveat=saveat,\n", " stepsize_controller=stepsize_controller,\n", " max_steps=None,\n", ")" ] }, { "cell_type": "code", "execution_count": 5, "id": "28185196-75f2-4465-ad59-ff45ec8c4d01", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAGiCAYAAACVh9NOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABESUlEQVR4nO3df3RU5Z0/8PedSTIh0oS6QCIYjb+VggTDkg3agqdZ08phyx/dpcACZgXXH9kDZKuAAhE5JdYfbGyLZsVNaffIgu3xx57C4qGpWdclhSWQU7WCi4ChrBOhfkkgQCaZ+3z/SBkdcz8PeSY3mcy97xdn/uCZ+2tm7nOf3Od+ns9jKaUUiIiIPCaQ7AMgIiIaCGzgiIjIk9jAERGRJ7GBIyIiT2IDR0REnsQGjoiIPIkNHBEReRIbOCIi8iQ2cERE5Els4IiIyJPYwBER0YB6++23MXPmTIwZMwaWZeH111+/5DoNDQ247bbbEAqFcP3112Pz5s3G+2UDR0REA6qjowMTJ07Exo0b+7T80aNHMWPGDNx5551obm7G0qVLsWjRIrz55ptG+7WYbJmIiAaLZVl47bXXMGvWLHGZ5cuXY/v27XjvvfdiZd/73vdw+vRp7Ny5s8/7SuvPgRIRUWq4cOECIpGIa9tTSsGyrLiyUCiEUCjU7203NjaitLQ0rqysrAxLly412g4bOCIij7tw4QKuuSYP4XCba9scPnw4zp49G1dWVVWFxx9/vN/bDofDyM3NjSvLzc1Fe3s7zp8/j2HDhvVpO2zgiIg8LhKJIBxuw5GP/wnZ2X1rHHTa28/j2quX4fjx48jOzo6Vu3H35iY2cEREPpGdPcyVBu7z7WXHNXBuycvLQ2tra1xZa2srsrOz+3z3BrCBIyLyDaW6oVS3K9sZSCUlJdixY0dc2a5du1BSUmK0HQ4TICLyCaWirr1MnD17Fs3NzWhubgbQMwygubkZLS0tAICVK1diwYIFseXvv/9+HDlyBI888ggOHjyI559/Hq+88gqWLVtmtF82cERENKD27duHSZMmYdKkSQCAyspKTJo0CWvWrAEAfPLJJ7HGDgCuueYabN++Hbt27cLEiRPx7LPP4qWXXkJZWZnRfjkOjojI49rb25GTk4PwqWddCzLJG/mPaGtrG5BncG7hMzgiIp9IlWdwbmEXJREReRLv4IiIfKInQMSNOzizIJNkYQNHROQTyu6Gsl1o4FzYxmBgFyUREXkS7+CIiPxCdfe83NhOCmADR0TkE4yiJCIi8gDewRER+YXdDdhd7mwnBbCBIyLyiZ4uyqAr20kF7KIkIiJP4h0cEZFf2N2A3f87OHZREhHR0OKzBo5dlERE5Em8gyMi8o2oS4O0mYuSiIiGEMvuhmX3v+POYhclERFR8vAOjojIL+xuwIU7uFQJMmEDR0TkFz5r4NhFSUREnsQ7OCIin7BUNyzlQpBJiqTqYgNHROQXtg3YLoT423b/tzEI2EVJRESexDs4IiKf6BkHZ7mynVTABo6IyC/sqEtRlKmRyYRdlERE5EmuN3Bvv/02Zs6ciTFjxsCyLLz++uuXXKehoQG33XYbQqEQrr/+emzevNntwyIiIrvbvVcKcL2B6+jowMSJE7Fx48Y+LX/06FHMmDEDd955J5qbm7F06VIsWrQIb775ptuHRkTka5Ydde2VClx/Bvftb38b3/72t/u8fG1tLa655ho8++yzAIBbbrkF77zzDv7pn/4JZWVlbh8eERH5RNKDTBobG1FaWhpXVlZWhqVLl4rrdHZ2orOzM/Z/27bx2Wef4c/+7M9gWf2PECIiSjalFM6cOYMxY8YgEHCps025FGSifHoHZyocDiM3NzeuLDc3F+3t7Th//jyGDRvWa53q6mqsXbt2sA6RiChpjh8/jiuvvNKVbVm27Ur3opUiA72T3sAlYuXKlaisrIz9v62tDVdddRWOtTyH7Oz4BvFC50nHbUTP/UHcfuCc8zqBc390LA+ePe28/Pmz4j6scx3Ob1w471x+7oJzecT5Ya8SFgcAFXH+C05FnE8H1eVcbncHxX3YwjqqW9h31HlbKiosr0k3pJTzXbySxv+4MC7okgLKsdiSyi3n8p73nC8uVlAqd76gWWnOywfS5QCCQJqwLWEdK0Mqly+QVqbwRoZwucoSVsjs/ccxAKisy8R928OGO5ZHh49wXj7rz4TyUeI+glnOjVVmKH6d9vbzKLhqCb7yla+I2yK9pDdweXl5aG1tjStrbW1Fdna2490bAIRCIYRCoV7l2dnDkJ2dFVeW0el88kfTeq9/USCY7lwecP66gnC+OAcCcgNgWcIFOiBdbIXyoHAx11yzlbAPsVzYh90lNzK20KWiglIDZ1iu6WYxb+AGYbRMQGh8EmnghMO1hNPNCgr7EGp/IF3+PgLpwudId/5urQyhPCSfoGIDJ62TJZQPE86dLM0fZsJ70eHCH3mXOV8r7Mvk60tQaJAzQ1mO5a4+drGj7vxB59cgE1MlJSXYsWNHXNmuXbtQUlKSpCMiIvKmnghINzKZpEYD5/qfrmfPnkVzczOam5sB9AwDaG5uRktLC4Ce7sUFCxbElr///vtx5MgRPPLIIzh48CCef/55vPLKK1i2bJnbh0ZERD7i+h3cvn37cOedd8b+f/FZ2cKFC7F582Z88sknscYOAK655hps374dy5Ytw3PPPYcrr7wSL730EocIEBG5jV2U/TN9+nQoJT8/cMpSMn36dBw4cMDtQ0k+3UkgRSGZRidJiycS5ORmcIb0nhAcIu1DetYmPWfTbUt61qbbllssad/SD6XrW5GeGQrP7SwpxNwW6qnmd5W+W8v0O9Sdn6bnrmldSpGL80BgFyUREZEHJD3IhIiIBgm7KImIyIssW7kySNuSureHGHZREhGRJ/EOjojIL+xoYgFoTttJAZ5v4JRKYN4i4ccTI4cG48c2jq40j3AUow9Noys17xnvI4F9m0ZLarflFinw0TS6Urct088nZMjQfR+m+xAj9rTfuRTdOQj5Dweh7id0TXKLcqmBS5Fky+yiJCIiT/L8HRwREfWwlG0+ZlHYTipgA0dE5Bc+ewbHLkoiIvIk3sEREfmFbbs00JtdlERENJSwgSPX6E4C4wSxUui02SEBmslCpZB1Mexe08PtVlLlBCYptcVZwAdhOIBARYUweiFBsu7ZgTSEQIkzoZolYdbOli587+JwB2l53YkrJhA3HD7gVkJzSlls4IiIfMKybVgutO9upPsaDGzgiIj8wrZdiqJMjQaOUZRERORJvIMjIvILn93BsYEjIvILNnA+YbuZhDmBH9ulaEkpY442YtA0EXICCXPdSqosRUtKkZK6bYnRo4bb0ZGiIkUB5x9Q9/mkd1RAiJZ0K0EyYHzuSNvSn5/C53ArulJDrMummTsSub6Q6/zbwBER+Y2Kyn8QGG2Hd3BERDSE+G2YAKMoiYjIk3gHR0TkFwwyISIiT2IDR65FQLkZXSkub1gOiJGMYoSjadSldh0hZ2F30LncNOpStw/Dz5cIBSGaUIhwFHNqCtGVgOazS9+hFNkplCfyu4qTaIrRseIuzC/AxnXGxYszoyWHNDZwRER+YSt3Gng3IjEHARs4IiK/sJVLXZSp0cAxipKIiDyJd3BERH7h2oSnvIMjIqKhxLbdexnauHEjCgoKkJmZieLiYuzdu1e7fE1NDW666SYMGzYM+fn5WLZsGS5cuGC0T97BucE0Tx1gnupGGc7o7WaEozTDs3ZGb+d9mM62LS6fwIzekkGZ6VvKGSpEXeo+gfSTB4JmeS2DUqSmbkZvcaZ2sxnZ9fkuxZ3L6zjv3Gx5ILG6TJe0bds2VFZWora2FsXFxaipqUFZWRkOHTqE0aNH91p+y5YtWLFiBerq6jB16lR8+OGHuOeee2BZFjZs2NDn/fIOjojIL2zl3svAhg0bsHjxYpSXl2PcuHGora1FVlYW6urqHJffvXs3br/9dsydOxcFBQW46667MGfOnEve9X0ZGzgiIr9QtnsvAO3t7XGvzs7OXruMRCJoampCaWlprCwQCKC0tBSNjY2Ohzl16lQ0NTXFGrQjR45gx44duPvuu40+Lhs4IiJKSH5+PnJycmKv6urqXsucOnUK0WgUubm5ceW5ubkIh8OO2507dy6eeOIJ3HHHHUhPT8d1112H6dOn49FHHzU6Pj6DIyLyC+XSOLg/PQ89fvw4srOzY8WhUMiFjQMNDQ1Yv349nn/+eRQXF+Pw4cNYsmQJ1q1bh9WrV/d5O2zgiIj8wuWB3tnZ2XENnJORI0ciGAyitbU1rry1tRV5eXmO66xevRrz58/HokWLAAATJkxAR0cH7rvvPjz22GMIBPrW+cgGzoBlGGGlXd6lGb3lGZY1J4D0nuHM3Xa3JtpO3IdZtJ20nURm9DYd/yNGj2pYwmRbUrQkhByVCc3oLeSWtCyzyEe7Ww4gCAZNZ303j4KVZ9WW6obZjN66eil9ctO6T5/LyMhAUVER6uvrMWvWLACAbduor69HRUWF4zrnzp3r1YgFgz25VpVBNC0bOCIiv0hSqq7KykosXLgQkydPxpQpU1BTU4OOjg6Ul5cDABYsWICxY8fGnuHNnDkTGzZswKRJk2JdlKtXr8bMmTNjDV1fsIEjIvKJLwRA9ns7JmbPno2TJ09izZo1CIfDKCwsxM6dO2OBJy0tLXF3bKtWrYJlWVi1ahVOnDiBUaNGYebMmfjBD35gtF82cERENOAqKirELsmGhoa4/6elpaGqqgpVVVX92icbOCIiv/DZbAJs4IiI/MKGSw2cC9sYBBzoTUREnsQ7OAfGIcFiWHMiyV6lbTkXJxISL4fkC+VRIWpJk6TYPKmy8z6k5bUJkqXPISUKdjHZslLOn8OSQviF31ubx9rwO7QMhyKISZghnwsq6FxnEkm2LH0+y3A4gEi3vOG2Um74gM/u4NjAERH5hYI82M90OymAXZRERORJvIMjIvIJZVtit7HZdlw4mEHABo6IyC989gyOXZRERORJ/r2Ds7sTWMfFP1sME8qKXQJSuS4yUIxMdCcRcs9xmSVJliMDDRNDA7BtIdLPONmyeVeOFC0pJVu2As7lAcjRebbwd2kg6HwyiNGS0rFqflclRX1KYZ/Sd6iNgpX27VzuWnSljnGkZgLXl8GgLOOk487b6f8mBoN/GzgiIp/x2zO4Aemi3LhxIwoKCpCZmYni4uLYtOOSmpoa3HTTTRg2bBjy8/OxbNkyXLhwYSAOjYiIfML1O7ht27ahsrIStbW1KC4uRk1NDcrKynDo0CGMHj261/JbtmzBihUrUFdXh6lTp+LDDz/EPffcA8uysGHDBrcPj4jIv2yXuij9ege3YcMGLF68GOXl5Rg3bhxqa2uRlZWFuro6x+V3796N22+/HXPnzkVBQQHuuusuzJkz55J3fUREZEhZ7r1SgKsNXCQSQVNTE0pLSz/fQSCA0tJSNDY2Oq4zdepUNDU1xRq0I0eOYMeOHbj77rvF/XR2dqK9vT3uRURE9EWudlGeOnUK0Wg0NondRbm5uTh48KDjOnPnzsWpU6dwxx13QCmF7u5u3H///Xj00UfF/VRXV2Pt2rVuHnq/SPkE3cx5J3YJGOaVBADVLUQZdhtGUWr+ihNzSAoReuLxCuVSpKRuW2Kkppt/jQrbCggRi9LvakP+fFKEpbKEz21JkatSRGQCeUylc0eYfVmly1GGYp5K026xBOqlWJc9gkEmg6yhoQHr16/H888/j/379+PVV1/F9u3bsW7dOnGdlStXoq2tLfY6fvz4IB4xEVGKsgPuvVKAq3dwI0eORDAYRGtra1x5a2sr8vLyHNdZvXo15s+fj0WLFgEAJkyYgI6ODtx333147LHH4qYxvygUCiEUCrl56ERE5DGuNsMZGRkoKipCfX19rMy2bdTX16OkpMRxnXPnzvVqxIJ/6tZQKkVGExIRpYKLUZRuvFKA68MEKisrsXDhQkyePBlTpkxBTU0NOjo6UF5eDgBYsGABxo4di+rqagDAzJkzsWHDBkyaNAnFxcU4fPgwVq9ejZkzZ8YaOiIi6j+lLFfmP0yVew/XG7jZs2fj5MmTWLNmDcLhMAoLC7Fz585Y4ElLS0vcHduqVatgWRZWrVqFEydOYNSoUZg5cyZ+8IMfuH1oRETkIwOSqquiogIVFRWO7zU0NMQfQFoaqqqqUFVVNRCHApVITjjlvI5rM30DYs5JOVrSuVjMG6jNJyjlBzTLJyhFYwLmOSeVsHwieSWjYqSmWY+8m7kopTPHEmbPDmpCBqUIS+nTKcMZvS1L03MizNwtnTtydKzu/BSOS6p/Yp1xL0eluG/hWqGT0DXJLXbApYHeqXELx1yUREQ+oWzzpOPO20mNBi41Yj2JiIgM8Q6OiMgvXJsux6dRlERENDS5F0WZGg0cuyiJiMiTeAfnBtPoSh1pgIlhrkYpKrFnF2YRb3ZUiGTU5iw0ndHb7JikSEndvsUIThf/GpW2JUVXStGEujNKirCUck7aUed9i7OPaxINSueCJcwmLv6umvNTzkUpRf+6GPDgZl0eitxKs5UiuSjZwBER+YR7yZbZRUlERJQ0vIMjIvIJvwWZsIEjIvILnz2DYxclERF5Eu/giIh8wm9BJv5t4BJIkmocQqxNtiy8JyVVlvq8pdB3XTeENLTAsFwK+de+Z7gtaTiAFK4OaIYDSJ9vEJ4nSCH5lpAI2XLxmIKW80klJlt2mGQ49l7Q7LuVzzXN+Wk4nMMyTaqcQLJl47qfyPVlEPjtGRy7KImIyJP8ewdHROQ3PgsyYQNHROQTfnsGxy5KIiLyJN7BERH5hN+CTNjAORCnp5dIUVmapLXilO/SKuIuDBPTatZRUlJlqVyX8FioALZtlrg5kUTPckJnKTpv4DsyLCGS0RISBQeE5MWA/J0EhEhN6TsPCBGc2t9VPEeEZMu2EE2YwPkpfIWaOiPUMW29NHu4ZHytSDbl0jO41JjQm12URETkTbyDIyLyCb8FmbCBIyLyCaXceX7m5hR8A4ldlERE5Em8gyMi8guXuih1QUJDCRs4E1LEVCI574T3xAAvMVpSiBjURakZRiYmkotS2oeYc9LwmLR5MIWoQduwUibSlSPlnATMIhl1pAhL6TuUWJbwPYmfAQhIuSilc0eMwNXkapSi/MS8ls6LWy7WS+NclEOUUgFXooZVivRRsouSiIg8iXdwRER+YVvudC+yi5KIiIYSv2UyYRclERF5Eu/giIh8ggO9fcLSRnENQsSUYS5KKT+gPEu1JsJRWqdbinw0zwcpRj+KOSqlmbuFciFSsuc9w0g/w3IdceZuMU+ktCWziEgdK+C8k0ACn1vMg9ktRAWnOdezhM5P4RwR80FKdcxNwr6115ckYhQlERGRB/j2Do6IyG/YRUlERJ7EKEoiIiIP4B0cEZFP+O0Ojg2cASm3nRzFpct5Z1pumOuvO4FZmU2jDDW5D8XoR8MclVK0ZFSXB9Pwc0iRnQkRtiXNti1FV7oZpWaac1JaHgCUkAfTOBdlAuennI9V2JBYLtdLqS6LeS1TjFIuPYNLkQaOXZRERORJvIMjIvIJv42DYwNHROQTfhsmwC5KIiLyJN7BERH5BKMoiYjIk9jAeY0ahKSn4jT3mtBi6SGt1LctnVDi8uZh9HJov1k5AEQHeDiAlJwZkMP+3UyqbEoZDhOQhhUkwoo6bysqDAeQjqlnW2bngpTQWXd+ulcHpDqmG74zCMMBBuOaRAD4DI6IyDeU/XmgSf9e5vveuHEjCgoKkJmZieLiYuzdu1e7/OnTp/HQQw/hiiuuQCgUwo033ogdO3YY7dP7d3BERAQgeV2U27ZtQ2VlJWpra1FcXIyamhqUlZXh0KFDGD16dK/lI5EI/vIv/xKjR4/GL3/5S4wdOxYff/wxRowYYbRfNnBERJSQ9vb2uP+HQiGEQqFey23YsAGLFy9GeXk5AKC2thbbt29HXV0dVqxY0Wv5uro6fPbZZ9i9ezfS09MBAAUFBcbHxy5KIiKfuDjQ240XAOTn5yMnJyf2qq6u7rXPSCSCpqYmlJaWxsoCgQBKS0vR2NjoeJz//u//jpKSEjz00EPIzc3F+PHjsX79ekSjZpNR8w6OiMgnbGW5knv14jaOHz+O7OzsWLnT3dupU6cQjUaRm5sbV56bm4uDBw86bv/IkSP4zW9+g3nz5mHHjh04fPgwHnzwQXR1daGqqqrPx+nfBk5KkAxd8mTDpMq2JhJOWEWM9OsWktYKEWS2LpmtuC1hHSnRs6aiSFGOcrnwOYTlo5ooSjla0mzfiTyrEKMiA1IUpZC8OIEoSmnfluX8uS3h/NR+bum7En4P6VzTnZ8BaR/SeSscr/T59PVSqstSEmbDa4XHZGdnxzVwbrFtG6NHj8aLL76IYDCIoqIinDhxAk8//TQbOCIicuBSqi5xaIaDkSNHIhgMorW1Na68tbUVeXl5jutcccUVSE9PRzD4+R81t9xyC8LhMCKRCDIyMvq07wF5BpeMcFAiItK7GEXpxquvMjIyUFRUhPr6+liZbduor69HSUmJ4zq33347Dh8+DPsLd9Qffvghrrjiij43bsAANHAXw0Grqqqwf/9+TJw4EWVlZfj0008dl78YDnrs2DH88pe/xKFDh7Bp0yaMHTvW7UMjIqIkqKysxKZNm/Czn/0MH3zwAR544AF0dHTEoioXLFiAlStXxpZ/4IEH8Nlnn2HJkiX48MMPsX37dqxfvx4PPfSQ0X5d76JMVjgoERHpJWsc3OzZs3Hy5EmsWbMG4XAYhYWF2LlzZyzwpKWlBYHA5/db+fn5ePPNN7Fs2TLceuutGDt2LJYsWYLly5cb7dfVBu5iOOgXW2KTcNA33ngDo0aNwty5c7F8+fK4/tcv6uzsRGdnZ+z/Xx6LQUREvSUzF2VFRQUqKioc32toaOhVVlJSgt/+9rfG+/kiVxu4wQoHra6uxtq1a9089L4xzVOnzXknlHdLUWpCuZDzUczbhwRyUYp5IoV9a7Yl5qg0jJaUjqnnPbOITHE7CVRiKYekdCoEAs77kKIudXQ5JE2Wjwbk7zYQdP4gVlSIBhVzUWrOT+EcEYMjhDoj1rFE8k0ORo5Kcl3SB3p/MRy0qKgIs2fPxmOPPYba2lpxnZUrV6KtrS32On78+CAeMRFRarJVwLVXKnD1Dm6wwkGldDBERCRTyqUZvVNkuhxXm+FkhoMSERF9kev3mckKByUiIr1kjINLJteHCSQrHJSIiPQ4o7cLkhEO6irTPHIJ5KKUJgwU80FK+RWliMEEIhyV6Yzeulm1hShHaZ1uMbrSbKZvQI68FHNUwr3KKv3kFoToSuGYggnMKCnNEy1Vcvk3kvctztwdMDx3NOenlItSrgPC7y3mj0wgF6W4vD9yTqYq5qIkIvIJt2cTGOrYwBER+YTfuihTYzADERGRId7BERH5hN/u4NjAERH5BJ/B+YRlSzFnunWkaEnDckCTJ89wNmPTvH2QI9ikGa/F/JG6SEYpF6VhbknT7QAQ0whJlXIw/hoVZ/oWoivF8yMBUWmm76g0+7jmdxXyVAbTnKMJxXNNc36K57RQB8TJNxPJRSm8J9Z9ge76Yp5llBLl2waOiMhvlHLnDzqVIq00GzgiIp/w2zM4RlESEZEn8Q6OiMgnlEtBJqlyB8cGjojIJ9hFSURE5AFGd3DTp09HYWEhampqBuhwhgZLSqAqlUuJcbVJXYVyw6TKYrkmma2coNksqbIU2q9fRxqKYDYcoFuzb+mvS9NhAol05QSkkHyhXAnltuXeX8iWJfwW0r61v6vziSsmYXbx/BSHFki/kzhMQFMvxbosJW5OrWTLfruDYxclEZFP+G2gd5+7KO+55x7853/+J5577jlYlgXLsnDs2LEBPDQiIqLE9fkO7rnnnsOHH36I8ePH44knngAAjBo1asAOjIiI3MUuSkFOTg4yMjKQlZWFvLy8gTwmIiIaAOyiJCIi8gDPB5ko5Zz0VPv3h/G09VLklbyKEqLOxKTKQgJhKRpNlwhZ2reYCFlaXpdsWYyWFJIqG0ZL6pItS+9J3SpudrdIP7kURSmVBwPuZVu2hO9cTACt+V2lRMzSdx6UInO156dhHRATkTsvb2mjm12q+xrSNWkwKFhQ+qtfn7eTCowauIyMDESjqRUWS0REPfz2DM6oi7KgoAB79uzBsWPHcOrUKXFMDBERUbIZNXDf//73EQwGMW7cOIwaNQotLS0DdVxEROSyi0EmbrxSgVEX5Y033ojGxsaBOhYiIhpA7KIkIiLyAM9HUYoSySEnRksK5bppb01zS3abRaPZmlx/Uq5BaR1xeU0ko7ROt7AP02hJXRSlmAdTiPwalFyUECIWkzg1ciAg5aKUn627du7ozk8xwlgol3JUivleE4iiNI6uHJrBeDZcGgeXIlGUvIMjIiJP8u8dHBGRz/jtGRwbOCIin7BhudK9yC5KIiKiJOIdHBGRX7jURSlOMjvEsIFzIM/obRhhpQu8MpzpWJwBWcq7KEWWQRd5aZZzUoqIBIBuIT+gLeQTFHNRSrkrhe3otiUFz7k5aDUqdN3IM307byeo6QIyvUBJ+5Z+I8uST9yAkCPTNOek7vwUz2njumEYXQkY1/FUm9GbswkQERF5AO/giIh8glGURETkSTb0PbQm20kF7KIkIiJP4h0cEZFPsIvSL2zNrLpCZJQcXSmF58m7kGYhNs1RKUap6WbbliITo86ngxjJqImilCIZu6RclMY5KuUKJkVYSpFfg1FZxdmzhXJdukQVcOd4xfyYuuhYMYrS+dyRJkhO05yfxpGXYv5WqY5pIh+FuizXfak8ebN269jKnQhI3aToQwm7KImIyJP8ewdHROQzChaUC2m23NjGYGADR0TkExzoTURE5AG8gyMi8omeIBN3tpMK2MAREfkEn8F5jZvhuoaJWJVu10Iftm2YUFZa3hbC63veM0uqLC2vTbYshvcL2zIcDtCl+XzSMIHBSLYsMU22bEtvAACkMHcpebKQbNkSlo/Kf54Hg877Nj13tOenS3UgKA0L0dRLyzSheiKG6BACL/J+A0dERAD8F2TCBo6IyCeU0icRMNlOKmAUJREReRLv4IiIfELBgs0gEyIi8homW/YJ7VTzYgJVKcLKvWTLqluIOpOi1ITybinRLICo8J4YLdntfJrYQrQiYB4t2SV8DikislvYvm7f0l+uyYyiDAgJjzVfrchWziecJURLSscUDMg7l86FNCm6UjjXdOdn0PBcl+qMnGxZ3LWmLkt13zAxOw0q3zZwRER+47coygEJMtm4cSMKCgqQmZmJ4uJi7N27t0/rbd26FZZlYdasWQNxWEREvqZcfKUC1xu4bdu2obKyElVVVdi/fz8mTpyIsrIyfPrpp9r1jh07hu9///v4+te/7vYhERGRD7newG3YsAGLFy9GeXk5xo0bh9raWmRlZaGurk5cJxqNYt68eVi7di2uvfZatw+JiIjweRelG69U4GoDF4lE0NTUhNLS0s93EAigtLQUjY2N4npPPPEERo8ejXvvvbdP++ns7ER7e3vci4iI9GwXX6nA1SCTU6dOIRqNIjc3N648NzcXBw8edFznnXfewb/8y7+gubm5z/uprq7G2rVr+3OoWsb56MRy+a8cJb0n5VEU8/NJEWeaKErhPSm3pBSV2KXLRSnuwzAXpfCXojYXpfDdSlGUUshzIpVY+otRygcpRVEmMlZJCfkrA0JuyYBh7koACFjO34r0ewcN80cC8jkt1QGxzkh1TFMvTeu4eK2gISGpmUzOnDmD+fPnY9OmTRg5cmSf11u5ciXa2tpir+PHjw/gURIRecPFcXBuvFKBq3dwI0eORDAYRGtra1x5a2sr8vLyei3/0Ucf4dixY5g5c2aszP7TX0RpaWk4dOgQrrvuul7rhUIhhEIhNw+diMjzOEygHzIyMlBUVIT6+vpYmW3bqK+vR0lJSa/lb775Zrz77rtobm6Ovf7qr/4Kd955J5qbm5Gfn+/m4RERkY+43kVZWVmJTZs24Wc/+xk++OADPPDAA+jo6EB5eTkAYMGCBVi5ciUAIDMzE+PHj497jRgxAl/5ylcwfvx4ZGRkuH14RES+lcxxcMkYH+16JpPZs2fj5MmTWLNmDcLhMAoLC7Fz585Y4ElLSwsCmlRAREQ0MJLVRXlxfHRtbS2Ki4tRU1ODsrIyHDp0CKNHjxbX6+/46AFJ1VVRUYGKigrH9xoaGrTrbt682f0DcqKdblsgRUx1C/kEo5ooSim6zDQXpRBNKC0PaKIohTyDUoSjbkZvKcpRLhciNQ2jK3XvyTkqB570awQDznu3pZyIkC8u6QEhL6I0c7fwFWpn9BaON2iYo1J3fpqe01KdkSOP5XPHEuqy8YzeiVxfPOyL46MBoLa2Ftu3b0ddXR1WrFjhuM4Xx0f/13/9F06fPm28X95KERH5hNvj4L48Hrmzs7PXPgdrfLQTNnBERD7h9jCB/Px85OTkxF7V1dW99qkbHx0Ohx2P8+L46E2bNvXr83I2ASIiSsjx48eRnZ0d+78bw7cSHR/thA0cEZFPKLjzzPnik8rs7Oy4Bs7JYI2PdsIuSiIin1BwqYvSIJVcMsdHe/8Ozu5KYB1hNl5hxmQoKfJK/vtBnIXYMHefNDNyNCr/tFI0oTgLt7RvTT5IacZtqVyKlpSX10SJGuaWFHNRJhBOLc2SLeV3lCJag5p8kOlCJKPpMcnlutm2nfedJkRwiuea5vyUzuk0MSpSKBfqmDYXpVjHzWb01krkmpTiKisrsXDhQkyePBlTpkxBTU1Nr/HRY8eORXV1dWx89BeNGDECAHqVX4r3GzgiIgIA2Krn5cZ2TCRrfDQbOCIin3BrNu5EtpGM8dF8BkdERJ7EOzgiIp/w22wCbOCIiHzCrdm4U2WaV3ZREhGRJ/n2Ds7ShfcaTlsv/TmjhHBnAFDKMKmyWC4NE5D/dnErqXKXZh8RYZ1OITxcGg7QKYWZ65ItG67j6l+jwj7EZMtCqL6tGQpga4ZIODOr5pYmhCAoDCHoDpol8I6mycmIpXPatG5IdUxXL8Wwf8Nrgvb6kkRuzcbtyxm9iYho6GIXJRERkQfwDo6IyCeUkhMvmW4nFbCBIyLyCRsWbIM8krrtpAJ2URIRkSfxDs6BZRgxJc1OrzRRhsbJlqVoNCEiTIp81L0nRkuKiZDlfZgmT5aiJd1MtixGUQ5Cd0tA+INXiqIMaqLUTJMtS6SkyrpEz13ClyUlQu4OmJ+f0jkt1QHTZMv6eulcLl0TxGvFEJWsXJTJwgaOiMgvXHoG50pCy0HALkoiIvIk3sEREfmE34JM2MAREfmE34YJsIuSiIg8yft3cFJYlDYXpWk+OmnXmlyUUg7JLuefRMzDJ0QrRoWcj4Ac/RgxzBMp5ZvUvSdFS0rLdwnRhNIxAUC38NelFEU5GH+NWmIUpXN5mrQC5DRJbk1hEtDlohQiOKVzJy3oXJd056d0Tkt1QKozcnSlLheleFBCuXSt0FxfpGvSIPBbqi7vN3BERATAf8ME2EVJRESexDs4IiKfUHBnCFuK3MCxgSMi8oueLkoXhgmkSAvHLkoiIvIk397BJTSjd7ewTrfwF5HmLyUpH54UKSbOwi3lotREism5KM0iHDs1uShNoyUjYs5J5+9Qiq4EdLkonZeX/qJNJLpSCn4MSDknhRWimnyQtjRTdcCdP6stIUclAASjznUjaDmXi+ea7vwUz2nnOpAuzegt5ZzU3cFIdVmq+yk3o7e/xsH5toEjIvIbvw0TYBclERF5Eu/giIh8gl2URETkSeyiJCIi8gDewTkR884535crKQpPEykm588T8u0Zztytm21byhsozcLdaZijsmcfZtGSnUK0pGl0JSBHSyZ3Rm/nfYszemtyUUaFaElb/Lva7O/YNE0Ep3QupBnmqMyw5XyM0jkt1QGxzkh1TJcjVjhHLOkkSbEZvZVLqbrYRUlEREOK3zKZsIuSiIg8iXdwREQ+4bfZBNjAERH5hN+GCbCLkoiIPIl3cEREPuG3cXDeb+CkcGRNMlQxUaqYhFkIPxfCnQFAiaHQzj+JLYTwSwloI0I5YJ5UWQrVv6D5fPI6ZsMEuhNIttwl/ExSt0q3C9OHXIoUem8JwwHSNX0r0ucQkzAb9icFNL+rJXyOdGGd9IBzXdKdnxmGdUCqM1Id09VLMdmyaVJlbTJ3eYjEQPPbMzh2URIRkSd5/w6OiIgA+G8cHBs4IiKfYBclERGRB/AOjojIJ/w2Ds63DZx2SnkxWlJItixFcQnJXgE5EawUFdklLN8lLC8lrAXkBLhiFKVQ3qlJtixFRYrlQnRll1CRdMmWxShKYfnuQYh5TgsISXyF5aWE0QAgfFVIF6NBhQhAIYLTsuTfNShFS1rO+0gPCMmWo3L9k85pqQ5IdUZMtqypl2Ii5u4u53LT6Mok89swAXZREhGRJw1IA7dx40YUFBQgMzMTxcXF2Lt3r7jspk2b8PWvfx1f/epX8dWvfhWlpaXa5YmIKDE2Pg806dcr2R+kj1xv4LZt24bKykpUVVVh//79mDhxIsrKyvDpp586Lt/Q0IA5c+bgrbfeQmNjI/Lz83HXXXfhxIkTbh8aEZGvKRdfqcD1Bm7Dhg1YvHgxysvLMW7cONTW1iIrKwt1dXWOy7/88st48MEHUVhYiJtvvhkvvfQSbNtGfX2924dGREQ+4moDF4lE0NTUhNLS0s93EAigtLQUjY2NfdrGuXPn0NXVhcsvv1xcprOzE+3t7XEvIiLSU4l0Rzq8fBlFeerUKUSjUeTm5saV5+bm4uDBg33axvLlyzFmzJi4RvLLqqursXbt2n4dq3aq+W4hV5wQ2qaEsDa7W/77wRbWiQpRXFJUZLeQn0+KfNS9F5G2ZZhXUveeFC15QQg6k3JOSpGSuvekyEQxt6O8C5H0jUSEjQWFwEdtFKXwnUSFPJFu/h0r5dRMF86d9KjzB9ednyEp76qYv1XKOWleL6W6LP4g0rVCd31JIqVcymSSIg3ckIqifPLJJ7F161a89tpryMzMFJdbuXIl2traYq/jx48P4lESEVEqcPUObuTIkQgGg2htbY0rb21tRV5ennbdZ555Bk8++SR+/etf49Zbb9UuGwqFEAqF+n28RER+wnFw/ZCRkYGioqK4AJGLASMlJSXiek899RTWrVuHnTt3YvLkyW4eEhER/UnPMzTlwivZn6RvXM9kUllZiYULF2Ly5MmYMmUKampq0NHRgfLycgDAggULMHbsWFRXVwMAfvjDH2LNmjXYsmULCgoKEA6HAQDDhw/H8OHD3T48IiLyCdcbuNmzZ+PkyZNYs2YNwuEwCgsLsXPnzljgSUtLCwKBz28cX3jhBUQiEXz3u9+N205VVRUef/xxtw+PiMi3OF2OCyoqKlBRUeH4XkNDQ9z/jx07NhCH8Dlh9tyEclFKOQ6FKC5pRuGewxLy5wl5Ik1n7pbyTQJApxClJs3QLZVL0ZWAebTkBTFHpfPyuvyRXUKIl9StootYdIsULSmkqBRnMgeALuFrD0k7ES9HzhuSjgmQf7+gEF2ZFnD+oTI056d0Tos5J4VtSXVMVy/FXJTS+ZZILspkz+jt0nZSwZCKoiQiInKLb2cTICLyG/Wnf25sJxWwgSMi8gl2URIREXkA7+CIiHzCbwO9/dvA6aKclPPPp4TgJykqy9bN6C3lnBTKu8RZuM1m5wbMoyUvRJ3D6s5rcvqdE9YRoyiF8ojQF6LPRem8TlSKrhS2oxJIuCfNki19U0Fh+XRN30q6GA3qvC1biK6Uv0J55wHh2UvQcj53pCjKUMB8xnmpDkh1RsxRqZvRW6gDUt23hGuF9vqSREq59AwuRZJRsouSiIg8yb93cEREPsMuSiIi8iR2URIREXkA7+CIiHxCwZ3uxdS4f2MDR0TkG7ZSsF1onuwU6aL0fgMnxffqwnilaei7hXBrIexYClPu2YUQCm2cVFkK7dckWxbXkRIkmyVO1r13XvhqLwgZjzulYQKaCtYlJMCNChVbqqyJPKuwIAwTEIYDBIXl05X89CBdSMQs/ExQwvABie65RdByflfK85wuDAe4EJDPz1DU+SSR6oBUZ8TkzJp6KQ4hEOq+eK3QDkNKXrLlZNq4cSOefvpphMNhTJw4ET/+8Y8xZcoUx2U3bdqEn//853jvvfcAAEVFRVi/fr24vITP4IiIfEK5+M/Etm3bUFlZiaqqKuzfvx8TJ05EWVkZPv30U8flGxoaMGfOHLz11ltobGxEfn4+7rrrLpw4ccJov2zgiIh8wnbxZWLDhg1YvHgxysvLMW7cONTW1iIrKwt1dXWOy7/88st48MEHUVhYiJtvvhkvvfQSbNtGfX290X7ZwBERUULa29vjXp2dnb2WiUQiaGpqQmlpaawsEAigtLQUjY2NfdrPuXPn0NXVhcsvv9zo+NjAERH5hA3l2gsA8vPzkZOTE3tVV1f32uepU6cQjUaRm5sbV56bm4twONyn416+fDnGjBkT10j2hfeDTIiICID7UZTHjx9HdnZ2rDwUCvV721/25JNPYuvWrWhoaEBmZqbRup5v4CxhenhLiLQDIE5Dr4QEwkpIOhzVJHXtFt7r6kp3LJcS0HYKkWJS4mTdexeE6DwpcbJUDgAdhtGSF6LO33mn8FtElByl1iU8IeiG8zpShXc1ilKIZEyD82+RLoVEAsgQEhtHhchLW+iosYVjhVgOCMGgCAaEaFDhXMsQkjADQKYYMWxWZ6Q6pquXUl2W6r50HdFdX6RrUirKzs6Oa+CcjBw5EsFgEK2trXHlra2tyMvL0677zDPP4Mknn8Svf/1r3HrrrcbHxy5KIiKfSEYUZUZGBoqKiuICRC4GjJSUlIjrPfXUU1i3bh127tyJyZMnJ/R5PX8HR0REPb74/Ky/2zFRWVmJhQsXYvLkyZgyZQpqamrQ0dGB8vJyAMCCBQswduzY2DO8H/7wh1izZg22bNmCgoKC2LO64cOHY/jw4X3eLxs4IiIaULNnz8bJkyexZs0ahMNhFBYWYufOnbHAk5aWFgQCn3covvDCC4hEIvjud78bt52qqio8/vjjfd4vGzgiIp9I1h0cAFRUVKCiosLxvYaGhrj/Hzt2LIGj6o0NHBGRTySShUTaTirwbwOnzUXp/J4ScthJ+etsIRoMkPPkSZFiEWHfnbaUV1IXRSlERQoRZOeF5aW8kj3bcq4A56PO3+154feICJGPnegS9x2xnN+LwvmAu62Bj2pLU86/a1CoghnKOTIQAELCe13Cbx5VwnkrlAc0UZRiTk0hR2W65XweZAg5KgH5nJbqgFRnpDqmq5dSXZbqvnSt0F5faND4t4EjIvIZ5VIXJe/giIhoSLEtG5bV/xnhbFdmlRt4HAdHRESexDs4IiKfsKFgJSmKMhnYwBER+cTFVMlubCcVeL+Bk6KZdLkopSjKrgznTUkzDWty3pnO3N0pRIqdF3NRyr3P54T3pGjJDiFg8awQKQkAHUK05DkhD9859J5mAwA6Lan8grjviBVxLI8KkZdRJUdkuiUYcI58DMK5PEM5n2sAEFLOCWdDyjnRbdR2LpdrgBzhKOWpDArRlekB53MtIyCfO9I5nZVmVmekOqarl1JdVl1SFKXzuaa9vjDCctB4v4EjIiIAPX/UuNNFmRrYwBER+QSjKImIiDyAd3BERD5hw4blwt1XqtzBsYEjIvIJNnA+YXVrIueE6EDTXJRRTc47aRbizm7n8guGM3dLkZIA0CHknJRm4ZaiJc90yzkcO4TIxA4h+vGcddaxvNM671geUc7lANClnPfRbTtHZEohz0oza7jEEmbbtoSnAWkB5wjHdDhHSgJAZ2CYY3lIOZd3K+f5s6K28z5sTR5M6ZIhzegt5ajURVFmSjPOSxHGQp2R6piuXprnonT+HNrrCw0a3zZwRER+w3FwRETkSYyiJCIi8gDewRER+YSC7crdF7soiYhoSFGIQrnQcaeEiYiHGnZREhGRJ3n+Ds4Skvtqk6FGhWECQghxd8Q5HLlbCFMG5NDmTiGE+ZyQbFkaDiANBQCAc8JX0tHl/LmlxMnSUAAAaBfC/juE8vNodyzvtJ2X77LPifuWhgNEbefEuPYgJFsOWEKy5YBzUmVp+AAApAeyHMsjlvNwgG7L+QePwnl5CMMKACAgJOROE5IRpwlJmDOEJMwAkBl0fm+YUAcuMxw+oKuXUl2W6r50rdBdX8Rr0iDo6Z70T5CJ5xs4IiLq0TOPmxsNXGrMB8cuSiIi8iTewRER+URPkIlzt7HpdlIBGzgiIp/w2zM4dlESEZEnef8OTpgeXpcMVUmBl0IklS1EkHUKEVkA0ClEcl0QIsXOC5FiHUIS2I5uuRtCTp7s/F2dEaIPT1vOkY8A0B74f47l52zn8ogQLdnZfcaxPGrLyZZtJSRVFqMlB6O7RUjCLERXBiw5ijIoJFvuSnOOLI0EnMu7Al91LLd1EcZ2tmNxoNs5GjQYcP7c6UJyZgDIDDqvkxWUEo4LUZRSQnNNvZTqspyE2Xk72mTLwjVpMDAXJREReZKNKODCMzg7RZ7BDUgX5caNG1FQUIDMzEwUFxdj79692uV/8Ytf4Oabb0ZmZiYmTJiAHTt2DMRhERGRj7jewG3btg2VlZWoqqrC/v37MXHiRJSVleHTTz91XH737t2YM2cO7r33Xhw4cACzZs3CrFmz8N5777l9aEREvnaxi9KNVypwvYHbsGEDFi9ejPLycowbNw61tbXIyspCXV2d4/LPPfccvvWtb+Hhhx/GLbfcgnXr1uG2227DT37yE7cPjYjI12wVde2VClx9BheJRNDU1ISVK1fGygKBAEpLS9HY2Oi4TmNjIyorK+PKysrK8Prrr4v76ezsRGfn54EEbW1tAID29t6BB+fPOD/sjXTIf4Go885BGJ2dzuucjTj/2Gd1M14LD6HPRZ0DOs4Lsxx32s7HGrHlfnYhIxe6hJO2WzkfUxTyg3RbePouVQxp9mylpNm25UwK8num5W4S0r8Jx6o0xyR/J87fofSdS79RVJO6TDoX5HPK+byNaG4ApHP6fNSszkh1LFNTL0NCXVZC3Q8J1wpLc325IFyTutrjg4EuXs905zrpudrAnTp1CtFoFLm5uXHlubm5OHjwoOM64XDYcflwOCzup7q6GmvXru1VXnDVkgSO2kSrYTnRRdIFz/liZys516YUhNdl+Ee1czyry+SP4XF/SOC91xxL//jHPyInJ6ffRwQwijIlrFy5Mu6u7/Tp07j66qvR0tLi2ongB+3t7cjPz8fx48eRne0c/k3x+J0lht+buba2Nlx11VW4/PLLXdtmTwPX/+5FXzZwI0eORDAYRGtr/B1Na2sr8vLyHNfJy8szWh4AQqEQQqHe44RycnJYeRKQnZ3N780Qv7PE8HszF9DMvEB6rn5zGRkZKCoqQn19fazMtm3U19ejpKTEcZ2SkpK45QFg165d4vJERJQYpWzYLryk58BDjetdlJWVlVi4cCEmT56MKVOmoKamBh0dHSgvLwcALFiwAGPHjkV1dTUAYMmSJZg2bRqeffZZzJgxA1u3bsW+ffvw4osvun1oRES+1tO16EayZZ82cLNnz8bJkyexZs0ahMNhFBYWYufOnbFAkpaWlrhb7qlTp2LLli1YtWoVHn30Udxwww14/fXXMX78+D7vMxQKoaqqyrHbkmT83szxO0sMvzdz/M76z1KMQSUi8rT29nbk5OQgJ3McLMt56IYJpaJou/B7tLW1DelnqikZRUlEROZs2LB81EXJ8BwiIvIk3sEREflET/SjC3dwfo2iJCKiocmNQd5ubmegsYuSiIg8KWUbuB/84AeYOnUqsrKyMGLEiD6to5TCmjVrcMUVV2DYsGEoLS3F//7v/w7sgQ4hn332GebNm4fs7GyMGDEC9957L86edZ5J+6Lp06fDsqy41/333z9IR5wcnM/QnMl3tnnz5l7nVGZm5iAe7dDw9ttvY+bMmRgzZgwsy9ImmL+ooaEBt912G0KhEK6//nps3rzZaJ9KKag/DdTu3ys1gu9TtoGLRCL467/+azzwwAN9Xuepp57Cj370I9TW1mLPnj247LLLUFZWhgsXLgzgkQ4d8+bNw/vvv49du3bhV7/6Fd5++23cd999l1xv8eLF+OSTT2Kvp556ahCONjk4n6E50+8M6EnZ9cVz6uOPPx7EIx4aOjo6MHHiRGzcuLFPyx89ehQzZszAnXfeiebmZixduhSLFi3Cm2++2ed9+m0+OKgU99Of/lTl5ORccjnbtlVeXp56+umnY2WnT59WoVBI/du//dsAHuHQ8Pvf/14BUP/zP/8TK/uP//gPZVmWOnHihLjetGnT1JIlSwbhCIeGKVOmqIceeij2/2g0qsaMGaOqq6sdl/+bv/kbNWPGjLiy4uJi9fd///cDepxDiel31tc66ycA1GuvvaZd5pFHHlFf+9rX4spmz56tysrKLrn9trY2BUANyyhQWaFr+/0allGgAKi2trb+fOwBl7J3cKaOHj2KcDiM0tLSWFlOTg6Ki4vFueq8pLGxESNGjMDkyZNjZaWlpQgEAtizZ4923ZdffhkjR47E+PHjsXLlSpw75805UC7OZ/jFc6Qv8xl+cXmgZz5DP5xTQGLfGQCcPXsWV199NfLz8/Gd73wH77///mAcbkpz41xTKuraKxX4Jory4vxypnPPeUU4HMbo0aPjytLS0nD55ZdrP//cuXNx9dVXY8yYMfjd736H5cuX49ChQ3j11VcH+pAH3WDNZ+gliXxnN910E+rq6nDrrbeira0NzzzzDKZOnYr3338fV1555WAcdkqSzrX29nacP38ew4YNu+Q23ArvT5VhAkPqDm7FihW9Hj5/+SVVGr8a6O/svvvuQ1lZGSZMmIB58+bh5z//OV577TV89NFHLn4K8pOSkhIsWLAAhYWFmDZtGl599VWMGjUK//zP/5zsQyOPGVJ3cP/4j/+Ie+65R7vMtddem9C2L84v19raiiuuuCJW3traisLCwoS2ORT09TvLy8vr9dC/u7sbn332mXbuvS8rLi4GABw+fBjXXXed8fEOZYM1n6GXJPKdfVl6ejomTZqEw4cPD8QheoZ0rmVnZ/fp7g1wL8VWqgSZDKkGbtSoURg1atSAbPuaa65BXl4e6uvrYw1ae3s79uzZYxSJOdT09TsrKSnB6dOn0dTUhKKiIgDAb37zG9i2HWu0+qK5uRkA4v5I8Iovzmc4a9YsAJ/PZ1hRUeG4zsX5DJcuXRor89N8hol8Z18WjUbx7rvv4u677x7AI019JSUlvYagmJ5rfuuiTNkoyo8//lgdOHBArV27Vg0fPlwdOHBAHThwQJ05cya2zE033aReffXV2P+ffPJJNWLECPXGG2+o3/3ud+o73/mOuuaaa9T58+eT8REG3be+9S01adIktWfPHvXOO++oG264Qc2ZMyf2/h/+8Ad10003qT179iillDp8+LB64okn1L59+9TRo0fVG2+8oa699lr1jW98I1kfYcBt3bpVhUIhtXnzZvX73/9e3XfffWrEiBEqHA4rpZSaP3++WrFiRWz5//7v/1ZpaWnqmWeeUR988IGqqqpS6enp6t13303WRxh0pt/Z2rVr1Ztvvqk++ugj1dTUpL73ve+pzMxM9f777yfrIyTFmTNnYtctAGrDhg3qwIED6uOPP1ZKKbVixQo1f/782PJHjhxRWVlZ6uGHH1YffPCB2rhxowoGg2rnzp2X3NfFKMr0YK7KSLui36/0YG5KRFGmbAO3cOFCBaDX66233ootA0D99Kc/jf3ftm21evVqlZubq0KhkPrmN7+pDh06NPgHnyR//OMf1Zw5c9Tw4cNVdna2Ki8vj/uD4OjRo3HfYUtLi/rGN76hLr/8chUKhdT111+vHn744SF/UvfXj3/8Y3XVVVepjIwMNWXKFPXb3/429t60adPUwoUL45Z/5ZVX1I033qgyMjLU1772NbV9+/ZBPuLkM/nOli5dGls2NzdX3X333Wr//v1JOOrkeuuttxyvYRe/q4ULF6pp06b1WqewsFBlZGSoa6+9Nu76pnOxgUsLjlLpabn9fqUFR6VEA8f54IiIPO7ifHDBwOWwrP7HFiplI2p/NuTngxtSUZRERERuGVJBJkRENJAU4EoEZGp0/LGBIyLyCffmg0uNBo5dlERE5Em8gyMi8omeAdou3MGxi5KIiIYWdxq4VHkGxy5KIiLyJN7BERH5hUtBJkiRIBM2cEREPuG3Z3DsoiQiIk9iA0d0CSdPnkReXh7Wr18fK9u9ezcyMjJQX1+fxCMjMmW7+Br62MARXcKoUaNQV1eHxx9/HPv27cOZM2cwf/58VFRU4Jvf/GayD4/IgOp5ftbfVwJdlBs3bkRBQQEyMzNRXFyMvXv3apf/xS9+gZtvvhmZmZmYMGFCr6mC+oINHFEf3H333Vi8eDHmzZuH+++/H5dddhmqq6uTfVhEKWHbtm2orKxEVVUV9u/fj4kTJ6KsrKzXJMwX7d69G3PmzMG9996LAwcOYNasWZg1axbee+89o/1yNgGiPjp//jzGjx+P48ePo6mpCRMmTEj2IRH1ycXZBIAg3BsHF+3zbALFxcX48z//c/zkJz8B0DMpbn5+Pv7hH/4BK1as6LX87Nmz0dHRgV/96lexsr/4i79AYWEhamtr+3yUvIMj6qOPPvoI//d//wfbtnHs2LFkHw5RghynoTN89Whvb497dXZ29tpbJBJBU1MTSktLY2WBQAClpaVobGx0PMLGxsa45QGgrKxMXF7CBo6oDyKRCP72b/8Ws2fPxrp167Bo0SKxe4VoqMnIyEBeXh6AqGuv4cOHIz8/Hzk5ObGXU7f9qVOnEI1GkZubG1eem5uLcDjseLzhcNhoeQnHwRH1wWOPPYa2tjb86Ec/wvDhw7Fjxw783d/9XVwXCtFQlZmZiaNHjyISibi2TaUULCu+uzMUCrm2fTewgSO6hIaGBtTU1OCtt96KPW/413/9V0ycOBEvvPACHnjggSQfIdGlZWZmIjMzc9D3O3LkSASDQbS2tsaVt7a2/umusre8vDyj5SXsoiS6hOnTp6Orqwt33HFHrKygoABtbW1s3IguISMjA0VFRXFjRm3bRn19PUpKShzXKSkp6TXGdNeuXeLyEt7BERHRgKqsrMTChQsxefJkTJkyBTU1Nejo6EB5eTkAYMGCBRg7dmzsGd6SJUswbdo0PPvss5gxYwa2bt2Kffv24cUXXzTaLxs4IiIaULNnz8bJkyexZs0ahMNhFBYWYufOnbFAkpaWFgQCn3coTp06FVu2bMGqVavw6KOP4oYbbsDrr7+O8ePHG+2X4+CIiMiT+AyOiIg8iQ0cERF5Ehs4IiLyJDZwRETkSWzgiIjIk9jAERGRJ7GBIyIiT2IDR0REnsQGjoiIPIkNHBEReRIbOCIi8qT/D4VCuJc0yamLAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(5, 5))\n", "plt.imshow(\n", " sol.ys.vals,\n", " origin=\"lower\",\n", " extent=(x0, x_final, t0, t_final),\n", " aspect=(x_final - x0) / (t_final - t0),\n", " cmap=\"inferno\",\n", ")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"t\", rotation=0)\n", "plt.clim(0, 1)\n", "plt.colorbar()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "26ba8fec-3ca9-4612-b6f9-83962333d96d", "metadata": {}, "source": [ "That worked!\n", "\n", "However, for more complicated PDEs then we may wish to define a custom solver. So as an example, here's how to solve the same PDE using the famous [Crank–Nicolson](https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method) scheme.\n", "\n", "(See the page on [abstract solvers](https://docs.kidger.site/diffrax/api/solvers/abstract_solvers/) for more details about how to define a custom solver.)" ] }, { "cell_type": "code", "execution_count": 6, "id": "059fed69-c042-4fec-bf36-60e365c98de8", "metadata": { "tags": [] }, "outputs": [], "source": [ "class CrankNicolson(diffrax.AbstractSolver):\n", " rtol: float\n", " atol: float\n", "\n", " term_structure = diffrax.ODETerm\n", " interpolation_cls = diffrax.LocalLinearInterpolation\n", "\n", " def order(self, terms):\n", " return 2\n", "\n", " def init(self, terms, t0, t1, y0, args):\n", " return None\n", "\n", " def step(self, terms, t0, t1, y0, args, solver_state, made_jump):\n", " del solver_state, made_jump\n", " δt = t1 - t0\n", " f0 = terms.vf(t0, y0, args)\n", "\n", " def keep_iterating(val):\n", " _, not_converged = val\n", " return not_converged\n", "\n", " def fixed_point_iteration(val):\n", " y1, _ = val\n", " new_y1 = y0 + 0.5 * δt * (f0 + terms.vf(t1, y1, args))\n", " diff = jnp.abs((new_y1 - y1).vals)\n", " max_y1 = jnp.maximum(jnp.abs(y1.vals), jnp.abs(new_y1.vals))\n", " scale = self.atol + self.rtol * max_y1\n", " not_converged = jnp.any(diff > scale)\n", " return new_y1, not_converged\n", "\n", " euler_y1 = y0 + δt * f0\n", " y1, _ = lax.while_loop(keep_iterating, fixed_point_iteration, (euler_y1, False))\n", "\n", " y_error = y1 - euler_y1\n", " dense_info = dict(y0=y0, y1=y1)\n", "\n", " solver_state = None\n", " result = diffrax.RESULTS.successful\n", " return y1, y_error, dense_info, solver_state, result\n", "\n", " def func(self, terms, t0, y0, args):\n", " return terms.vf(t0, y0, args)" ] }, { "cell_type": "code", "execution_count": 7, "id": "da4511b8-f112-4839-94f5-dfc7728da8ea", "metadata": { "tags": [] }, "outputs": [], "source": [ "solver = CrankNicolson(rtol=rtol, atol=atol)\n", "sol = diffrax.diffeqsolve(\n", " term,\n", " solver,\n", " t0,\n", " t_final,\n", " δt,\n", " y0,\n", " saveat=saveat,\n", " stepsize_controller=stepsize_controller,\n", " max_steps=None,\n", ")" ] }, { "cell_type": "code", "execution_count": 8, "id": "6667e3c7-5b45-4740-9caf-3e0aa4b1d7a9", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAGiCAYAAACVh9NOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABENElEQVR4nO3df3RU5Z0/8PedSTIh0oS6QCIYjb+VggTDkg3agqdZ08phyx/dpcACZgXXH9kDZKuAAhE9JdYfbGyLZsVNaffIgu3xx57C4qGpWdclhSWQU7WCi4ChrBOhfkkgQCaZ+3z/SBkZcz8PeSY3mcy97xdn/uCZ+ysz97nP3Od+ns9jKaUUiIiIPCaQ7AMgIiIaCGzgiIjIk9jAERGRJ7GBIyIiT2IDR0REnsQGjoiIPIkNHBEReRIbOCIi8iQ2cERE5Els4IiIyJPYwBER0YB65513MHPmTIwZMwaWZeGNN9645DoNDQ247bbbEAqFcP3112PTpk3G+2UDR0REA6qjowMTJ07Ehg0b+rT8kSNHMGPGDNx5551obm7G0qVLsWjRIrz11ltG+7WYbJmIiAaLZVl4/fXXMWvWLHGZ5cuXY9u2bXj//fdjZd/73vdw6tQp7Nixo8/7SuvPgRIRUWo4f/48IpGIa9tTSsGyrLiyUCiEUCjU7203NjaitLQ0rqysrAxLly412g4bOCIijzt//jyuuSYP4XCba9scPnw4zpw5E1dWVVWFxx9/vN/bDofDyM3NjSvLzc1Fe3s7zp07h2HDhvVpO2zgiIg8LhKJIBxuw+FP/gnZ2X1rHHTa28/h2quX4dixY8jOzo6Vu3H35iY2cEREPpGdPcyVBu6L7WXHNXBuycvLQ2tra1xZa2srsrOz+3z3BrCBIyLyDaW6oVS3K9sZSCUlJdi+fXtc2c6dO1FSUmK0HQ4TICLyCaWirr1MnDlzBs3NzWhubgbQMwygubkZLS0tAICVK1diwYIFseXvv/9+HD58GI888ggOHDiAF154Aa+++iqWLVtmtF82cERENKD27t2LSZMmYdKkSQCAyspKTJo0CWvWrAEAfPrpp7HGDgCuueYabNu2DTt37sTEiRPx3HPP4eWXX0ZZWZnRfjkOjojI49rb25GTk4PwyedcCzLJG/mPaGtrG5BncG7hMzgiIp9IlWdwbmEXJREReRLv4IiIfKInQMSNOzizIJNkYQNHROQTyu6Gsl1o4FzYxmBgFyUREXkS7+CIiPxCdfe83NhOCmADR0TkE4yiJCIi8gDewRER+YXdDdhd7mwnBbCBIyLyiZ4uyqAr20kF7KIkIiJP4h0cEZFf2N2A3f87OHZREhHR0OKzBo5dlERE5Em8gyMi8o2oS4O0mYuSiIiGEMvuhmX3v+POYhclERFR8vAOjojIL+xuwIU7uFQJMmEDR0TkFz5r4NhFSUREnsQ7OCIin7BUNyzlQpBJiqTqYgNHROQXtg3YLoT423b/tzEI2EVJRESexDs4IiKf6BkHZ7mynVTABo6IyC/sqEtRlKmRyYRdlERE5EmuN3DvvPMOZs6ciTFjxsCyLLzxxhuXXKehoQG33XYbQqEQrr/+emzatMntwyIiIrvbvVcKcL2B6+jowMSJE7Fhw4Y+LX/kyBHMmDEDd955J5qbm7F06VIsWrQIb731ltuHRkTka5Ydde2VClx/Bvftb38b3/72t/u8fG1tLa655ho899xzAIBbbrkF7777Lv7pn/4JZWVlbh8eERH5RNKDTBobG1FaWhpXVlZWhqVLl4rrdHZ2orOzM/Z/27bx+eef48/+7M9gWf2PECIiSjalFE6fPo0xY8YgEHCps025FGSifHoHZyocDiM3NzeuLDc3F+3t7Th37hyGDRvWa53q6mqsXbt2sA6RiChpjh07hiuvvNKVbVm27Ur3opUiA72T3sAlYuXKlaisrIz9v62tDVdddRWOtjyP7Oz4BvF85wnHbUTP/kHcfuCs8zqBs390LA+eOeW8/Lkz4j6ssx3Ob5w/51x+9rxzecT5Ya8SFgcAFXH+BacizqeD6nIut7uD4j5sYR3VLew76rwtJYzZUba8b6WkdQzv7hMZLxRQRotbwvKWJW/HCjhfoMRtBYXl05wvUoF0OYAgkCZsS1jHypDK5QuklSm8kSFcrrKEFTJ7/zgGAJV1mbhve9hwx/Lo8BHOy2f9mVA+StxHMMu5scoMxa/T3n4OBVctwVe+8hVxW6SX9AYuLy8Pra2tcWWtra3Izs52vHsDgFAohFAo1Ks8O3sYsrOz4soyOp1P/mha7/UvCATTncsDzh9XEM4X20BAvghbltBNEJAuqkJ5ULiYa67NStiHWC7sw+6SuzpsoUtFBaUGTiiXGjhhecDNBi6BrpyA2S/bhBq4oLCO2MAJ5ULtD6TLf3cg3fnvs9KdP1srQygPyd+F2MBJ62QJ5cOEcypL88NMeC86XPiRd5nztcK+TL6+BIUGOTOU5Vju6mMXO5rYDzen7aSApDdwJSUl2L59e1zZzp07UVJSkqQjIiLypp4ISDcymaRGA+f6MIEzZ86gubkZzc3NAHqGATQ3N6OlpQVAT/figgULYsvff//9OHz4MB555BEcOHAAL7zwAl599VUsW7bM7UMjIiIfcf0Obu/evbjzzjtj/7/wrGzhwoXYtGkTPv3001hjBwDXXHMNtm3bhmXLluH555/HlVdeiZdffplDBIiI3MYuyv6ZPn06lJKfHzhlKZk+fTr279/v9qEkn+4kkKKQTKOTpMUTCXIyfXalqyjSe8JcVPKzNiH4RPOQUT5eYd+6B5amhOOVnqkp6YvS9a0I+wCk4BNhY7ZQTzXfq/TZWqafoe78ND13TetSilycBwK7KImIiDwg6UEmREQ0SNhFSUREXmTZypVB2pbUvT3EsIuSiIg8iXdwRER+YUcTC0Bz2k4K8HwDp1QC8xYJX54YOTQYX7ZxdKV5hKMYfZhAZhA5xZbhPhLJSmIYLWmc4SQRQl+JJR2r7iokbcv07xMyZOg+D9N9iBF72s9ciu4chPyHg1D3E7omuUW51MClSLJldlESEZEnef4OjoiIeljKNh+zKGwnFbCBIyLyC589g2MXJREReRLv4IiI/MK2XRrozS5KIiIaStjAkWt0J4FxglgpdNrskABASRN5SiHrYti9pofbNKmyuG+zxMkAYGsmQ3Xch5vJliXC96SEiWx1f4E0hEBJk+gKiZ6lJMy671X6nsThDtLyuhNXTCBuOHzArYTmlLLYwBER+YRl27BcaN/dSPc1GNjAERH5hW27FEWZGg0coyiJiMiTeAdHROQXPruDYwNHROQXbOB8wnYzCXMCX7ZL0ZJSxhxtZKBpIuQEEua6lVRZipbURUrKSZUHPrrSEiIWpdmzrIDzF6j7+6R3VECIlnQrQTJgfO5I29Kfn8Lf4VZ0pYZYl00zdyRyfSHX+beBIyLyGxWVfxAYbYd3cERENIT4bZgAoyiJiMiTeAdHROQXDDIhIiJPYgNHrkVADcZJIEaWadaRItgMo+3ESEntOkLOwu6gc7lp1KVuH4Z/XyKk3JKWEOEo5tQUoisBzd8ufYZCZKeUozKR71WcRFOMjhV34c4FWMfNesloySGNDRwRkV/Yyp0G3o1IzEHABo6IyC9s5VIXZWo0cIyiJCIiT+IdHBGRX7g24Snv4IiIaCixbfdehjZs2ICCggJkZmaiuLgYe/bs0S5fU1ODm266CcOGDUN+fj6WLVuG8+fPG+2Td3BuMM1TB8ipbqQTRxnO6O1mhKM0w7N2Rm/nfUg5FqXIQHF5TV5J02jJZM7oLf7E1M1YLm0qaJbXMihFaupm9BZnajebkV2f71LcubC8YaJWnUTqMl3S1q1bUVlZidraWhQXF6OmpgZlZWU4ePAgRo8e3Wv5zZs3Y8WKFairq8PUqVPx0Ucf4Z577oFlWVi/fn2f98s7OCIiv7CVey8D69evx+LFi1FeXo5x48ahtrYWWVlZqKurc1x+165duP322zF37lwUFBTgrrvuwpw5cy551/dlbOCIiPxC2e69ALS3t8e9Ojs7e+0yEomgqakJpaWlsbJAIIDS0lI0NjY6HubUqVPR1NQUa9AOHz6M7du34+677zb6c9nAERFRQvLz85GTkxN7VVdX91rm5MmTiEajyM3NjSvPzc1FOBx23O7cuXPxxBNP4I477kB6ejquu+46TJ8+HY8++qjR8fEZHBGRXyiXxsH96XnosWPHkJ2dHSsOhUIubBxoaGjAunXr8MILL6C4uBiHDh3CkiVL8OSTT2L16tV93g4bOCIiv3B5oHd2dnZcA+dk5MiRCAaDaG1tjStvbW1FXl6e4zqrV6/G/PnzsWjRIgDAhAkT0NHRgfvuuw+PPfYYAoG+dT6ygTNgGUZYaZd3aUZveYZlzQkgvWc4c7fdrYm2E/dhFm0nbSeRGb1Nx/+I0aMaljDZlpSjUgmRgdLM4IBmRm9hHcsyi3y0u+V9B4Oms76bR8HKs2pLdcMsulJXL8WZ1xldmbCMjAwUFRWhvr4es2bNAgDYto36+npUVFQ4rnP27NlejVgw2JNrVaozTtjAERH5RZJSdVVWVmLhwoWYPHkypkyZgpqaGnR0dKC8vBwAsGDBAowdOzb2DG/mzJlYv349Jk2aFOuiXL16NWbOnBlr6PqCDRwRkU9cFADZ7+2YmD17Nk6cOIE1a9YgHA6jsLAQO3bsiAWetLS0xN2xrVq1CpZlYdWqVTh+/DhGjRqFmTNn4gc/+IHRftnAERHRgKuoqBC7JBsaGuL+n5aWhqqqKlRVVfVrn2zgiIj8wmezCbCBIyLyCxsuNXAubGMQcKA3ERF5Eu/g3CCGNSeS7FXalnNxIiHxcki+UB4VopY0SYpNkyq7Vd6zc2kds2NKhFLOn5UU9i+FxGvzWLv0WRknYYZ8Lqigcxh9IsmWpeO1DIcDiHTLuzHb9VDmszs4NnBERH6hIA/2M91OCmAXJREReRLv4IiIfELZlthtbLYdFw5mELCBIyLyC589g2MXJREReZJ/7+DsbvEtMbGqmxFWhgllxS4BqVwXGShG27mTCLnnuMySJNtCdJ6YVFnTzWLbQqSfcbJl864cKVpSSrZsBZzLA5CT+9ou/S4NSseq+V6VFPUphX1Kn6E2Clbat3O5a9GVOqaJmzXXl6RSlnHSceft9H8Tg8G/DRwRkc/47RncgHRRbtiwAQUFBcjMzERxcXFs2nFJTU0NbrrpJgwbNgz5+flYtmwZzp8/PxCHRkREPuH6HdzWrVtRWVmJ2tpaFBcXo6amBmVlZTh48CBGjx7da/nNmzdjxYoVqKurw9SpU/HRRx/hnnvugWVZWL9+vduHR0TkX7ZLXZR+vYNbv349Fi9ejPLycowbNw61tbXIyspCXV2d4/K7du3C7bffjrlz56KgoAB33XUX5syZc8m7PiIiMqQs914pwNUGLhKJoKmpCaWlpV/sIBBAaWkpGhsbHdeZOnUqmpqaYg3a4cOHsX37dtx9993ifjo7O9He3h73IiIiupirXZQnT55ENBqNTWJ3QW5uLg4cOOC4zty5c3Hy5EnccccdUEqhu7sb999/Px599FFxP9XV1Vi7dq2bh94vUj5BV3PeiVGUZnklAUB1C7klTXNUan7FiXkRhQg98XiFcilSUrctMSejm79GhW0FhIhF6Xu1If99UoSlsoS/25JycEoRkeZ5TMVzRzjXVLouilnah7iKsLx5vRTrskcwyGSQNTQ0YN26dXjhhRewb98+vPbaa9i2bRuefPJJcZ2VK1eira0t9jp27NggHjERUYqyA+69UoCrd3AjR45EMBhEa2trXHlrayvy8vIc11m9ejXmz5+PRYsWAQAmTJiAjo4O3HfffXjsscfipjG/IBQKIRQKuXnoRETkMa42wxkZGSgqKkJ9fX2szLZt1NfXo6SkxHGds2fP9mrEgsGebg2lUmQ0IRFRKrgQRenGKwW4PkygsrISCxcuxOTJkzFlyhTU1NSgo6MD5eXlAIAFCxZg7NixqK6uBgDMnDkT69evx6RJk1BcXIxDhw5h9erVmDlzZqyhIyKi/lPKcmX+w1S593C9gZs9ezZOnDiBNWvWIBwOo7CwEDt27IgFnrS0tMTdsa1atQqWZWHVqlU4fvw4Ro0ahZkzZ+IHP/iB24dGREQ+MiCpuioqKlBRUeH4XkNDQ/wBpKWhqqoKVVVVA3EoUInkhFPO64h55yTaKEqzqDo5P5+Uq1GXT9CdmbvFaEwkMKO3lKMygbySUTFS06xH3s1clNKZYwmzZwc1IYNShKX016mA8zFJ35FlaXpOpJm7DWf61p+fwnGJeR+FDbmYo1Lct3Ct0EnomuQWO+DSQO/UuIVjLkoiIp9QtnnSceftpEYDlxqxnkRERIZ4B0dE5BeuTZfj0yhKIiIamtyLokyNBo5dlERE5Em8g3ODaXSljjTAxDQHoC4XpeHM3dJs2/qchWYzepsekxQpqdu3GMHp4q9RaVtSdKUUTag7o6QISynnpB113rc4+7gm0aB0LlhBs5m+deennItSmh3cxYAHN+vyUORWmq0UyUXJBo6IyCfcS7bMLkoiIqKk4R0cEZFP+C3IhA0cEZFf+OwZHLsoiYjIk3gHR0TkE34LMvFvA5dAklTjEGJtsmXhPSmpstTnLSYv1iTMNRxaIA8fkDsAxPcMtyUNB5DC1QHNcADp7xuE5wlSSL4lJEK2XDymoOV8UonJlh0mGY69FzQclmKa2BsQ66Y4BMM0qXICyZaN634i15dB4LdncOyiJCIiT/LvHRwRkd/4LMiEDRwRkU/47RkcuyiJiMiTeAdHROQTfgsyYQPnQJyeXiJFZWmS1opTvkuriLswTEyrWUeMbJOSF+sSHgsVwLbNEjcnkuhZTuhsVimlRME6lhCxKEYAComCA0LyYkD+TAJCpKb0mQeECE7d9yqeC9L3ZAvRhAmcn8JHq6kzQh3T1kuzh0vG14pkUy49g0uNCb3ZRUlERN7EOzgiIp/wW5AJGzgiIp9Qyp3nZ25OwTeQ2EVJRESexDs4IiK/cKmLUhckNJSwgTMhRUwlkvNOeE8M8BKjJaUIR02UmmFkohSVqMtFKe1DzDnp0jHp1pH2YbodHcsyWycYdD6ndH+fFGFp+vdZlvAdCdGYuuOy0swic8XoSkCO8hNzpTovbrlYL41zUQ5RSgUSig7uvZ3U6KNkFyUREXkS7+CIiPzCttzpXmQXJRERDSV+y2TCLkoiIvIk3sEREfkEB3r7hKWN4hqEiCnDXJRSfkB5lmpNlKHpjN7SDNmabgox+lHMUWkWqSnlV+x5z72/w5SYc1KITJQiH6U8kYmwAs4nVcDF79X8nErk/BQiOMXo5kGI9BP2rb2+JBGjKImIiDzAt3dwRER+wy5KIiLyJEZREhEReQDv4IiIfMJvd3Bs4AxIue3kKC5dzjvTcsNcf93meSLl/JFmy/esI0VFmuWolKIlownkojSN7HSTNNu2FF3pZpSaac5JaXkAsKNCRKZQrqTyBM5POR+rsCGxXK6XUl0W81qmGKVcegaXIg0cuyiJiMiTeAdHROQTfhsHxwaOiMgn/DZMgF2URETkSbyDIyLyCUZREhGRJ7GB8xo1CElPxWnuNaHF0kNaqW9bOqHE5c3D6OXQfiHRs+Ykl5IIm25LGg4gJWcG5LD/wUi2LFGGwwSkYQU6luV8vkmfeVQYDhAIyuet6TkiJXTWnZ/u1QGpjumG7wzCcIDBuCYRAD6DIyLyDWV/EWjSv5f5vjds2ICCggJkZmaiuLgYe/bs0S5/6tQpPPTQQ7jiiisQCoVw4403Yvv27Ub79P4dHBERAUheF+XWrVtRWVmJ2tpaFBcXo6amBmVlZTh48CBGjx7da/lIJIK//Mu/xOjRo/HLX/4SY8eOxSeffIIRI0YY7ZcNHBERJaS9vT3u/6FQCKFQqNdy69evx+LFi1FeXg4AqK2txbZt21BXV4cVK1b0Wr6urg6ff/45du3ahfT0dABAQUGB8fGxi5KIyCcuDPR24wUA+fn5yMnJib2qq6t77TMSiaCpqQmlpaWxskAggNLSUjQ2Njoe57//+7+jpKQEDz30EHJzczF+/HisW7cO0ajZZNS8gyMi8glbWa7kXr2wjWPHjiE7OztW7nT3dvLkSUSjUeTm5saV5+bm4sCBA47bP3z4MH7zm99g3rx52L59Ow4dOoQHH3wQXV1dqKqq6vNx+reBkxIkQ5c82TCpsq2JhBNWESP9pGTLQrmtS2bbLSRPliITpX1oEh5LUY5SuRh1KS2viaKUoyWlYxqEZMsBKYpSSEacQBSl9BlKLOH81H2vgYB0HgrRscK5pjs/A4bnuvR9S3+fvl5KdVlKwmx4rfCY7OzsuAbOLbZtY/To0XjppZcQDAZRVFSE48eP45lnnmEDR0REDlxK1SUOzXAwcuRIBINBtLa2xpW3trYiLy/PcZ0rrrgC6enpCAa/+IF0yy23IBwOIxKJICMjo0/7HpBncMkIByUiIr0LUZRuvPoqIyMDRUVFqK+vj5XZto36+nqUlJQ4rnP77bfj0KFDsC+6o/7oo49wxRVX9LlxAwaggbsQDlpVVYV9+/Zh4sSJKCsrw2effea4/IVw0KNHj+KXv/wlDh48iI0bN2Ls2LFuHxoRESVBZWUlNm7ciJ/97Gf48MMP8cADD6CjoyMWVblgwQKsXLkytvwDDzyAzz//HEuWLMFHH32Ebdu2Yd26dXjooYeM9ut6F2WywkGJiEgvWePgZs+ejRMnTmDNmjUIh8MoLCzEjh07YoEnLS0tCAS+uN/Kz8/HW2+9hWXLluHWW2/F2LFjsWTJEixfvtxov642cBfCQS9uiU3CQd98802MGjUKc+fOxfLly+P6Xy/W2dmJzs7O2P+/PBaDiIh6S2YuyoqKClRUVDi+19DQ0KuspKQEv/3tb433czFXG7jBCgetrq7G2rVr3Tz0vjHNU6fNeSeUdwuRYlLkoxQ5pzkBjXNR2s770EXtifkPXYqWlI6p5z13clEmEk4t5ZCUJoi0LOd9SFGXOlJeS9PlowH5s5XyVFpRIRpUzEWpOT+Fc0SqA1KdEetYIvkmByNHJbku6QO9Lw4HLSoqwuzZs/HYY4+htrZWXGflypVoa2uLvY4dOzaIR0xElJpsFXDtlQpcvYMbrHBQKR0MERHJlHJpRu8UmS7H1WY4meGgREREF3P9PjNZ4aBERKSXjHFwyeT6MIFkhYMSEZEeZ/R2QTLCQV1lmkcugVyU0oSBYj5IwxyVdgIRjspwtmbtrNpClKNb0ZLSTN+AZnZwzfE6bgfmlViKlpQEA1L0YQK5KA2jKOXvyHxG74t/tPZled35KeWilOuAcD6L+SMTyEUpLu+PnJOpirkoiYh8wu3ZBIY6NnBERD7hty7K1BjMQEREZIh3cEREPuG3Ozg2cEREPsFncD5h2d0JrCNFSxqWA5o8eYaRYlLePk22AimCTZrxWswfqYlklKIcpUjGaNT5VDTdDgAxjZBUKQfj16iU91H6ngLCTN96UnV2PtetqPMxSRGRgJynMpjmHE0onmua81M+p81mnE8oF6Xwnlj3Bbrri3l8LCXKtw0cEZHfKOXOD7oERrEkBRs4IiKf8NszOEZREhGRJ/EOjojIJ5RLQSapcgfHBo6IyCfYRUlEROQBRndw06dPR2FhIWpqagbocIYGS0qgKpVLmZO1SV2FcimUvVtKKCssr0lmKydoNkuqLCXr1a8jDRMwTMKsmVFYXsdsmEAiXTkBYTiAOExASpCs++lpOILAsoTPQ9i3/nt13rmYhNnF81OqA1KdkYcJaOqlWJelxM2plWzZb3dw7KIkIvIJvw307nMX5T333IP//M//xPPPPw/LsmBZFo4ePTqAh0ZERJS4Pt/BPf/88/joo48wfvx4PPHEEwCAUaNGDdiBERGRu9hFKcjJyUFGRgaysrKQl5c3kMdEREQDgF2UREREHuD5IBOlhESzupWMp62XIq/kVZQQdaa6zRIhS9FoukTI0r7F6ENpeV2yZTFa0vl4peW7hYg+XbJl6T2pW8XN7hbpK5eiKKVy3S/kYMDseC3hM5f2HdB8r1IiZukzD0qRudrz07AOSHVG2IeljW52qe5rSNekwaBgQemvfn3eTiowauAyMjIQjaZWWCwREfXw2zM4oy7KgoIC7N69G0ePHsXJkyfFMTFERETJZtTAff/730cwGMS4ceMwatQotLS0DNRxERGRyy4EmbjxSgVGXZQ33ngjGhsbB+pYiIhoALGLkoiIyAM8H0UpSiSHnBgtKZTrpr2V8u0Z5uETIwM1UYZSrkFbisiUlk9gH93CPkyjJXVRlGIeTCHya1ByUUKIWEzi1MiBgJSLUn627ta5ozs/xXPaMB+rnO81gShK4+jKoRmMZ8OlcXApEkXJOzgiIvIk/97BERH5jN+ewbGBIyLyCRuWK92L7KIkIiJKIt7BERH5hUtdlOIks0MMGzgH8ozehhFWusArw5mOxXIhD58u1584c7dhzkkpIhIAuqXjEvIJirkopdyVCczoLQXPuTloNSp03cgzfTtvJ6jpAjK9QEn7lr4jy5JP3EDA+T3TnJPaXJRibkmzcqmO6eulWR1PtRm9OZsAERGRB/AOjojIJxhFSUREnmRD30Nrsp1UwC5KIiLyJN7BERH5BLso/cI2n1VXjq6UwvPkbUmRYqY5KqVoNGk2Y0Azq3bU+XQQIxk1UZRSJGOXlIvSMFqyW/f3CetIkV+DUVnF2bOFcm2KSuljNwzoE/Nj6qJjxShK53NHmiA5LYEoXznnpJS7Uqpjmg9KqMvG0ZIJXF8Gg63ciYDUTYo+lLCLkoiIPMm/d3BERD6jYEG5kGbLjW0MBjZwREQ+wYHeREREHsA7OCIin+gJMnFnO6mADRwRkU/wGZzXJBKu61KyZaXbtdCHbRsmlJWWt2051NsWQvilpMrS8tpky8L+xaTK4vLOn1OX5u+Twv6j0mc+CM8TpOEAQaHclrIwa6iA8zrSEIVuS0i2HJV/ngeDznXD9NzRnp8u1YGgNCxEUy8t44TqCSRbHqJDCLzI+w0cEREB8F+QCRs4IiKfUOoSSQQMtpMKGEVJRESexDs4IiKfULBgM8iEiIi8hsmWfUKbPNU4ijKBZMuGiWOlBLRSebeUaBZAVHhPjq4UlheSGgMDHy3ZLWy/57ikKEqzJMxukqIobcv5JJGW13M+by0hWlKM7AyYJ+q2pehK4VzTnZ9Bw3PdNEG5djIzsS6bRVEaJ2emAeHbBo6IyG/8FkU5IEEmGzZsQEFBATIzM1FcXIw9e/b0ab0tW7bAsizMmjVrIA6LiMjXlIuvVOB6A7d161ZUVlaiqqoK+/btw8SJE1FWVobPPvtMu97Ro0fx/e9/H1//+tfdPiQiIvIh1xu49evXY/HixSgvL8e4ceNQW1uLrKws1NXVietEo1HMmzcPa9euxbXXXuv2IREREb7oonTjlQpcbeAikQiamppQWlr6xQ4CAZSWlqKxsVFc74knnsDo0aNx77339mk/nZ2daG9vj3sREZGe7eIrFbgaZHLy5ElEo1Hk5ubGlefm5uLAgQOO67z77rv4l3/5FzQ3N/d5P9XV1Vi7dm1/DlXLPB+dVC7/ylFCFJmUb08ulyLONFGUwntSbkmpvEuXi1LclvPxmkZL6nJRShGZ0vgfKeQ5kUos/WKU8kFGhWMKBgb+KUcAQi5KTQRnQIj6DEr5IA3zRwLyOW1cN6RITU29NK3j4rWChoSkZjI5ffo05s+fj40bN2LkyJF9Xm/lypVoa2uLvY4dOzaAR0lE5A0XxsG58UoFrt7BjRw5EsFgEK2trXHlra2tyMvL67X8xx9/jKNHj2LmzJmxMvtPv4jS0tJw8OBBXHfddb3WC4VCCIVCbh46EZHncZhAP2RkZKCoqAj19fWxMtu2UV9fj5KSkl7L33zzzXjvvffQ3Nwce/3VX/0V7rzzTjQ3NyM/P9/NwyMiIh9xvYuysrISGzduxM9+9jN8+OGHeOCBB9DR0YHy8nIAwIIFC7By5UoAQGZmJsaPHx/3GjFiBL7yla9g/PjxyMjIcPvwiIh8K5nj4JIxPtr1TCazZ8/GiRMnsGbNGoTDYRQWFmLHjh2xwJOWlhYENKmAiIhoYCSri/LC+Oja2loUFxejpqYGZWVlOHjwIEaPHi2u19/x0QOSqquiogIVFRWO7zU0NGjX3bRpk/sH5EQ73bZAipjqdv49o6KaKEopuqzbLA+fNDOymLcPchSlGF0pRThqIuGkKEepPCJFakq5KzUVTHpPjJYchChKKe+jFLFoSzkRoTleZZb/UJo0XDujd8D5U0kzPKd056fpOS3VGTnqUj53LKEui3Vfksj1xcMuHh8NALW1tdi2bRvq6uqwYsUKx3UuHh/9X//1Xzh16pTxfnkrRUTkE26Pg/vyeOTOzs5e+xys8dFO2MAREfmE28ME8vPzkZOTE3tVV1f32qdufHQ4HHY8zgvjozdu3Nivv5ezCRARUUKOHTuG7Ozs2P/dGL6V6PhoJ2zgiIh8QsGdNFsXnlRmZ2fHNXBOBmt8tBN2URIR+YSCS12UQoo5J8kcH+39Ozi7K4F1hGg0Jfz2UVLkleb3gzQLsWHuPmlm5GhU/mrF2bZNc1Rq8kFKOSRNc0tKy+uiKLvFGb3NoivdJEVLBoVypckHqXvPiRTBKZfrZtt2rgPdAed10sRZ4uXzUzqnpUhNqc6IOSd1uSjFOm42o7dWItekFFdZWYmFCxdi8uTJmDJlCmpqanqNjx47diyqq6tj46MvNmLECADoVX4p3m/giIgIAGCrnpcb2zGRrPHRbOCIiHzCrdm4E9lGMsZH8xkcERF5Eu/giIh8wm+zCbCBIyLyCbdm406VaV7ZRUlERJ7k2zs4SxfeazhtvfRzRgnhzoAc9i8mlBXLpWEC8m+X7m7nr900qXKXZh9S8uROITxcHj4glOuSLRsOLRiMX6PSJyUNE4gKSY0BIN1wmIBpNbc0IQRBYQhBd1AaxuK872ianIx4oOuGrl6KYf+G1wTt9SWJ3JqN25czehMR0dDFLkoiIiIP4B0cEZFPKCUnXjLdTipgA0dE5BM2LNgGeSR120kF7KIkIiJP4h2cA8swYkqanV5pogxVt/N7YrJlKRpNiAiTIh9174nRkmIiZHkfpsmTO42TM8ufrRQtKUZRDkJ3S0D4wStFUQY1UWpKiLA0ffAvJVWWjgkAuoQPS0qELCVh1p2f0jkt1QGpzsh1TFcvncula4J4rRiikpWLMlnYwBER+YVLz+BcSWg5CNhFSUREnsQ7OCIin/BbkAkbOCIin/DbMAF2URIRkSd5/w5OCovS5qI0zUcnbUcTreVWvj0hWjEq5HwEzHNOSpGMUr5J3XtStKS0vJRzUjomAOgWfl1KUZSD8WvUEqMoncvTpBWgOd1cyg8Y0OWiFCI4xXNHyFGpOz+lc9q0bkh1TFcv5Q9XqvvStUJzfZGuSYPAb6m6vN/AERERAP8NE2AXJREReRLv4IiIfELBnSFsKXIDxwaOiMgverooXRgmkCItHLsoiYjIk3x7B5fQjN7dwjrdQnSeLued8F60S5jxWpqFW8pFqZm1WIp4kyIZ5YhIeR+m0ZIRMeek82erndFbzEXpvLz0izaR6Eop+DEg5ZwUVohq8kHaSogaDLjzs9oSclQCQDDqXDeClnO5mKNSd36K57SQj1WoM1Id09VLqS6LdT/lZvT21zg43zZwRER+47dhAuyiJCIiT+IdHBGRT7CLkoiIPIldlERERB7AOzgnYt455/tyJUXhaSLFTCO/TGfu1s22HRHyAEqzcHcKy+vyQZpGS3YK0ZLS8lJEJCBHXkrfk5S70k1SbklLiJZM1/z0jArbssXf1Wa/Y9M0EZzSuZAm5KiUzrUMW87HKJ3TUh0wjUjW1UvpHLGkgV8pNqO3cilVF7soiYhoSPFbJhN2URIRkSfxDo6IyCf8NpsAGzgiIp/w2zABdlESEZEn8Q6OiMgn/DYOzvsNnBSOrEmGKiZKFZMwC2HbQrgzACgxFNr5K7GFEH4pAW2XZt/dQli1aWj/ec0+5HXMhgl0J5RsWSiXRn9A+P4S6IYJSLl6hbizoLCClDAaANKFMH4xCbNhf1JA872KwxqEddIDznVJe34a1gGpzkh1TFcvxWTLpkmVtcnc5SESA81vz+DYRUlERJ7k/Ts4IiIC4L9xcGzgiIh8gl2UREREHsA7OCIin/DbODjfNnDaKeXFaEkh2bIUxSUkewXkRLBSBFm3kCC2S4qiFMoBOQGuGEUplHdqki1LUZFyUmWzcimhcs97zuVSnewehJjnNCFa0hL2rU22LLyXYXhMlpgAWt55UIqWFP6Q9ICQbLlbrn9dQbM6INUZMdmypl6KiZi7u5zLTaMrk8xvwwTYRUlERJ40IA3chg0bUFBQgMzMTBQXF2PPnj3ishs3bsTXv/51fPWrX8VXv/pVlJaWapcnIqLE2Pgi0KRfr2T/IX3kegO3detWVFZWoqqqCvv27cPEiRNRVlaGzz77zHH5hoYGzJkzB2+//TYaGxuRn5+Pu+66C8ePH3f70IiIfE25+EoFrjdw69evx+LFi1FeXo5x48ahtrYWWVlZqKurc1z+lVdewYMPPojCwkLcfPPNePnll2HbNurr690+NCIi8hFXG7hIJIKmpiaUlpZ+sYNAAKWlpWhsbOzTNs6ePYuuri5cfvnl4jKdnZ1ob2+PexERkZ5KpDvS4eXLKMqTJ08iGo0iNzc3rjw3NxcHDhzo0zaWL1+OMWPGxDWSX1ZdXY21a9f261i1U813C7nihCSHSghrs7vl3w+2sE7UMFqyW8jPJ0U+6t6LSNsyzCupe69TzEXpvB05ilLctfielKNSqqyJPGeQPpGIsLGgEAwqHWvPe2a5M5Umr6WpNCkXpXDupAsJQHXnZ0jYllQHpDoj1TFdvZTqsviFSNcK3fUliZRyKZNJijRwQyqK8qmnnsKWLVvw+uuvIzMzU1xu5cqVaGtri72OHTs2iEdJRESpwNU7uJEjRyIYDKK1tTWuvLW1FXl5edp1n332WTz11FP49a9/jVtvvVW7bCgUQigU6vfxEhH5CcfB9UNGRgaKioriAkQuBIyUlJSI6z399NN48sknsWPHDkyePNnNQyIioj/peYamXHgl+y/pG9czmVRWVmLhwoWYPHkypkyZgpqaGnR0dKC8vBwAsGDBAowdOxbV1dUAgB/+8IdYs2YNNm/ejIKCAoTDYQDA8OHDMXz4cLcPj4iIfML1Bm727Nk4ceIE1qxZg3A4jMLCQuzYsSMWeNLS0oJA4IsbxxdffBGRSATf/e5347ZTVVWFxx9/3O3DIyLyLU6X44KKigpUVFQ4vtfQ0BD3/6NHjw7EIXxBmD03oVyUUo5DIYpLmlG457CEiDAhT6SUby8ilQvbAYBOIUpNmqFbKpeiKwE5WvK88LGfl3JXCsvr8kcKKUMRFUK/dBGLbpGiJYNCPkjd39clBfoJO5Fn9HbekDQrOSB/f0EhujIt4PyHZGjOT+mcFnNOCtuS6piuXoq5KKXvI5FclMme0dul7aSCIRVFSURE5BbfziZAROQ36k//3NhOKmADR0TkE+yiJCIi8gDewRER+YTfBnr7t4HTRTkp569PCcFPUlSWrZvRW5ydWMjDJ87CbTY7N2AeLXk+6hxWd06T0++ssI4YRSmUR4S+EF0uStNoye5BSKyXJkRLStGHUnQlAKSLUaJCjkohulL+COXvNSA8ewlazueOFEUZCuhypQqz1IsRxs7bEnNU6mb0FuqAVPct4Vqhvb4kkVIuPYNLkWSU7KIkIiJP8u8dHBGRz7CLkoiIPIldlERERB7AOzgiIp9QcKd7MTXu39jAERH5hq0UbBeaJztFuii938BJ8b26MF5pGvpuIdxaCDuWwpR7diGEQhsnVZZC+zXJlsV1pATJZomTde+dEz7a80IMf6c0TEBTwbqEBLhRoWJLlTWRZxUWnM+RgDRMQFg+PSA/PUgXElMLXxOUMHxAontuEbSc35WSSacLwwHOB+TzMxR1PkmkOiDVGTE5s6ZeikMIhLovXiu0w5CSl2w5mTZs2IBnnnkG4XAYEydOxI9//GNMmTLFcdmNGzfi5z//Od5//30AQFFREdatWycuL+EzOCIin1Au/jOxdetWVFZWoqqqCvv27cPEiRNRVlaGzz77zHH5hoYGzJkzB2+//TYaGxuRn5+Pu+66C8ePHzfaLxs4IiKfsF18mVi/fj0WL16M8vJyjBs3DrW1tcjKykJdXZ3j8q+88goefPBBFBYW4uabb8bLL78M27ZRX19vtF82cERElJD29va4V2dnZ69lIpEImpqaUFpaGisLBAIoLS1FY2Njn/Zz9uxZdHV14fLLLzc6PjZwREQ+YUO59gKA/Px85OTkxF7V1dW99nny5ElEo1Hk5ubGlefm5iIcDvfpuJcvX44xY8bENZJ94f0gEyIiAuB+FOWxY8eQnZ0dKw+FQv3e9pc99dRT2LJlCxoaGpCZmWm0rm8bOEuItAMgTkOvhATCSkg6HNUkde0W3uvqSncslxLQdgqRYlLiZN1754XoPClxslQOyNGSZ4VoyfNR58+8U/guIkqOUusSnhB0w3kdqcK7GkUpRDKmwfm7SJdCIgFkCImNo8p5HVvsqJG+P/l7lXJABwNCNKhwrmUISZgBIFOMGDarM1Id09VLqS5LdV+6jmivLx6SnZ0d18A5GTlyJILBIFpbW+PKW1tbkZeXp1332WefxVNPPYVf//rXuPXWW42Pj12UREQ+kYwoyoyMDBQVFcUFiFwIGCkpKRHXe/rpp/Hkk09ix44dmDx5ckJ/r2/v4IiI/Obi52f93Y6JyspKLFy4EJMnT8aUKVNQU1ODjo4OlJeXAwAWLFiAsWPHxp7h/fCHP8SaNWuwefNmFBQUxJ7VDR8+HMOHD+/zftnAERHRgJo9ezZOnDiBNWvWIBwOo7CwEDt27IgFnrS0tCBwUXKDF198EZFIBN/97nfjtlNVVYXHH3+8z/tlA0dE5BPJuoMDgIqKClRUVDi+19DQEPf/o0ePJnBUvbGBIyLyiUSykEjbSQWeb+AsO5FclM7vKSGHnZS/zhaiwQA5T54UKRYR9t1pS3kldVGUQlSkEEF2TlheipQEgDPdzhXgXNT5sz0nfB8RIfKxE13iviOW83tROB9wtzXwuQHTlPP3GhSqYIZyjgwEgJDwXpfwnUeVcN4K5booSjGnppCjMt1yPg8yhByVgHxOS3VAqjNSHdPVS6kuS3Vfulbori/iNYlc5/kGjoiIeiiXuih5B0dEREOKbdmwrP6P0bNdmVVu4HEcHBEReRLv4IiIfMKGgpWkKMpkYANHROQTF1Ilu7GdVOD9Bk6KZtLlipOiKLsynDclzTSsyXlnOnN3pxApdk7MRSn3Pp8V3pOiJTuEgEUpUhIAOoRoybNCBNk5RBzLz1vnHcs7hXIAiFjO24oKkZdRJUdkuiUYcI58DMK5PEM5n2sAEFLOCWczhfKoLZy34h7kCEcpwjIoRFdKM5NnBORzRzqns9LM6oxUx3T1UqrLqkuKonQ+17TXF10EN7nK+w0cEREB6PlR404XZWpgA0dE5BOMoiQiIvIA3sEREfmEDRuWC3dfqXIHxwaOiMgn2MD5hNWtiZwTogNNc1FGNTnvpFmIO7udy88bztwtRUoCQIeQc7JDSJEnRUue7pZz6nUIkYkdQvTjOeusY/l5q8OxPKLOifvuUs776LY7HculkGelmTVcYgmzbVvC04C0QMixPB3OEZEA0BkY5lgeUZc5lnepLMfyqO28D1uTB1O6ZEgzeks5KnVRlJnSjPNShLFQZ6Q6pquX5rkonf8O7fWFBo1vGzgiIr/hODgiIvIkRlESERF5AO/giIh8QsF25e6LXZRERDSkKEShXOi4U8JExEMNuyiJiMiTPH8HJ04Pr0uGGhWGCQghxN0R53DkbiFMGZBDmzuFEOazQrJlaTiANBQAAM4KH0lHl/PfLSVOloYCAEC7dcZ5HaH8HNodyztt5+W7bOdhBYA8HCBqOyfGtQch2XLAEpItB5wTIUvDBwAgPeAc9h+xnIdOdFnZjuVRDHfegRLKAQSEhNxpQjLiNCEJc4aQhBkAMoPO7w0T6sBlhsMHdPVSqstS3ZeuFbrri3hNGgQ93ZP+CTLxfANHREQ9euZxc6OBS4354NhFSUREnsQ7OCIin+gJMnHuNjbdTipgA0dE5BN+ewbHLkoiIvIk79/BCdPD65KhKinwUoiksoUIsk4hIgsAOoVIrvNCpNg5IVKsQ0gC29Etd0PIyZOdP6vTQvThKcs58hEA2gP/z7H8rO1cHhGiJTu7TzuWR2052bKthKTKYrTkYHS3CEmYhejKgCVHUQaFZMtdac6RpZGAc3lX4KuO5bYuwth2jsgMdDtHgwYDzn93upCcGQAyg87rZAWlhONCFKWU0FxTL6W6LCdhdt6ONtmycE0aDMxFSUREnmQjCrjwDM5OkWdwA9JFuWHDBhQUFCAzMxPFxcXYs2ePdvlf/OIXuPnmm5GZmYkJEyZg+/btA3FYRETkI643cFu3bkVlZSWqqqqwb98+TJw4EWVlZfjss88cl9+1axfmzJmDe++9F/v378esWbMwa9YsvP/++24fGhGRr13oonTjlQpcb+DWr1+PxYsXo7y8HOPGjUNtbS2ysrJQV1fnuPzzzz+Pb33rW3j44Ydxyy234Mknn8Rtt92Gn/zkJ24fGhGRr9kq6torFbj6DC4SiaCpqQkrV66MlQUCAZSWlqKxsdFxncbGRlRWVsaVlZWV4Y033hD309nZic7OLwIJ2traAADt7b0DD86ddn7YG+mQf4Goc85BGJ2dzuuciTh/2Wd0M14LD6HPRp0DOs4Jsxx32s7HGrHlfnYhIxe6hJO2WzkfUxTyg3RbePouVQxp9mylpNm25UwK8num5W4S0r8Jx6o0xyR/Js6fofSZS99RVJO6TDoX5HPK+byNaG4ApHP6XNSszkh1LFNTL0NCXVZC3Q8J1wpLc305L1yTutrjg4EuXM905zrpudrAnTx5EtFoFLm5uXHlubm5OHDggOM64XDYcflwOCzup7q6GmvXru1VXnDVkgSO2kSrYTnRBdIFz/liZys516YUhNdl+KPaOZ7VZfKf4XF/SOC91x1L//jHPyInJ6ffRwQwijIlrFy5Mu6u79SpU7j66qvR0tLi2ongB+3t7cjPz8exY8eQne0c/k3x+Jklhp+buba2Nlx11VW4/PLLXdtmTwPX/+5FXzZwI0eORDAYRGtr/B1Na2sr8vLyHNfJy8szWh4AQqEQQqHe44RycnJYeRKQnZ3Nz80QP7PE8HMzF9DMvEB6rn5yGRkZKCoqQn19fazMtm3U19ejpKTEcZ2SkpK45QFg586d4vJERJQYpWzYLryk58BDjetdlJWVlVi4cCEmT56MKVOmoKamBh0dHSgvLwcALFiwAGPHjkV1dTUAYMmSJZg2bRqee+45zJgxA1u2bMHevXvx0ksvuX1oRES+1tO16EayZZ82cLNnz8aJEyewZs0ahMNhFBYWYseOHbFAkpaWlrhb7qlTp2Lz5s1YtWoVHn30Udxwww144403MH78+D7vMxQKoaqqyrHbkmT83MzxM0sMPzdz/Mz6z1KMQSUi8rT29nbk5OQgJ3McLMt56IYJpaJoO/97tLW1DelnqikZRUlEROZs2LB81EXJ8BwiIvIk3sEREflET/SjC3dwfo2iJCKiocmNQd5ubmegsYuSiIg8KWUbuB/84AeYOnUqsrKyMGLEiD6to5TCmjVrcMUVV2DYsGEoLS3F//7v/w7sgQ4hn3/+OebNm4fs7GyMGDEC9957L86ccZ5J+4Lp06fDsqy41/333z9IR5wcnM/QnMlntmnTpl7nVGZm5iAe7dDwzjvvYObMmRgzZgwsy9ImmL+goaEBt912G0KhEK6//nps2rTJaJ9KKag/DdTu3ys1gu9TtoGLRCL467/+azzwwAN9Xufpp5/Gj370I9TW1mL37t247LLLUFZWhvPnzw/gkQ4d8+bNwwcffICdO3fiV7/6Fd555x3cd999l1xv8eLF+PTTT2Ovp59+ehCONjk4n6E5088M6EnZdfE59cknnwziEQ8NHR0dmDhxIjZs2NCn5Y8cOYIZM2bgzjvvRHNzM5YuXYpFixbhrbfe6vM+/TYfHFSK++lPf6pycnIuuZxt2yovL08988wzsbJTp06pUCik/u3f/m0Aj3Bo+P3vf68AqP/5n/+Jlf3Hf/yHsixLHT9+XFxv2rRpasmSJYNwhEPDlClT1EMPPRT7fzQaVWPGjFHV1dWOy//N3/yNmjFjRlxZcXGx+vu///sBPc6hxPQz62ud9RMA6vXXX9cu88gjj6ivfe1rcWWzZ89WZWVll9x+W1ubAqCGZRSorNC1/X4NyyhQAFRbW1t//uwBl7J3cKaOHDmCcDiM0tLSWFlOTg6Ki4vFueq8pLGxESNGjMDkyZNjZaWlpQgEAti9e7d23VdeeQUjR47E+PHjsXLlSpw96805UC7MZ3jxOdKX+QwvXh7omc/QD+cUkNhnBgBnzpzB1Vdfjfz8fHznO9/BBx98MBiHm9LcONeUirr2SgW+iaK8ML+c6dxzXhEOhzF69Oi4srS0NFx++eXav3/u3Lm4+uqrMWbMGPzud7/D8uXLcfDgQbz22msDfciDbrDmM/SSRD6zm266CXV1dbj11lvR1taGZ599FlOnTsUHH3yAK6+8cjAOOyVJ51p7ezvOnTuHYcOGXXIbboX3p8owgSF1B7dixYpeD5+//JIqjV8N9Gd23333oaysDBMmTMC8efPw85//HK+//jo+/vhjF/8K8pOSkhIsWLAAhYWFmDZtGl577TWMGjUK//zP/5zsQyOPGVJ3cP/4j/+Ie+65R7vMtddem9C2L8wv19raiiuuuCJW3traisLCwoS2ORT09TPLy8vr9dC/u7sbn3/+uXbuvS8rLi4GABw6dAjXXXed8fEOZYM1n6GXJPKZfVl6ejomTZqEQ4cODcQheoZ0rmVnZ/fp7g1wL8VWqgSZDKkGbtSoURg1atSAbPuaa65BXl4e6uvrYw1ae3s7du/ebRSJOdT09TMrKSnBqVOn0NTUhKKiIgDAb37zG9i2HWu0+qK5uRkA4n4keMXF8xnOmjULwBfzGVZUVDiuc2E+w6VLl8bK/DSfYSKf2ZdFo1G89957uPvuuwfwSFNfSUlJryEopuea37ooUzaK8pNPPlH79+9Xa9euVcOHD1f79+9X+/fvV6dPn44tc9NNN6nXXnst9v+nnnpKjRgxQr355pvqd7/7nfrOd76jrrnmGnXu3Llk/AmD7lvf+paaNGmS2r17t3r33XfVDTfcoObMmRN7/w9/+IO66aab1O7du5VSSh06dEg98cQTau/everIkSPqzTffVNdee636xje+kaw/YcBt2bJFhUIhtWnTJvX73/9e3XfffWrEiBEqHA4rpZSaP3++WrFiRWz5//7v/1ZpaWnq2WefVR9++KGqqqpS6enp6r333kvWnzDoTD+ztWvXqrfeekt9/PHHqqmpSX3ve99TmZmZ6oMPPkjWn5AUp0+fjl23AKj169er/fv3q08++UQppdSKFSvU/PnzY8sfPnxYZWVlqYcfflh9+OGHasOGDSoYDKodO3Zccl8XoijTg7kqI+2Kfr/Sg7kpEUWZsg3cwoULFYBer7fffju2DAD105/+NPZ/27bV6tWrVW5urgqFQuqb3/ymOnjw4OAffJL88Y9/VHPmzFHDhw9X2dnZqry8PO4HwZEjR+I+w5aWFvWNb3xDXX755SoUCqnrr79ePfzww0P+pO6vH//4x+qqq65SGRkZasqUKeq3v/1t7L1p06aphQsXxi3/6quvqhtvvFFlZGSor33ta2rbtm2DfMTJZ/KZLV26NLZsbm6uuvvuu9W+ffuScNTJ9fbbbztewy58VgsXLlTTpk3rtU5hYaHKyMhQ1157bdz1TedCA5cWHKXS03L7/UoLjkqJBo7zwRERedyF+eCCgcthWf2PLVTKRtT+fMjPBzekoiiJiIjcMqSCTIiIaCApwJUIyNTo+GMDR0TkE+7NB5caDRy7KImIyJN4B0dE5BM9A7RduINjFyUREQ0t7jRwqfIMjl2URETkSbyDIyLyC5eCTJAiQSZs4IiIfMJvz+DYRUlERJ7EBo7oEk6cOIG8vDysW7cuVrZr1y5kZGSgvr4+iUdGZMp28TX0sYEjuoRRo0ahrq4Ojz/+OPbu3YvTp09j/vz5qKiowDe/+c1kHx6RAdXz/Ky/rwS6KDds2ICCggJkZmaiuLgYe/bs0S7/i1/8AjfffDMyMzMxYcKEXlMF9QUbOKI+uPvuu7F48WLMmzcP999/Py677DJUV1cn+7CIUsLWrVtRWVmJqqoq7Nu3DxMnTkRZWVmvSZgv2LVrF+bMmYN7770X+/fvx6xZszBr1iy8//77RvvlbAJEfXTu3DmMHz8ex44dQ1NTEyZMmJDsQyLqkwuzCQBBuDcOLtrn2QSKi4vx53/+5/jJT34CoGdS3Pz8fPzDP/wDVqxY0Wv52bNno6OjA7/61a9iZX/xF3+BwsJC1NbW9vkoeQdH1Ecff/wx/u///g+2bePo0aPJPhyiBDlOQ2f46tHe3h736uzs7LW3SCSCpqYmlJaWxsoCgQBKS0vR2NjoeISNjY1xywNAWVmZuLyEDRxRH0QiEfzt3/4tZs+ejSeffBKLFi0Su1eIhpqMjAzk5eUBiLr2Gj58OPLz85GTkxN7OXXbnzx5EtFoFLm5uXHlubm5CIfDjscbDoeNlpdwHBxRHzz22GNoa2vDj370IwwfPhzbt2/H3/3d38V1oRANVZmZmThy5AgikYhr21RKwbLiuztDoZBr23cDGziiS2hoaEBNTQ3efvvt2POGf/3Xf8XEiRPx4osv4oEHHkjyERJdWmZmJjIzMwd9vyNHjkQwGERra2tceWtr65/uKnvLy8szWl7CLkqiS5g+fTq6urpwxx13xMoKCgrQ1tbGxo3oEjIyMlBUVBQ3ZtS2bdTX16OkpMRxnZKSkl5jTHfu3CkuL+EdHBERDajKykosXLgQkydPxpQpU1BTU4OOjg6Ul5cDABYsWICxY8fGnuEtWbIE06ZNw3PPPYcZM2Zgy5Yt2Lt3L1566SWj/bKBIyKiATV79mycOHECa9asQTgcRmFhIXbs2BELJGlpaUEg8EWH4tSpU7F582asWrUKjz76KG644Qa88cYbGD9+vNF+OQ6OiIg8ic/giIjIk9jAERGRJ7GBIyIiT2IDR0REnsQGjoiIPIkNHBEReRIbOCIi8iQ2cERE5Els4IiIyJPYwBERkSexgSMiIk/6/5sex0hmMCxEAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(5, 5))\n", "plt.imshow(\n", " sol.ys.vals,\n", " origin=\"lower\",\n", " extent=(x0, x_final, t0, t_final),\n", " aspect=(x_final - x0) / (t_final - t0),\n", " cmap=\"inferno\",\n", ")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"t\", rotation=0)\n", "plt.clim(0, 1)\n", "plt.colorbar()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b4b4ced9-0602-4354-a1b9-277ddf70245c", "metadata": {}, "source": [ "Some final notes.\n", "\n", "1. We wrote down the general Crank–Nicolson method, which uses a fixed point iteration to solve the implicit problem. If you know something about the structure of your problem (e.g. that it is linear) then it is often possible to more specialised solvers, which run faster. (E.g. linear solvers.)\n", "\n", "2. To keep this example brief, we didn't worry about doing a von Neumann stability analysis." ] } ], "metadata": { "kernelspec": { "display_name": "py38", "language": "python", "name": "py38" }, "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.16" } }, "nbformat": 4, "nbformat_minor": 5 }